program field_check ! reduced version of tracking_master for calculating fields of elements. use bmad use muon_mod use muon_interface use parameters_bmad use runge_kutta_mod use materials_mod use quad_scrape_parameters use random_mod implicit none type (lat_struct), target:: lat, low_e_lat,high_e_lat, matchlat type (lat_struct) latb type (branch_struct), pointer :: branch ! ring type (branch_struct), pointer:: from_branch ! injection line type (ele_struct), pointer :: fork_ele, ele type (ele_struct) ele_at_s, ele_save type (ele_struct), pointer :: lord_ele type (ele_struct), allocatable :: ele_time_bins(:) type (coord_struct), allocatable :: from_orbit(:) type (coord_struct), allocatable, save :: co(:), all_orb(:) type (coord_struct), allocatable :: coSteps(:),co_electron(:) type (coord_struct), allocatable :: low_e_orb(:), high_e_orb(:) type (coord_struct) to_orbit type (coord_struct) start_orb, orb_at_s, co_start, co_end type (coord_struct) opt_orb type (coord_struct) co_muon_end type (coord_struct) co_elec_end type (em_field_struct) field type (track_struct) track type (sumBetaCrossE_struct), allocatable :: sumBetaCrossE(:) procedure(track1_custom_def):: track1_custom procedure(em_field_custom_def) :: em_field_custom procedure(wall_hit_handler_custom_def) :: wall_hit_handler_custom procedure(track1_preprocess_def) :: track1_preprocess procedure(track1_postprocess_def) :: track1_postprocess ! procedure(make_mat6_custom_def) :: make_mat6_custom integer pgopen, istat1, istat2 integer i, j, k,l,steps/1000/ !steps: number of steps for tracking the muon decay integer esteps/100/,tot_track/7/,n_row !used in the track_electron_s subroutine integer ix, i_dim/4/ integer ie_from integer j1 integer lun,lun32 integer nargs, iargc, ios integer end, nturns, nmuons, nmuon_first integer nmu integer bin,bins/0/ integer n integer track_state integer unit integer turn/0/ integer ix_ele integer tot_inf_scatter, tot integer ix_ele_b, ix_ele_c, ix_ele_k integer status/0/ integer nbranch, inj_branch/0/ integer vparam_id/0/, vparam_cbo_turns/10/ integer ix_end integer seed/1/ integer i2,i3 integer lunclose integer datetime_values(8) integer system integer imin/0/ integer write_phase_space_every_n_turns/0/ integer ndum real(rp) xmean_m, xsigma_m, ymean_m, ysigma_m, tmean, tsigma, pxmean, pxsigma, pymean, pysigma, pzmean, pzsigma real(rp) epsx, epsy, tlength, pz real(rp), allocatable :: NN(:) real(rp) circumference, ele_len,circumference2 real(rp) distance real(rp) kickermagnitude_0, delta_kick real(rp) angle real(rp) l_angle_b, l_angle_c, delta_angle real(rp) ring_theta/0./ real(rp) delta_e/1.e-3/,chrom_x,chrom_y,pz_ref/0.0/ real(rp) vparam_max_tot/0./, cbo_min/100./,vparam_cbo_min/0./ real(rp) vparam, vparam_max, delta_vparam, vparam_min real(rp) mat_inv(6,6), mat(6,6) real(rp) s_end real(rp) s1,s2 real(rp) inflector_angle real(rp) s_rel,t_rel real(rp) Q(3) real(rp) dt, dt1, dt2 real(rp) ele_limit(100,4) real(rp) s_start, delta_s real(rp) time_bin_width/0./ real(rp) xmax/0./, ymax/0./ real(rp) kick(10)/10*0./ !upstream steering real(rp) muon_lifetime real(rp) DeltaB_onB, Delta_Energy, E_tot real(rp) pz_save real(rp) BetaCrossE_start_time/0/ integer tot_max/0/ real(rp) t_save real(rp) r_save real(rp) time_binning_start/10.e-9/ !time when start to be by time real(rp) peak_to_peak real(rp) coef(4),chisq real(rp) x,y,s,dx,t integer ndf character*120 line, lat_file, lat_file_name/''/,new_file/' '/, muon_file/' '/ character*120 string character*16 ele_name character*20 vparam_label(23)/'beta x [m]','beta y [m]','eta [m]','kick [G]','kick width [ns]', & 'initial offset [mm]','initial angle [mrad]','kick rise time [ns]','kick fall time [ns]','inflector [Bmagic]', & 'energy offset [%]', ' etap ',' alpha x','kick_tStart [ns]','kick(1) G','kick(2) G','kick(3) G','quad [n]', & 'inf angle','rf_quad%amp_h','rf_quad%amp_v','rf_quad%phi_h','rf_quad%phi_v' / character*16 date, time, zone character *120 azimuthal_exp_z, azimuthal_exp_r character *120 azimuthal_exp_z_nom/'BzFourier20170628_LogID983.dat'/, azimuthal_exp_r_nom/'BrFourier2016.dat'/ character*16 steering(10)/10*' '/ !name of upstream steering character*3 char_turn integer vparam_column_index(23)/17,20,18,19,21,24,25,22,23,26,27,28,29,30,1,2,3,31,32,1,5,9,13/ logical err_flag, inf_end_us/.false./, inf_end_ds/.false./ logical first/.true./, twiss_file/.false./ logical create_new_distribution/.false./ logical loop/.false./,gnu_input_only/.false./ logical inTracker1/.false./ logical first_loss/.true./ logical enerloss/.false./, quad_p/.false./ logical spin_tracking_on/.true./,muon_decay_on/.true./ logical local_ref/.false./,calcd/.false./ !why was this removed logical opt_incident/.false./,inj_matrix_tracking/.false./ logical start_tracking_at_inflector_exit/.false./ logical make_movie/.false./ logical use_lattice_twiss/.false./ logical write_phase_space_file/.false./ logical save_element_by_element_info/.true./ logical FiberScatter0/.false./, FiberEnergyLoss0/.false./ logical error_fields/.false./ logical err logical int_betaCrossE/.false./, betaCrossE_verbose logical initial_offset_from_co/.false./ logical itexists logical save_scraping_status logical opt_inf_exit_angle/.true./ logical em_field_calc_verbose0 type (muon_struct), allocatable :: muons(:), muons_ele(:,:), muons_temp(:), muons_ele_inj(:,:) type (electron_struct),allocatable::electrons_s(:,:) type (g2twiss_struct) twiss, twiss_opt, twiss_match, twiss_start_inj_line type (g2moment_struct) moment, moment_ele type (initial_offsets_struct) initial_offsets,initial_offsets_ref, inflector_end_target type (kicker_params_struct) kicker_params, kicker_params_0 type (field_file_struct) fringe_file, inflector_file type (rf_quad_struct) rf_quad_0(4) real(rp), dimension(100) :: fnal_w_integral real(rp) start, finish character*16 epsdistr, tdistr, pzdistr, inf_aperture, twiss_ref, ring_twiss, input_call/''/ ! Input namelist structure ("/g-2/test/input.dat") namelist /input/ lat_file_name, nmuons, nturns, & ! GENERAL nmuon_first, & ! when reading from a list start here create_new_distribution, new_file, muon_file, & ! BEAM (write and/or read from file) tdistr, tlength, tsigma, & ! BEAM (longitudinal) pzdistr, pz, pzsigma, & ! BEAM (energy) epsdistr, epsx, epsy, twiss, twiss_ref, & ! BEAM (transverse) inf_aperture, inf_end_us, inf_end_ds, initial_offsets, initial_offsets_ref, enerloss, inflector_angle, inflector_field, & ! INFLECTOR inflector_end_target, & !opt_incident finds initial offsets so that on momentum muon exits inflector with these target values ring_theta, & ! ELEMENT POSITIONING kickerPlates, kickerCircuit, kickerPulseFile, kickerFieldType, kicker_sextupole_scale_factor, & ! KICKER kicker_params, & ! KICKER kicker_tStart, & ! KICKER quad_plate, & ! QUAD SCATTERING quadPlates, quadCircuit, quadFieldType, & ! QUAD quad_params, & ! QUAD twiss_file, & loop, vparam_max, vparam_id, delta_vparam, vparam_min, vparam_cbo_turns, & gnu_input_only,spin_tracking_on, muon_decay_on, inflector_width, opt_incident, & inj_matrix_tracking, & fringe_file, inflector_file, & ! names of field maps for fringe and inflector and a few other details about the maps start_tracking_at_inflector_exit, & !track distribution from inflector exit rather than start of injection line ring_twiss, & ! ring_twiss = 'open' use input twiss as starting point for propagation around ring. ring_twiss='closed' compute closed ring seed, & make_movie, & use_lattice_twiss, & ! to compute twiss parameters through injection line, if true, start with values defined in lattice file write_phase_space_file, & ! to write a file with muon phase space at each element on first turn save_element_by_element_info, & !if true then allocate array muons_ele(:,:) which saves distribution at each element rf_quad, & !parameters of rf quads em_field_calc_verbose, & !locigal. if true write quad dipole field and phase space coordinate for each step along trajectory, starting after initialization time_bin_width, & !bin for tbt moments vs time scraping_on, init_quad_focus, init_quad_steer, quad_ramp_start_time, quad_ramp_end_time, & B_radial, & !turn on quad scraping, scraping focus fraction, scraping steering fraction, start time and end time FiberScatter0, FiberEnergyLoss0, & !if true turn on scattering in harp fibers, default = false error_fields, & ! if true call field_errors to introduce errors in g-2 dipole field, default = false azimuthal_exp_z, azimuthal_exp_r, & ! fourier expansion of z and r fields steering, kick, & !name and strength of upstream steering int_betaCrossE, betaCrossE_verbose, & !if int_ true integrate BetaCrossEfield along trajectory of each muon, if verbose save trajectory DeltaB_onB, no_energy_change, & ! fractional field offset and fractional momentum offset. Always both at the same time no_energy_change, & ! if true then set vec(6) on quad exit ot vec(6) at entrance initial_offset_from_co, & ! if true, offset from closed orbit (rather than zero) quad_fringe_energy_change, & ! if true, offset vec(6) at entrance and exit of quad proportional to Efield input_call,& ! more info BetaCrossE_start_time, & ! time to start integrating beta X E. 35 microseconds write_phase_space_every_n_turns, & ! write phase space every n turns. If n=0, do not write bad_resistors, extra_inf_scatter_factor, & !bad resistors T or F. Bad resistors requires scraping. extra inflector scatter factor normally 0. bigger number more scatter horizontal_scrape, quad_scrape_rf_verbose, & !if true standard scraping. If false no horizontal scraping, only vertical. Default is horizontal_scrape = scraping_on, if verbose write field vs time from track1_preprocess multipole_field_file, & ! file with 2d multipoles for magnetic field at 8600 azimuthal angles pbin_width, & ! if >0 bin spin according to time and momentum. this is the width of the momentum bin time_binning, & !'position_slice': binned by time at the end of the lattice,'momentum_slice': binned by momentum set_polarization, & ! if true then when reading phase_space file to track, all polarizations are set to s/|p| = -p/|p| polarization !if set_polarization is T, then 'antiparallel' or 'spintune' ! Read the lattice file. The default name is "bmad." The lattice file contains the layout of the ! ring elements, including length, locations and apertures of kicker and quads. ! The lattice file also defines beam energy, field of the main dipole (in terms of bending radius), ! and kicker amplitude and width call cpu_time(start) multipole_field_file = '' extra_inf_scatter_factor =0. quad_plate = .false. scraping_on=.false. horizontal_scrape= .true. quad_scrape_rf_verbose=.false. em_field_calc_verbose=.false. bad_resistors = .false. fixed_quad_time = .false. kicker_params%dtRise = 0. kicker_params%dtFall = 0. quad_params%short_quad_field_index = -99. quad_params%long_quad_field_index = -99. quad_params%short_quad_plate_index%inner = 0. quad_params%short_quad_plate_index%bottom = 0. quad_params%short_quad_plate_index%outer = 0. quad_params%short_quad_plate_index%top = 0. quad_params%long_quad_plate_index = quad_params%short_quad_plate_index rfquad_params_reset=.true. rf_quad(1:4)%amp_h=0 rf_quad(1:4)%amp_v=0 betaCrossE_verbose = .false. DeltaB_onB=0 Delta_Energy=0 quad_fringe_energy_change = .false. set_polarization = .false. ix_end_flag = 0 !do not start to accumulate in compute_time_dep_moments_calo until ix_end_flag = 1 call date_and_time(date,time,zone,datetime_values) print '(a10,a10,a1,a10)',' date = ',date,' time = ', time directory = trim(date)//'_'//trim(time(1:6)) status= system('mkdir ' // directory) if(status == 0)print '(a20,a)',' create directory = ', trim(directory) OPEN (UNIT=5, FILE='input.dat', STATUS='old', IOSTAT=ios) READ (5, NML=input, IOSTAT=ios) print *, 'ios=', ios rewind(unit=5) ! READ (5, NML=input) CLOSE(5) if(input_call /= '')then open(unit=5,file=input_call, status='old', iostat=ios) read(5,nml=input, iostat=ios) string = 'cp '//input_call//' '//trim(directory)//'/.' call execute_command_line(string) endif OPEN (UNIT=5, FILE='input.dat', STATUS='old', IOSTAT=ios) lun=lunget() open(unit=lun, file=trim(directory)//'/input.dat') do while(ios == 0) read(5,'(a)', IOSTAT=ios)string if(ios == 0)write(lun,'(a)')string end do close(unit=lun) WRITE(6,NML=input) print '(a,es16.8)',' muon anomalous magnetic moment = ', anomalous_moment_of(antimuon$) print '(a,es15.8)',' anomalous moment of muon ',anomalous_moment_of(antimuon$) ! Check consistency of input parameters number_turns = nturns bin_width = time_bin_width tilt = inflector_angle ! save initial kicks kicker_params_0%kicker_field(1:4) = kicker_params%kicker_field(1:4) ! save initial quads quad_params_0%short_quad_field_index(1:4) = quad_params%short_quad_field_index(1:4) quad_params_0%long_quad_field_index(1:4) = quad_params%long_quad_field_index(1:4) quad_params_0%short_quad_plate_index(1:4) = quad_params%short_quad_plate_index(1:4) quad_params_0%long_quad_plate_index(1:4) = quad_params%long_quad_plate_index(1:4) rf_quad_0(1:4)%amp_h = rf_quad(1:4)%amp_h rf_quad_0(1:4)%amp_v = rf_quad(1:4)%amp_v rf_quad_0(1:4)%phi_h = rf_quad(1:4)%phi_h rf_quad_0(1:4)%phi_v = rf_quad(1:4)%phi_v em_field_calc_verbose0 = em_field_calc_verbose em_field_calc_verbose=.false. !!! rf_quad(1:4)%amp_h=0 !!! rf_quad(1:4)%amp_v=0 !turn on spin tracking if (spin_tracking_on) bmad_com%spin_tracking_on = .true. ! set names for field_files field_file(1) = fringe_file field_file(2) = inflector_file lat_file = lat_file_name print *, ' lat_file = ', lat_file !Copy lattice file to working directory string = 'cp'//' '//trim(lat_file)//' '//trim(directory)//'/.' print *,string call execute_command_line(string) inquire(file='quad_plate_misalign.bmad' , exist=itexists) if(itexists)then string = 'cp'//' '//'quad_plate_misalign.bmad'//' '//trim(directory)//'/.' print *,string call execute_command_line(string) endif ! Construct the lattice from the file ! There are two branches, the injection line (branch=0) from somewhere upstream to the inflector exit ! and the ring (branch=1) track1_custom_ptr => track1_custom em_field_custom_ptr => em_field_custom wall_hit_handler_custom_ptr => wall_hit_handler_custom track1_preprocess_ptr => track1_preprocess track1_postprocess_ptr => track1_postprocess !make_mat6_custom_ptr => make_mat6_custom bmad_com%auto_bookkeeper=.false. bmad_com%min_ds_adaptive_tracking = 0.00001 ! bmad_com%rel_tol_adaptive_tracking = 1d-10 ! Runge-Kutta tracking relative tolerance. ! bmad_com%abs_tol_adaptive_tracking = 1d-12 ! Runge-Kutta tracking absolute tolerance. bmad_com%rel_tol_tracking= 1d-8 bmad_com%abs_tol_tracking= 1d-10 ! Scattering in inflector etc. is off while parser computes transfer matrices call initializeMaterials() save_scraping_status = scraping_on scraping_on = .false. call bmad_parser (lat_file, lat) !initialize error fields lat%param%ixx = 0 if(error_fields)lat%param%ixx = 1 scraping_on = save_scraping_status nbranch = size(lat%branch) -1 ! ring branch print *,' nbranch = ', nbranch ! Read the input namelist ("/g-2/files/input.dat") tilt = inflector_angle ! Turn off scattering until reference orbit and transfer matrices are established if(no_energy_change .and. (enerloss .or. FiberEnergyLoss0))then print *, 'NO_ENERGY_CHANGE and ELOSS or FIBERENERGYLOSS are true. This cannot be' stop endif lun=lunget() open(unit=lun,file=trim(directory)//'/'//'input_data.dat') write(lun,nml=input) close(lun) branch => lat%branch(nbranch) ! this is the ring branch from_branch => lat%branch(0) ! injection line branch ie_from = branch%ix_from_ele fork_ele => from_branch%ele(ie_from) circumference = branch%ele(branch%n_ele_track)%s !initialize error fields branch%param%ixx = 0 if(error_fields)then branch%param%ixx = 1 if(azimuthal_exp_z == ' ') azimuthal_exp_z=azimuthal_exp_z_nom if(azimuthal_exp_r == ' ') azimuthal_exp_r=azimuthal_exp_r_nom call init_and_write_generalized_cyl_exp(azimuthal_exp_z, azimuthal_exp_r) endif print *,' Fork_ele assigned: fork_ele%name = ', fork_ele%name call set_fringe_params(lat, 0, DeltaB_onB) call set_inflector_params(lat, 0,inflector_angle) call set_quad_params(lat, nbranch) call set_dipole_params(lat, nbranch, DeltaB_onB) lat%branch(0)%ele(0)%value(p0c$)=lat%branch(0)%ele(0)%value(p0c$)*(1.+Delta_Energy) lat%branch(1)%ele(0)%value(p0c$)=lat%branch(0)%ele(0)%value(p0c$) call set_flags_for_changed_attribute(lat%branch(0)%ele(0),lat%branch(0)%ele(0)%value(p0c$)) call set_flags_for_changed_attribute(lat%branch(1)%ele(0),lat%branch(1)%ele(0)%value(p0c$)) call lattice_bookkeeper(lat) rf_quad(:)%amp_h=0 rf_quad(:)%amp_v=0 rfquad_params_reset=.true. ! Calculate the ring's Twiss parameters, coupling, and dispersion parameters scraping_on = .false. call closed_orbit_calc(lat, co, i_dim=4, ix_branch=nbranch) do i=1,branch%n_ele_track call lat_make_mat6(lat,ix_ele = i,ref_orb = co, ix_branch=nbranch) end do call twiss_at_start(lat, status=status, ix_branch=nbranch) open(lun, file=trim(directory)//'/'//'co_tunes_scraping_off.dat') write(lun,'(a)')' scraping off' write(lun, '(2(a5,es12.4,1x),3(a8,es12.4,1x))')'Qx = ',lat%a%tune/twopi, 'Qy = ',lat%b%tune/twopi,' Betax = ', lat%branch(1)%ele(0)%a%beta,& ' Betay = ', lat%branch(1)%ele(0)%b%beta,' eta = ', lat%branch(1)%ele(0)%x%eta call update_quad_params(vparam,vparam_min) call set_quad_params(lat, nbranch) print '(a,4es12.4)',' short_quad_index ',quad_params%short_quad_field_index(1:4) print '(a,4es12.4)',' long_quad_index ',quad_params%long_quad_field_index(1:4) ! Set the kicker fields to zero to compute ring twiss parameters call set_kicker_params (lat,kicker_params) call set_steering(lat,steering,kick) ! call write_header(lat,nbranch,twiss, kicker_params,tlength, pz, initial_offsets) bmad_com%aperture_limit_on = .true. ! bmad_com%aperture_limit_on = .false. bmad_com%max_aperture_limit = 1. open(unit=11,file = trim(directory)//'/'//'bfield_vs_xyst.dat') do i=1,branch%n_ele_max ele=>branch%ele(i) if(index(ele%name,'KICKER1')==0)cycle x=0 s=ele%s - 0.1 dx=0.001 dt=1.e-9 do j=-100,100 t=j*dt write(11,'(/,a,es12.4,/)')'# t=',t do n = -44,44 y = n*dx call fields(ele,x,y,s,t,field) write(11, '(4es12.4, 1x,3es12.4)')x,y,s,y,field%B end do end do end do close(unit=11) end