*$ CREATE MAGFLD.FOR
*COPY MAGFLD
*
*===magfld=============================================================*
*
      SUBROUTINE MAGFLD ( X, Y, Z, BTX, BTY, BTZ, B, NREG, IDISC )

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1988-2010      by Alberto Fasso` & Alfredo Ferrari *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *
*     Created  in     1988         by    Alberto Fasso`                *
*                                                                      *
*                                                                      *
*     Last change on 06-Nov-10     by    Alfredo Ferrari               *
*                                                                      *
*     Input variables:                                                 *
*            x,y,z = current position                                  *
*            nreg  = current region                                    *
*     Output variables:                                                *
*            btx,bty,btz = cosines of the magn. field vector           *
*            B = magnetic field intensity (Tesla)                      *
*            idisc = set to 1 if the particle has to be discarded      *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(RTDFCM)'
      INCLUDE '(LTCLCM)'


* Ovde Andrija pocinje njegovo magnetno polje
* Ovde je definisana varijabla BFIRST za prvo skaniranje polja
      LOGICAL     BFIRST
      DATA        BFIRST / .TRUE. /
      SAVE        BFIRST
* Podaci magnetnog polja (T)
      DATA DIPFLD / 0.741485316D+00 /
      SAVE DIPFLD

      DATA B1234E_F / -1.646253575 /
      SAVE B1234E_F
      DATA B12B1_F / 0.005D+00 /
      SAVE B12B1_F 
      DATA B12B2_F / -1.642961068D+00/
      SAVE B12B2_F
      DATA B1G2_F / 1.509417504D+00 /
      SAVE B1G2_F
      DATA B2G2_F / -1.556509D+00 /
*      DATA B2G2_F / 0.0D+00 /
      SAVE B2G2_F

* Jacina kvadrupola (T/m)


      
* Magnetni region 

      SAVE NUM_135
      SAVE NUM_30
      SAVE NUM_45
*      SAVE NUMQUAD
      
* array for roatations
      DIMENSION LAT_ROT(100)
* array for magnetic field  values
      DOUBLE PRECISION FIELD_VA(100)
* array for informatiio is it qaud or dipole magnet
      DIMENSION LAT_VER(100)
      DATA LAT_ROT / 100*0 /
      DATA FIELD_VA / 100*0. /
      DATA LAT_VER /100*0 /
      SAVE LAT_ROT
      SAVE FIELD_VA
      SAVE LAT_VER

* Standardne vrednosti
      IDISC = 0
      BTX = ZERZER
      BTY = ONEONE
      BTZ = ZERZER
      B   = ZERZER
      
*----------------------------------------------------------------------*
*     first call: initialisation
*----------------------------------------------------------------------*

      IF ( BFIRST ) THEN
*           some useful printout
            WRITE (LUNOUT,*) ''
            WRITE (LUNOUT,*) ' dedicated magfld.f routine'
            WRITE (LUNOUT,*) ' initialisation started'
            WRITE (LUNOUT,*) ''
            WRITE (LUNOUT,*) ''
            WRITE (LUNOUT,*) ''
            CALL GEON2R ("PIPE_135", NUM_135, IERR )
            
            
            WRITE (LUNOUT,*) ''
            WRITE (LUNOUT,*) ' "PIPE_135" region je broj: ',NUM_135
            WRITE (LUNOUT,*) ' B2G2 pipe ',B2G2_F 
            WRITE (LUNOUT,*) ''
            
*            CALL GEON2R("Q___PIPE", NUMQUAD, IERR )
*            WRITE (LUNOUT,*) 'Quad cevka je broj: ', NUMQUAD

            CALL GEON2R("D30_PIPE", NUM_30, IERR )
            WRITE (LUNOUT,*) 'Dipol 30 cevka je broj: ', NUM_30

************************
************************
* 
* THIS IS THE PROBLEMATIC MAGNET
*
*************************
*************************

            CALL GEON2R("PIPE__45", NUM_45, IERR )
            WRITE (LUNOUT,*) 'Dipol pipe: ', NUM_45
            
            
            

********************
* DIPOL
********************
            CALL GEON2L("R_B1234A", NLATT, IRTLAT, IERR)
            IF ( IERR .GT. 0 ) THEN
                  CALL FLABRT('MAGFLD','no "R_B1234A" lattice')
            ENDIF
            LAT_ROT( NLATT ) = IRTLAT
            LAT_VER( NLATT ) = 2
            FIELD_VA( NLATT ) = B1234E_F
            WRITE (LUNOUT,*) 'DIPOLE B1234A je latica broj: ', NLATT

            CALL GEON2L("R_B1234B", NLATT, IRTLAT, IERR)
            IF ( IERR .GT. 0 ) THEN
                  CALL FLABRT('MAGFLD','no "R_B1234B" lattice')
            ENDIF
            LAT_ROT( NLATT ) = IRTLAT
            LAT_VER( NLATT ) = 2
            FIELD_VA( NLATT ) = B1234E_F
            WRITE (LUNOUT,*) 'DIPOLE B1234B je latica broj: ', NLATT

            CALL GEON2L("R_B1234C", NLATT, IRTLAT, IERR)
            IF ( IERR .GT. 0 ) THEN
                  CALL FLABRT('MAGFLD','no "R_B1234C" lattice')
            ENDIF
            LAT_ROT( NLATT ) = IRTLAT
            LAT_VER( NLATT ) = 2
            FIELD_VA( NLATT ) = B1234E_F
            WRITE (LUNOUT,*) 'DIPOLE B1234C je latica broj: ', NLATT
            
            CALL GEON2L("R_B1234D", NLATT, IRTLAT, IERR)
            IF ( IERR .GT. 0 ) THEN
                  CALL FLABRT('MAGFLD','no "R_B1234D" lattice')
            ENDIF
            LAT_ROT( NLATT ) = IRTLAT
            LAT_VER( NLATT ) = 2
            FIELD_VA( NLATT ) = B1234E_F
            WRITE (LUNOUT,*) 'DIPOLE B1234D je latica broj: ', NLATT




**********************
* DIPOL 30
**********************
            CALL GEON2L("R_B12B1A", NLATT, IRTLAT, IERR)
            IF ( IERR .GT. 0 ) THEN
                  CALL FLABRT('MAGFLD','no "R_B12B1A" lattice')
            ENDIF
            LAT_ROT( NLATT ) = IRTLAT
            LAT_VER( NLATT ) = 2
            FIELD_VA( NLATT ) = B12B1_F
            WRITE (LUNOUT,*) 'DIPOLE R_B12B1A je latica broj: ', NLATT



**********************
* DIPOL 30
**********************
            CALL GEON2L("R_B12B2A", NLATT, IRTLAT, IERR)
            IF ( IERR .GT. 0 ) THEN
                  CALL FLABRT('MAGFLD','no "R_B12B2A" lattice')
            ENDIF
            LAT_ROT( NLATT ) = IRTLAT
            LAT_VER( NLATT ) = 2
            FIELD_VA( NLATT ) = B12B2_F
            WRITE (LUNOUT,*) 'DIPOLE R_B12B2A je latica broj: ', NLATT
            
**********************
* DIPOL 30
**********************
            CALL GEON2L("R_B12B2B", NLATT, IRTLAT, IERR)
            IF ( IERR .GT. 0 ) THEN
                  CALL FLABRT('MAGFLD','no "R_B12B2B" lattice')
            ENDIF
            LAT_ROT( NLATT ) = IRTLAT
            LAT_VER( NLATT ) = 2
            FIELD_VA( NLATT ) = B12B2_F
            WRITE (LUNOUT,*) 'DIPOLE R_B12B2B je latica broj: ', NLATT 

**********************
* DIPOL 45
**********************

************************
************************
* 
* THIS IS THE PROBLEMATIC MAGNET
*
*************************
*************************

            CALL GEON2L("R_B1G2__", NLATT, IRTLAT, IERR)
            IF ( IERR .GT. 0 ) THEN
                  CALL FLABRT('MAGFLD','no "R_B1G2__" lattice')
            ENDIF
            LAT_ROT( NLATT ) = IRTLAT
* information that it is a dipole
            LAT_VER( NLATT ) = 2
* magnetic field 
            FIELD_VA( NLATT ) = B1G2_F
            WRITE (LUNOUT,*) 'DIPOLE R_B1G2__ je latica broj: ', NLATT




         CALL FFLUSH
         CALL FLUSH(LUNOUT)      

            
*           procisti  .out fajl
            BFIRST = .FALSE.
      END IF
      
*----------------------------------------------------------------------*
* core
*----------------------------------------------------------------------*

      XNEW = X
      YNEW = Y
      ZNEW = Z
      MLAT = MLATTC
      IF (MLAT .GT. 0 ) THEN
* Ovde radim transformaciju koordinata iz lattice u original
            CALL DOTRSF(1, XNEW, YNEW, ZNEW, LAT_ROT(MLAT))
      ENDIF
      
      IF (NREG .EQ. NUM_135 ) THEN
* Ovo je 135 degree dipol
            BTX = ZERZER
            BTY = B2G2_F
            BTZ = ZERZER
            WRITE (LUNOUT,*) ' 135 Degree Magnet '
      ELSEIF (LAT_VER(MLAT) .EQ. 2) THEN
            BTX = ZERZER
            BTY = FIELD_VA(MLAT)
            BTZ = ZERZER
      ELSE  
*	 ni jedan drugi region nije magnetni
            CALL FLABRT('MAGFLD','NO MAGNETIC FIELD IN THIS REGION !')
      ENDIF
      
      B = SQRT(BTX**2 + BTY**2 + BTZ**2)
      
*     check against numerical precision
      IF ( B .GT. 1.0D-12 ) THEN
*     normalise field components
            BTX = BTX / B
            BTY = BTY / B
            BTZ = BTZ / B
      ELSE
*     restore default values
            BTX = ZERZER
            BTY = ZERZER
            BTZ = ONEONE
            B   = ZERZER
      ENDIF
      
      IF (MLAT .GT. 0) THEN
            CALL UNDRTO (1, BTX, BTY, BTZ, LAT_ROT(MLAT))
      ENDIF


      RETURN
*=== End of subroutine Magfld =========================================*
      END

