MODULE parameters_bmad

use bmad_struct
use precision_def
use magfield
IMPLICIT NONE

! FUNDAMENTAL CONSTANTS (PDG 2013) AND EXPERIMENTAL PARAMETERS (E821 LITERATURE)
REAL, PARAMETER :: me     = 0.510998918e6  ! mass of electron (eV)
REAL, PARAMETER :: mmu    = 105.6583715e6  ! mass of muon (eV)
REAL, PARAMETER :: tauMu  = 2.1969811e-6   ! lifetime of muon (sec) in rest frame
REAL, PARAMETER :: amu    = 11659208.9e-10 ! anomalous magnetic moment of muon
REAL, PARAMETER :: rMagic = 7.112          ! design radius (meters)
! DERIVED QUANTITIES
REAL, PARAMETER :: gammaMagic = SQRT(1.+1./amu)           ! dimensionless
REAL, PARAMETER :: betaMagic  = SQRT(1.-1./gammaMagic**2) ! dimensionaless
REAL, PARAMETER :: EMagic     = gammaMagic*mmu            ! eV
REAL, PARAMETER :: pMagic     = gammaMagic*mmu*betaMagic  ! eV
REAL, PARAMETER :: BMagic     = pMagic/(c_light*rMagic)   ! Tesla

! CUSTOM ATTRIBUTES
INTEGER, PARAMETER :: kick_width$  = custom_attribute1$
INTEGER, PARAMETER :: field_index$ = custom_attribute2$
INTEGER, PARAMETER :: x_angle$ = custom_attribute3$
INTEGER,PARAMETER::kicker_field$ = custom_attribute4$
INTEGER,PARAMETER::dtRise_and_fall$ = custom_attribute5$

! INJECTION LINE
type(field_file_struct) field_file(2)

! INFLECTOR PARAMETERS
REAL(RP) inflector_width/ 0.009/    ! horizontal half-width
REAL(RP), PARAMETER :: inflector_height = 0.028    !  vertical  half-width
REAL(rp), PARAMETER :: inf_length       = 1696.e-3 ! inflector coil-to-coil length
REAL(rp), PARAMETER :: map_tilt = 0.00235    ! angle of inflector field map with respect to tangent line. Rotate about downstream end.
REAL(rp) tilt     ! angle of inflector with respect to tangent line. Rotate about downstream end. If tilt = map_tilt the inflector map is tilted 2.35mrad. 
real(rp) inflector_field  ! units of Bmagic
LOGICAL eloss, us_scatter, ds_scatter

! KICKER PARAMETERS
INTEGER  kickerPlates/1/     ! Controls shape of B-field: 1=uniform, 821, 989
INTEGER  kickerCircuit/1/    ! Controls pulse shape: 1=perfect(square kick), 821=RLC, 989=Blumlein
INTEGER  kickerFieldType/1/  ! 1=analytic, 2=mapped
REAL(RP) kicker_tStart(3)    ! kicker start times

! E821 KICKER RLC CIRCUIT PARAMETERS (E821)
REAL(RP), PARAMETER :: kickerR = 11.9    ! Ohm
REAL(RP), PARAMETER :: kickerL =  1.6e-6 ! Henry
REAL(RP), PARAMETER :: kickerC = 10.1e-9 ! Farad
REAL(RP)  kickerE821_omega, kickerE821_tau, kickerE821_tStartToPeak, kickerE821_normalization

! QUAD PARAMETERS
INTEGER quadPlates/821/   ! Controls shape of E-field (Q1): 821, 989
INTEGER quadCircuit/1/    ! 0=off at injection => ramped up (tauRC=5us); 1=storage at injection, 2=scraping at injection => storage (tauRC=5us)
INTEGER quadFieldType/1/  ! 1=analytic, 2=mapped

! MULTIPOLE COEFFICIENTS
REAL(RP), DIMENSION(:,:),   ALLOCATABLE :: KMC ! kicker
REAL(RP), DIMENSION(:,:,:), ALLOCATABLE :: QMC ! quad

! These data are the normalized integral of the "Fermilab W" shown in
! CDR Figure 7.5, lower left pane.  These data are used to generate
! the E989 longitudinal beam profile if the user selects tdistr="e989".
! The assumption here is that protons =~ pions =~ muons for the distr's.
REAL(RP), PARAMETER, DIMENSION(100) :: fnalw_integral = [          &
  0.000000000, 0.000825196, 0.004630266, 0.009031312, 0.018887819, &
  0.030623940, 0.043643699, 0.056617613, 0.068124513, 0.078118553, &
  0.085866227, 0.092697016, 0.099023518, 0.105029111, 0.110897171, &
  0.116627699, 0.122266538, 0.127767845, 0.133177463, 0.138403704, &
  0.143446569, 0.148260212, 0.152936322, 0.157474900, 0.161875946, &
  0.166231147, 0.170632192, 0.175124926, 0.179709348, 0.184431302, &
  0.189244946, 0.194150277, 0.199376519, 0.205382112, 0.212212901, &
  0.220877458, 0.231375785, 0.243295283, 0.256315041, 0.271260258, &
  0.287580800, 0.305551735, 0.324852152, 0.346032183, 0.368358318, &
  0.391463806, 0.415256957, 0.439325173, 0.463530922, 0.487828359, &
  0.512171641, 0.536469078, 0.560674827, 0.584743043, 0.608536194, &
  0.631641682, 0.653967817, 0.675147848, 0.694448265, 0.712419200, &
  0.728739742, 0.743684959, 0.756704717, 0.768624215, 0.779122542, &
  0.787787099, 0.794617888, 0.800623481, 0.805849723, 0.810755054, &
  0.815568698, 0.820290652, 0.824875074, 0.829367808, 0.833768853, &
  0.838124054, 0.842525100, 0.847063678, 0.851739788, 0.856553431, &
  0.861596296, 0.866822537, 0.872232155, 0.877733462, 0.883372301, &
  0.889102829, 0.894970889, 0.900976482, 0.907302984, 0.914133773, &
  0.921881447, 0.931875487, 0.943382387, 0.956356301, 0.969376060, &
  0.981112181, 0.990968688, 0.995369734, 0.999174804, 1.000000000  &
]

! Integral of the Landau-Vavilov distribution (centered at x=0 and with a FWHM of deltax=4.01865).
! This integral is used for randomly generating energy fluctuations during scattering.  See PDG(2013)
! Chapter 27, "Passage of Particles Through Matter" -- in particular the text surrounding Fig. 27.7.
REAL(RP), PARAMETER :: LandauVavilov_integral(-5:100) = [ &
  0.000000, 0.000000, 0.000012, 0.006529, 0.082049, &
  0.246723, 0.417628, 0.551239, 0.647017, 0.714968, &
  0.764007, 0.800302, 0.827872, 0.849327, 0.866387, &
  0.880213, 0.891605, 0.901130, 0.909194, 0.916100, &
  0.922072, 0.927283, 0.931866, 0.935925, 0.939542, &
  0.942785, 0.945708, 0.948355, 0.950763, 0.952961, &
  0.954977, 0.956830, 0.958540, 0.960123, 0.961592, &
  0.962959, 0.964233, 0.965424, 0.966540, 0.967591, &
  0.968572, 0.969500, 0.970375, 0.971203, 0.971986, &
  0.972728, 0.973432, 0.974101, 0.974738, 0.975828, &
  0.976385, 0.976916, 0.977425, 0.977912, 0.978378, &
  0.978826, 0.979255, 0.979668, 0.980065, 0.980446, &
  0.980814, 0.981168, 0.981509, 0.981839, 0.982157, &
  0.982464, 0.982760, 0.983048, 0.983325, 0.983594, &
  0.983854, 0.984107, 0.984352, 0.984589, 0.984819, &
  0.985042, 0.985259, 0.985470, 0.985675, 0.985875, &
  0.986069, 0.986257, 0.986441, 0.986620, 0.986794, &
  0.986964, 0.987129, 0.987291, 0.987448, 0.987602, &
  0.987752, 0.987898, 0.988041, 0.988181, 0.988317, &
  0.988450, 0.988581, 0.988708, 0.988833, 0.988955, &
  0.989074, 0.989191, 0.989305, 0.989417, 0.989527, &
  1.000000  &
]


! TODO: Get rid of all of the STUFF we don't need & aren't using!


! QUADRUPOLE MULTIPOLE EXPANSION COEFFICIENTS
! See Table 5 of the E821 quad NIM paper
! REAL(RP), PARAMETER :: neff = 0.141 ! Effective field index for E989
 REAL(RP), PARAMETER :: Qb2  = -(1./0.141)*20177.8/0.045**2 ! Reasonable approximation
 REAL(RP), PARAMETER :: Qb4  = -(1./0.141)*   33.0/0.045**4 ! Reasonable approximation
 REAL(RP), PARAMETER :: Qb6  =  (1./0.141)*   45.9/0.045**6 ! Reasonable approximation

! Physical constants:
REAL(RP), PARAMETER :: muonmass_MeV = 105.658367
REAL(RP), PARAMETER :: muonmass_kg = 1.883531475e-28
REAL(RP), PARAMETER :: muoncharge_C = 1.602176487e-19
REAL(RP), PARAMETER :: speedoflight = 299792458. ! in meters/second
REAL(RP), PARAMETER :: omega_a = 1439311.81495     !radians/sec
REAL(RP), PARAMETER :: lifetime = 64.4e-6       !dilated lifetime at magic momentum, seconds
! Assuming times measured in nanoseconds, here are conversions for the units used in the program
!     All inputs must be converted to these units before being used by derivs.f90, because the
!     differential equations in that subroutine are written in dimensionless form
REAL(RP), PARAMETER :: timeunit_s = 1e-9
REAL(RP), PARAMETER :: lengthunit_m = speedoflight*timeunit_s
REAL(RP), PARAMETER :: momentumunit_MeVperc = muonmass_MeV
REAL(RP), PARAMETER :: Bfieldunit_T = muonmass_kg/(muoncharge_C*timeunit_s)
REAL(RP), PARAMETER :: voltageunit_V = muonmass_kg*speedoflight**2/muoncharge_C


! Coordinates, dimensions and constants, converted to my units:
REAL(RP), PARAMETER :: storageradius  = 7.112
REAL(RP), PARAMETER :: magicmomentum  = 3094.35/momentumunit_MeVperc


REAL(RP), PARAMETER :: Q1platethickness = 0.0001/lengthunit_m ! Q1 plate is 0.1 mm thick


! MATERIALS
REAL(RP), PARAMETER :: radlength_Al   = 0.08897 ! radiation length for aluminum, 8.897 cm
REAL(RP), PARAMETER :: radlength_Cu   = 0.014   ! radiation length for copper
REAL(RP), PARAMETER :: radlength_Nb   = 0.012   ! radiation length for Nb
REAL(RP), PARAMETER :: radlength_Ti   = 0.036   ! radiation length for Ti
REAL(RP), PARAMETER :: thickness_Al   = 0.00458 ! window + flange + superconductor
REAL(RP), PARAMETER :: thickness_Cu   = 0.00039 ! thickness of copper
REAL(RP), PARAMETER :: thickness_NbTi = 0.00043 ! thickness of NbTi


CONTAINS

SUBROUTINE initializeKickerMultipoleCoeffs(KMC)
  IMPLICIT NONE
  REAL(RP), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: KMC ! kicker multipole coefficients
  ALLOCATE( KMC(2,0:24) )

  ! E821 (r0=4.5cm)
  KMC(1,0:24) = [ &
    1.00011,      -0.00013228,    0.0944713,    0.0000887371, -0.135094,     &
   -0.000121285,   0.0660964,     0.0000865873, 0.0146219,     0.0000463802, &
   -0.0212747,    -0.0000811448, -0.000991705,  0.0000202727,  0.00593946,   &
   -0.0000395268, -0.000607024,   0.000150024, -0.0013922,    -0.00011031,   &
   0.000265688,   -0.0000477355,  0.000519422, -0.000025101,  -0.00027104    &
  ]

  ! E989 (r0=4.0cm)
  KMC(2,0:24) = [ &
    0.999646,      0.000807006,  0.382047,     0.00115502,  -0.184812,    &
    0.000246815,  -0.205932,    -0.000904902, -0.0423524,   -0.000429464, &
    0.0620312,     0.0000345979, 0.0572378,    0.00057415,   0.00785289,  &
    0.0000509476, -0.0222868,   -0.000321909, -0.0180536,   -0.000164196, &
   -0.00178202,   -0.000122991,  0.00753034,  -0.000265512,  0.00630098   &
  ]
END SUBROUTINE initializeKickerMultipoleCoeffs


SUBROUTINE initializeQuadMultipoleCoeffs(QMC)
  IMPLICIT NONE
  REAL(RP), DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: QMC ! quad multipole coefficients
  integer, dimension(:), allocatable :: quad_multipoles
  integer lun, open_status
  namelist/quad_multipole_use_list/quad_multipoles
  ALLOCATE( QMC(0:4,2,14) )
  allocate(quad_multipoles(14))
  quad_multipoles(:)=1
  
lun=lunget()
open(lun, file='quad_multipoles.dat', iostat=open_status, status='old', action='read')
if(open_status == 0)read(lun,nml=quad_multipole_use_list)
if(open_status /= 0)print *,' No quad multipole file. Use default'
close(unit=lun)

  ! QUAD storage (all quads) -- Table 5 of the E821 quad NIM paper
  QMC(0,1,:) = [ 0.0, 20177.8, 0.0, 33.0, 0.0, -45.9, 0.0, -5.5, 0.0, -391.3, 0.0, -6.5, 0.0, 52.3 ] ! cos
  QMC(0,2,:) = [ 0.0,     0.1, 0.0,  0.1, 0.0,   0.1, 0.0, -0.2, 0.0,    0.1, 0.0,  0.0, 0.0, -0.1 ] ! sin

  QMC(0,1,:) = QMC(0,1,:)*quad_multipoles(:)
  
  ! QUAD1 scraping -- Table 7 of the E821 quad NIM paper
  QMC(1,1,:) = [     0.0, 18706.0,   0.0, 375.0,    0.0, -41.0,   0.0, -76.0,  0.0, -364.0,  0.0, -3.0,  0.0, 49.0 ] ! cos
  QMC(1,2,:) = [ -2155.0,     0.0, 837.0,   0.0, -117.0,   0.0, -57.5,   0.0, 52.0,    0.0, -8.0,  0.0, -5.0,  0.0 ] ! sin

  ! QUAD2 scraping -- Table 8 of the E821 quad NIM paper (with theta => theta+pi for odd n-values in cosine terms)
  QMC(2,1,:) = [  2150.0, 17235.0, 840.0, 28.0,  121.0, -39.0,  57.0, -5.0, -53.0, -334.0, -9.0, -6.0,  5.0, 45.0 ] ! cos
  QMC(2,2,:) = [ -2156.0,     0.0, 837.0,  0.0, -117.0,   0.0, -57.0,  0.0,  52.0,    0.0, -8.0,  0.0, -5.0,  0.0 ] ! sin

  ! QUAD3 scraping -- Table 7 of the E821 quad NIM paper
  QMC(3,1,:) = QMC(0,1,:) ! cos
  QMC(3,2,:) = QMC(0,2,:) ! sin

  ! QUAD4 scraping -- Table 8 of the E821 quad NIM paper
  QMC(4,1,:) = [ -2150.0, 17235.0, -840.0, 28.0, -121.0, -39.0, -57.0, -5.0, 53.0, -334.0, 9.0, -6.0, -5.0, 45.0 ] ! cos
  QMC(4,2,:) = QMC(2,2,:) ! sin


print '(a)',' Quad multipole coefficients initialized'
print '(a)',' quad_multipoles = '
print '(14es12.4)',QMC(0,1,:)
END SUBROUTINE initializeQuadMultipoleCoeffs


END MODULE parameters_bmad
