MODULE parameters_bmad use bmad_struct use precision_def use magfield use muon_mod !, dummy => track1_custom, dummy=>track1_preprocess IMPLICIT NONE type (ele_struct), allocatable, target :: ele_ref(:) CHARACTER*16 directory !where all output is written integer lun_quad_params logical no_energy_change, quad_fringe_energy_change logical set_polarization character*16 polarization !'antiparallel' or 'spintune' integer muon_number ! 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 ! 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 real(rp) extra_inf_scatter_factor !multiplies multiple scattering in end of inflector ! 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 character*60 kickerPulseFile/''/ INTEGER kickerFieldType/1/ ! 1=analytic, 2=mapped REAL(RP) kicker_tStart(4) ! 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 !FIBER PARAMETERS logical FiberScatter, FiberEnergyLoss ! 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 type (quad_params_struct) quad_params, quad_params_0 type (rf_quad_struct) rf_quad(4) logical horizontal_scrape, quad_scrape_rf_verbose logical rfquad_params_reset ! MULTIPOLE COEFFICIENTS REAL(RP), DIMENSION(:,:), ALLOCATABLE :: KMC ! kicker REAL(RP), DIMENSION(:,:,:), ALLOCATABLE :: QMC ! quad ! DIPOLE RADIAL FIELD REAL(RP) B_RADIAL character*120 multipole_field_file !IBMS logical ibms/.false./ real(rp) bin_width, pbin_width character*100 time_binning !'position_slice','momentum_slice' integer number_turns, ix_end_flag ! 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