program tracking_master 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 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 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 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 !if (muon_decay_on .and. .not. spin_tracking_on) then ! print *, "ERROR: With muon_decay_on = T, you must also set spin_tracking_on = T" ! call out_io (s_fatal$, 'tracking_master', 'ERROR IN INPUT FILE!') ! if (global_com%exit_on_error) call err_exit !end if !global_com%exit_on_error=.false. 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 nargs = command_argument_count() if (nargs == 1) then call get_command_argument(1, lat_file) !print *, 'Using ', trim(lat_file) else if(lat_file_name /= '')then lat_file = lat_file_name else lat_file = 'bmad.' print '(a,$)',' Lattice file name? (Default = bmad.) ' read(5,'(a)') line call string_trim(line, line, ix) lat_file = line if (ix==0) lat_file = 'bmad.' endif 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 ! call check_for_m5_to_inj_patch(lat, initial_offsets_ref) ! call check_for_m5_to_inj_patch(lat, initial_offsets) ! lat%ele(1:lat%n_ele_track)%alias= 'No' ! nbranch = ubound(lat%branch,1) !refers to the ring branch nbranch = size(lat%branch) -1 ! ring branch print *,' nbranch = ', nbranch ! Read the input namelist ("/g-2/files/input.dat") call create_gnu_input(tot_max,vparam_max_tot, cbo_min, vparam_cbo_min) !Initialize random number generator that is used to create phase space distribution and multiple scattering call ran_seed_put(seed) call ran_seed_get(seed) print *,' seed from ran_seed_get =',seed 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 eloss = .false. us_scatter = .false. ds_scatter = .false. quad_p = quad_plate quad_plate = .false. FiberScatter = .false. FiberEnergyLoss = .false. ! if(index(inflector_ape,'989')/=0)inflector_width=0.018 lun=lunget() open(unit=lun,file=trim(directory)//'/'//'input_data.dat') write(lun,nml=input) close(lun) ! Rotate ring (including all elements) with respect to inflector exit. ! A rotation of zero places the center of the central kicker 90 degrees from the inflector exit if(ring_theta /= 0.)then bmad_com%auto_bookkeeper=.false. ! call rotate_ring(lat%branch(nbranch),ring_theta) !this needs to be checked as to whether pass lat% branch will work endif 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 ! Set the Twiss parameters of lat%ele(0)=BEGINNING ! to the values in "input.dat" !Note we are asssigning these twiss parameters to first element in branch nbranch, namely the ringi ! These values are set in the input file branch%ele(0)%a%beta = abs(twiss%betax ) branch%ele(0)%b%beta = abs(twiss%betay) branch%ele(0)%a%alpha = twiss%alphax branch%ele(0)%b%alpha = twiss%alphay branch%ele(0)%x%eta = twiss%etax branch%ele(0)%x%etap = twiss%etapx branch%ele(0)%y%eta = twiss%etay branch%ele(0)%y%etap = twiss%etapy ! Set variables for determining Twiss parameters, coupling, dispersion, etc. of the ring branch 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) !Initialize ele_struct for subroutine time_bins if(trim(time_binning) == 'momentum_slice')then allocate(ele_time_bins(1:branch%n_ele_track)) ele_time_bins(1:branch%n_ele_track) = branch%ele(1:branch%n_ele_track) endif ! Set the kicker fields to zero to compute ring twiss parameters do i=1,branch%n_ele_max if(index(branch%ele(i)%name,'KICKER')/= 0) then if(branch%ele(i)%slave_status /= super_slave$)then call set_ele_real_attribute(branch%ele(i), 'KICKER_FIELD', 0.0_rp, err_flag) elseif(branch%ele(i)%slave_status == super_slave$)then lord_ele => pointer_to_lord(branch%ele(i),1) call set_ele_real_attribute(lord_ele, 'KICKER_FIELD', 0.0_rp, err_flag) endif endif end do 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 do i=1, branch%n_ele_track write(lun,'(es12.4,1x,6es12.4)')branch%ele(i)%s, co(i)%vec end do close(lun) scraping_on = save_scraping_status print *,' TRACKING_MASTER: call closed_orbit_calc, line 311' call closed_orbit_calc(lat, co, i_dim=4, ix_branch=nbranch) lun=lunget() open(unit=lun,file=trim(directory)//'/'//'closed_orbit.dat') write(lun,'(a12,1x,6a12)')' s[m]','x[m]','xp[rad]','y[m]','yp[rad]','z[m]','pz[dp/p]' s_start=0. delta_s = 0.01 co_start =co(0) do while(s_start+delta_s <= branch%ele(branch%n_ele_track)%s) call track_from_s_to_s(lat, s_start,s_start+delta_s,co_start, co_end, ix_branch = nbranch, all_orb = all_orb ) write(lun,'(es12.4,1x,6es12.4)')s_start+delta_s, co_end%vec s_start = s_start+delta_s co_start=co_end end do co_start = co(0) do i=1,branch%n_ele_track print '(i10,1x,a16,1x,3es12.4,1x,es12.4,1x,i10)',i, branch%ele(i)%name, branch%ele(i)%s,co(i)%vec(1),co(i)%vec(3), value_of_attribute(branch%ele(i), 'FIELD_INDEX'), branch%ele(i)%slave_status ! write(99,'(i10,1x,2es12.4,1x,2es12.4,es12.4)')i, co(i)%vec(1),co(i)%vec(3), all_orb(i)%vec(1),all_orb(i)%vec(3), branch%ele(i)%s call lat_make_mat6(lat,ix_ele = i,ref_orb = co, ix_branch=nbranch) end do close(unit=lun) if(index(ring_twiss,'close')/= 0) call twiss_at_start(lat,status, ix_branch = nbranch) !if initial twiss values are set in input file comment out this call if(status/=0)print '(a,a)',' Status from twiss_at_start = ', matrix_status_name(status) WRITE(*,'(a,2(3x,a,f8.5))') 'Fractional Tune: ', 'Qx =', lat%a%tune/twopi , 'Qy =', lat%b%tune/twopi WRITE(*,'(3(3x,a,f8.5))') 'beta_x =', lat%branch(1)%ele(0)%a%beta , 'beta_y =', lat%branch(1)%ele(0)%b%beta,'eta =', lat%branch(1)%ele(0)%x%eta ! rf_quad(1:4)%freq_h = 2* abs(twopi*rf_quad(1)%freq_h - lat%a%tune) / branch%ele(branch%n_ele_track)%s * c_light ! rf_quad(1:4)%freq_v = 2 * abs(twopi*rf_quad(1)%freq_v - lat%b%tune) / branch%ele(branch%n_ele_track)%s * c_light ! rf_quad(1:4)%freq_h = rf_quad(1:4)%freq_h *(twopi- lat%a%tune) / branch%ele(branch%n_ele_track)%s * c_light ! rf_quad(1:4)%freq_v = rf_quad(1:4)%freq_v * lat%b%tune / branch%ele(branch%n_ele_track)%s * c_light call twiss_propagate_all(lat, ix_branch=nbranch) call chrom_calc(lat, delta_e,chrom_x,chrom_y, err_flag, ix_branch=nbranch) print '(2(1x,a,es12.4))',' chrom_x = ', chrom_x,'chrom_y = ', chrom_y call write_log(lat, chrom_x, chrom_y) ! if(err_flag)print *,' Error from chrom calc' print *,' here' if (twiss_file) then call twiss_propagate_cm(lat,nbranch) call write_transfermatrix(lat, nbranch) endif matchlat = lat if(index(ring_twiss,'open')/=0)call MatchToRing(matchlat, initial_offsets_ref, kicker_params, twiss_match) !for many turn spin tracking, more accuracy is necessary ! bmad_com%min_ds_adaptive_tracking = 0.0 ! bmad_com%rel_tol_adaptive_tracking = 1d-11 ! Runge-Kutta tracking relative tolerance. ! bmad_com%abs_tol_adaptive_tracking = 1d-13 ! Runge-Kutta tracking absolute tolerance. rf_quad(1:4)%amp_h = rf_quad_0(1:4)%amp_h rf_quad(1:4)%amp_v = rf_quad_0(1:4)%amp_v if(any(rf_quad_0(1:4)%amp_h /= 0) .or. any(rf_quad_0(1:4)%amp_v /= 0)) then do i=1,4 print '(a13,es12.4,a13,es12.4,a,es12.4)',' rf_quad_h = ',rf_quad(i)%freq_h,' rf_quad_v = ', rf_quad(i)%freq_v,' rf_quad_amp_h = ',rf_quad_0(i)%amp_h end do endif call write_rfquad ix_end_flag = 1 vparam = 0 do while (vparam <= (vparam_max-vparam_min)) if(loop .and.vparam_id /= 0)then if(vparam_id == 1) twiss%betax = vparam +vparam_min if(vparam_id == 2) twiss%betay = vparam +vparam_min if(vparam_id == 3) twiss%etax = vparam +vparam_min if(vparam_id == 4) kicker_params%kicker_field(1:4) = kicker_params_0%kicker_field(1:4) * (vparam+vparam_min) ! all three kickers scaled by vparam if(vparam_id == 5) kicker_params%kick_width(1:4) = vparam +vparam_min if(vparam_id == 6) initial_offsets%x_mean = vparam +vparam_min if(vparam_id == 7) initial_offsets%pxmean = vparam +vparam_min if(vparam_id == 8) kicker_params%dtRise(1:4) = vparam +vparam_min if(vparam_id == 9) kicker_params%dtFall(1:4) = vparam +vparam_min if(vparam_id == 10)then inflector_field = vparam +vparam_min call set_inflector_params(lat, 0,inflector_angle) endif if(vparam_id == 11) initial_offsets%pzmean = vparam +vparam_min if(vparam_id == 12) twiss%etapx = vparam +vparam_min if(vparam_id == 13) twiss%alphax = vparam +vparam_min if(vparam_id == 14) kicker_tStart(1:3) = vparam +vparam_min if(vparam_id == 15) kicker_params%kicker_field(1) = vparam +vparam_min if(vparam_id == 16) kicker_params%kicker_field(2) = vparam +vparam_min if(vparam_id == 17) kicker_params%kicker_field(3) = vparam +vparam_min if(vparam_id == 18)then 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 do i=1,branch%n_ele_max if(index(branch%ele(i)%name,'KICKER')/= 0) then call set_ele_real_attribute(branch%ele(i), 'KICKER_FIELD', 0.0_rp, err_flag) endif end do rf_quad(1:4)%amp_h = 0 rf_quad(1:4)%amp_v = 0 co(0)%vec=0 call closed_orbit_calc(lat, co, i_dim=4, ix_branch=nbranch) do i=1,branch%n_ele_track print '(i10,1x,a16,1x,3es12.4,1x,es12.4,1x,i10)',i,branch%ele(i)%name, & branch%ele(i)%s,co(i)%vec(1),co(i)%vec(3), value_of_attribute(branch%ele(i), 'FIELD_INDEX'), branch%ele(i)%slave_status call lat_make_mat6(lat,ix_ele = i,ref_orb = co, ix_branch=nbranch) end do if(index(ring_twiss,'close')/= 0) call twiss_at_start(lat,status, ix_branch = nbranch) if(status/=0)print '(a,a)',' Status from twiss_at_start = ', matrix_status_name(status) WRITE(*,'(a,2(3x,a,f8.5))') 'Fractional Tune: ', 'Qx =', lat%a%tune/twopi , 'Qy =', lat%b%tune/twopi ! rf_quad(1:4)%freq_h = lat%a%tune * branch%ele(branch%n_ele_track)%s/c_light ! rf_quad(1:4)%freq_v = lat%b%tune * branch%ele(branch%n_ele_track)%s/c_light if((lat%a%tune == 0) .or. (lat%b%tune == 0))then vparam = vparam + delta_vparam cycle endif call twiss_propagate_all(lat, ix_branch=nbranch) ! rf_quad(1:4)%amp_h = rf_quad_0(1:4)%amp_h ! rf_quad(1:4)%amp_v = rf_quad_0(1:4)%amp_v endif if(vparam_id == 19)then inflector_angle = vparam +vparam_min call set_inflector_params(lat, 0,inflector_angle) endif if(vparam_id >= 20 .and.vparam_id <=26)then !rf quad params if(vparam_id == 20)rf_quad(1:4)%amp_h= vparam + vparam_min if(vparam_id == 21)rf_quad(1:4)%amp_v= vparam + vparam_min if(vparam_id == 22)rf_quad(1:4)%phi_h= rf_quad_0(1:4)%phi_h + vparam_min + vparam ! all four quads scale with the same delta phi if(vparam_id == 23)rf_quad(1:4)%phi_v= rf_quad_0(1:4)%phi_v + vparam_min + vparam ! all four quads scale with the same delta phi if(vparam_id == 24)rf_quad(1:4)%start_h= vparam_min + vparam if(vparam_id == 25)rf_quad(1:4)%start_v= vparam_min + vparam if(vparam_id == 26)rf_quad(1:4)%freq_h= vparam_min + vparam endif xmax = 0. ymax = 0. endif print *, ' BETA x =', twiss%betax print '(a,es12.4,a,es12.4)', 'kicker_field=', kicker_params%kicker_field(1), ' kicker_rise = ',kicker_params%dtRise(1) ! Set kicker parameters (now that ring's Twiss parameters, etc., have been calculated) 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. print '(a,es16.8,a,es16.8)',' circumference = ',circumference,' radius = ',circumference/twopi if(abs(circumference/twopi - 7.112)>1.e-5) then print *,' circumference inconsistent with radius' stop endif circumference2 = 0. !test variable ! print*, circumference ! branch => lat%branch(1) call reallocate_coord (from_orbit, from_branch%n_ele_max) ! call init_coord(from_orbit) ! allocate(co_electron(nmuons)) call reallocate_coord (co_electron, nmuons) call reallocate_coord (co, branch%n_ele_max) call reallocate_coord (all_orb, branch%n_ele_max) ! call reallocate_coord (coSteps, steps) ! call init_coord(co(0)) ! call init_coord(coSteps(0)) if(.not. initial_offset_from_co)co_start%vec=0 start_orb%vec = 0 start_orb%t = initial_offsets%tmean start_orb%s = 0 start_orb%vec(1) = initial_offsets%x_mean +co_start%vec(1) start_orb%vec(2) = initial_offsets%pxmean +co_start%vec(2) start_orb%vec(3) = initial_offsets%y_mean +co_start%vec(3) start_orb%vec(4) = initial_offsets%pymean +co_start%vec(4) start_orb%vec(5) = initial_offsets%tmean +co_start%vec(5) start_orb%vec(6) = initial_offsets%pzmean +co_start%vec(6) co(0) = start_orb print '(a,6es12.4)','start_orb',co(0)%vec co(0)%t = initial_offsets_ref%tmean co(0)%s = 0 co(0)%vec(1) = initial_offsets_ref%x_mean co(0)%vec(2) = initial_offsets_ref%pxmean co(0)%vec(3) = initial_offsets_ref%y_mean co(0)%vec(4) = initial_offsets_ref%pymean co(0)%vec(5) = initial_offsets_ref%tmean co(0)%vec(6) = initial_offsets_ref%pzmean print '(a,6es12.4)','start_orb ref',co(0)%vec lat%ele(1:lat%n_ele_track)%alias= 'Record' !write magnetic field and position along trajectory to unit 12 in em_field_custom_inj.f90 call track_all(lat,co, ix_branch = 0, track_state = track_state, err_flag = err_flag) !compute trajectory through injection line lat%ele(1:lat%n_ele_track)%alias= 'No' if(opt_incident)call optimize_incident(start_orb, lat, inflector_end_target, opt_orb, inflector_angle) !routine to find incident offset and angle to minimize offset and angle at inflector exit print *,' trajectory through backleg and inflector' call track_write_backleg_trajectory(lat,co) call lat_make_mat6(lat, -1, co, ix_branch=0) ! Create transfer matrices for injection line elements referred to trajectory if(twiss_file)call write_transfermatrix(lat, 0) if(start_tracking_at_inflector_exit) then call mat_make_unit(mat_inv) else call get_branch_inv_matrix(lat,mat_inv,mat, inj_branch, err) ! Multiply matrices for injection line together then compute inverse to pass to CREATE_PHASE_SPACE if(.not. use_lattice_twiss .and. twiss%betax /= 0 .and. (.not. err))then !use twiss parameters propagated backwards from inflector call prop_phase_space(mat_inv,twiss, twiss_start_inj_line) call copy_twiss_to_lat(twiss_start_inj_line,lat,0) elseif (lat%branch(0)%ele(0)%a%beta /= 0)then ! use twiss parameters in lattice file call copy_lat_to_twiss(lat,twiss,0) ! copy twiss parameters specified in lattice file to twiss that are used to create distributon call mat_make_unit(mat_inv) elseif(lat%branch(0)%ele(0)%a%beta == 0)then print '(a)',' BETA ZERO at START. INITIALIZE in LATTICE file or in INPUT.DAT' ! stop endif endif call create_phase_space(lat, nmuons, create_new_distribution, new_file, muon_file, & epsdistr, epsx, epsy, tdistr, tlength, tsigma, pzdistr, pz, pzsigma, twiss,muons, twiss_ref, mat, mat_inv, nmuon_first ) !mat_inv propagates from end of inflector back to start of injection line ! call compute_beam_params(muons, 1, 'PropagatedBackToStart') do i=1, nmuons ! call init_coord(muons(i)%coord, muons(i)%coord%vec, lat%ele(0), element_end = upstream_end$, & ! particle=lat%param%particle, t_offset =muons(i)%coord%t ,spin = muons(i)%coord%spin ) !remove t_offset from call 10/17/2018 t_save = muons(i)%coord%t r_save = muons(i)%coord%r call init_coord(muons(i)%coord, muons(i)%coord%vec, lat%ele(0), element_end = upstream_end$, & particle=lat%param%particle, spin = muons(i)%coord%spin ) muons(i)%coord%t=t_save muons(i)%coord%r=r_save end do !turn on scattering now that reference orbit and transfer matrices through injection line are computed eloss = enerloss us_scatter = inf_end_us ds_scatter = inf_end_ds FiberScatter = FiberScatter0 FiberEnergyLoss = FiberEnergyLoss0 n_row = 1 if(save_element_by_element_info)allocate(muons_ele(0:nmuons, 0:branch%n_ele_track)) allocate(muons_ele_inj(0:nmuons, 0:lat%branch(0)%n_ele_track)) allocate(muons_temp(0:nmuons)) allocate(electrons_s(0:n_row,0:esteps*(tot_track+1))) from_orbit(0)%t=0. from_orbit(0)%s=0. from_orbit(0)%vec=0. forall (n=1:nmuons) muons(n)%coord%vec(1:6) = muons(n)%coord%vec(1:6) + from_orbit(0)%vec(1:6) muons(:)%coord%t = muons(:)%coord%t + from_orbit(0)%t muons(:)%coord%s = muons(:)%coord%s + from_orbit(0)%s ! muons(:)%lost = .false. ! call write_phase_space(nmuons,muons,'', tot_inf_scatter) call init_moment(moment) call init_moment(moment_ele) call compute_moments(muons,0.d0,moment, nmuons,nmu, 0, 'BEGINNING',0) call compute_moments(muons,-lat%branch(0)%ele(lat%branch(0)%n_ele_track)%s/circumference,moment, nmuons,nmu, 1,'BEGINNING',0) ! write header lines to average_spin_tbt.dat and average_spin_by_element.dat if (spin_tracking_on) then call compute_spins(muons(1:nmuons),0.d0,nmuons,0,'BEGINNING',0,0.d0) call compute_spins(muons(1:nmuons),-lat%branch(0)%ele(lat%branch(0)%n_ele_track)%s/circumference, & nmuons,1,'BEGINNING',0,-lat%branch(0)%ele(lat%branch(0)%n_ele_track)%s) endif end=lat%n_ele_track ibms=.true. i=0 if(nmuons == 1)then call write_single_particle_tbt(i,muons,Q,nturns, lat%ele(end)%ref_time) call write_single_particle_by_element(i,lat%branch(0)%ele(0)%name,muons(1)%coord,lat%ele(0)%s/circumference, & ((i-1)*lat%ele(lat%n_ele_track)%s +lat%ele(0)%s)/circumference,Q) call write_injection_line_trajectory(lat%branch(0)%ele(0)%s,start_orb%vec(1:6), start_orb%vec(1:6), lat%branch(0)%ele(0)%name, 1, 0.,0.,0.) call write_injection_line_twiss(lat%branch(0)%ele(0)) xmax = max(xmax, abs(muons(1)%coord%vec(1))) ymax = max(ymax, abs(muons(1)%coord%vec(3))) endif !to write kicker pulse height vs time along with number of muons vs time set alias to 'trajectory' ! branch%ele(1:branch%n_ele_track)%alias= 'trajectory' !lat%branch(0)%ele(1:lat%branch(0)%n_ele_track)%alias = 'trajectory' branch%ele(1:branch%n_ele_track)%alias= '' lat%branch(0)%ele(1:lat%branch(0)%n_ele_track)%alias = '' em_field_calc_verbose= em_field_calc_verbose0 DO i=1, nturns if(i == 1 ) then !propagate through injection line to start of ring DO n=1, nmuons !propagate all muons through branch 0 if(allocated(muons_ele))muons_ele(n,:)%coord%state = lost$ ! call init_coord(muons(n)%coord, particle = lat%param%particle) IF (muons(n)%coord%state /= alive$) CYCLE if(start_tracking_at_inflector_exit)then ! print *,' line 630' do j=1,lat%branch(0)%n_ele_track from_orbit(j) = muons(n)%coord from_orbit(j)%vec = muons(n)%coord%vec + start_orb%vec from_orbit(j)%t = muons(n)%coord%t !+ lat%branch(0)%ele(j)%s/c_light muons_ele_inj(n,j)%coord = from_orbit(j) ! print '(i10,1x,6es12.4)',j,from_orbit(j)%vec end do else from_orbit(0) = muons(n)%coord from_orbit(0)%vec = muons(n)%coord%vec+ start_orb%vec from_orbit(0)%t = muons(n)%coord%t from_orbit(0)%spin = muons(n)%coord%spin do j=1,lat%branch(0)%n_ele_track ele_save%tracking_method = lat%branch(0)%ele(j)%tracking_method if(inj_matrix_tracking .and.lat%branch(0)%ele(j)%key /= marker$)lat%branch(0)%ele(j)%tracking_method = linear$ call track1(from_orbit(j-1),lat%branch(0)%ele(j), lat%param, from_orbit(j), err_flag =err) if (err .or. .not. particle_is_moving_forward(from_orbit(j))) then call set_orbit_to_zero (from_orbit, j+1, lat%branch(0)%n_ele_track) if (from_orbit(j)%location == upstream_end$) from_orbit(j)%vec = 0 ! But do not reset orbit(n)%state exit endif lat%branch(0)%ele(j)%tracking_method = ele_save%tracking_method muons_ele_inj(n,j)%coord = from_orbit(j) end do endif IF (nmuons==1 .and. twiss_file)then ! Use injection channel with lattice exteneded 10cm upstream of IBMS1 - twiss measuemernt reference. call first_turn(lat,spin_tracking_on, from_orbit, circumference, make_movie, start_tracking_at_inflector_exit) ! call injection_channel(lat,spin_tracking_on, from_orbit, circumference, make_movie) endif rfquad_params_reset=.true. !passed in parameters_bmad reset after first turn ! print '(a,1x,i10,1x,6es12.4)',' from_orbit ',j, from_orbit(j)%vec to_orbit = from_orbit(ie_from) muons(n)%coord%vec = to_orbit%vec muons(n)%coord%t = to_orbit%t muons(n)%coord = to_orbit ! print *,' line 751' ! print '(i10,1x,6es12.4)',j,to_orbit%vec end do ! end of loop propagating all muons through branch 0 ! generate Decay time(lifetime) if(muon_decay_on)call compute_decayTime(muons(1:nmuons),nmuons,11) if(write_phase_space_file)call write_phase_space( nmuons, muons,lat%branch(0)%ele(lat%branch(0)%n_ele_track)%name, tot_inf_scatter,turn) endif !for first turn only do n=1, nmuons muon_number = n if(allocated(muons_ele))muons_ele(n,:)%coord%state = lost$ !write(6,'(a,1x,i10,2x,6es12.4,1x,i10,1x,a)')'tracking_master: muon = ',n, muons(n)%coord%vec, track_state,coord_state_name(muons(n)%coord)%state) if(muons(n)%coord%r < BetaCrossE_start_time)muons(n)%coord%state = lost$ IF (muons(n)%coord%state /= alive$) CYCLE co(0)%vec = muons(n)%coord%vec co(0)%t = muons(n)%coord%t co(0) = muons(n)%coord co(0)%s=0. if(int_betaCrossE)then call track_all_int_efield(lat, co , BetaCrossE_start_time, sumBetaCrossE, nbranch, track_state, err_flag, verbose = betaCrossE_verbose) else if(opt_inf_exit_angle .and. nmuons == 1 )then print '(a,7es12.4)','co(0)%vec, co(0)%t',co(0)%vec, co(0)%t call opt_inf_angle(lat, co,nbranch, track_state, err_flag, kicker_params) opt_inf_exit_angle=.false. endif call track_all(lat, co, nbranch, track_state, err_flag) endif if(track_state >= 0)write(99,'(a,1x,i10,1x,i10,1x,a)')'tracking_master: muon = ',n, track_state,coord_state_name(co(track_state)%state) if(nmuons == 1)then do j=1,branch%n_ele_track s_start=0. if(make_movie)call step_floor_around_ring(lat,j,nbranch,co(j-1) ,branch%ele(j-1)%floor%r, s_start) end do endif DO j=1, branch%n_ele_track if(allocated(muons_ele))muons_ele(n,j)%coord = co(j) IF (nmuons == 1)then if(spin_tracking_on) Q = co(j)%spin call write_single_particle_by_element(i,branch%ele(j)%name,co(j), branch%ele(j)%s/circumference, & ((i-1)*branch%ele(branch%n_ele_track)%s +branch%ele(j)%s)/circumference,Q) xmax = max(xmax, abs(muons(1)%coord%vec(1))) ymax = max(ymax, abs(muons(1)%coord%vec(3))) endif if (err_flag .and. allocated(muons_ele)) muons_ele(n,j)%lost = .true. if (j == track_state .and. track_state /= moving_forward$)then call write_lost_muons(n,branch%ele(j)%s, co(j), branch%ele(j)%name,i) do k=j,branch%n_ele_track if(allocated(muons_ele))muons_ele(n,k)%coord%state = lost$ enddo endif ENDDO ! track through lattice elements if(track_state == moving_forward$ .and. co(1)%t > time_binning_start & .and.bin_width > 0 .and. pbin_width >0 .and. trim(time_binning) == 'momentum_slice')call time_bins(co, branch%n_ele_track, nmuons, time_binning_start, ele_time_bins) if (track_state /= moving_forward$ .or. err_flag) then muons(n)%coord%state = co(track_state)%state muons(n)%lost = .true. cycle endif !*********************************************** !see if muon has decayed. If so,retrack to see where it decayed and obtain the corresponding electron phase space vector if(muon_decay_on) then muon_lifetime = co(0)%r if (co(branch%n_ele_track)%t >= muon_lifetime .and. co(0)%t < muon_lifetime) then do j = 1, branch%n_ele_track if (co(j)%t >= muon_lifetime) then track%n_pt = -1 call track1(co(j-1),branch%ele(j),lat%param,co(j), track=track, err_flag=err_flag, ignore_radiation=.true.) if(allocated(muons_ele))muons_ele(n,j:branch%n_ele_track)%coord%state = lost$ ! if(track%n_pt <= 0 .or. track%pt(0)%orb%t > muon_lifetime)pause ix=bracket_index(muon_lifetime, track%pt(0:track%n_pt)%orb%t, imin) !bracket index where muon decayed ! if(ix < 0 )pause co_muon_end = track%pt(ix)%orb if(ix < track%n_pt)then !interpolate dt=track%pt(ix+1)%orb%t- track%pt(ix)%orb%t dt1 = muon_lifetime - track%pt(ix)%orb%t dt2 = dt-dt1 if(dt > 0)then co_muon_end%s = (dt2*track%pt(ix)%orb%s+dt1*track%pt(ix+1)%orb%s)/dt co_muon_end%t = (dt2*track%pt(ix)%orb%t+dt1*track%pt(ix+1)%orb%t)/dt co_muon_end%vec = (dt2*track%pt(ix)%orb%vec+dt1*track%pt(ix+1)%orb%vec)/dt else co_muon_end%s = track%pt(ix)%orb%s co_muon_end%t = track%pt(ix)%orb%t co_muon_end%vec = track%pt(ix)%orb%vec endif endif ! save info about decay muon muons(n)%coord = co_muon_end muons(n)%coord%state = lost$ muons(n)%lost = .true. muons(n)%decay = .true. muons(n)%turn = i-1 !muon has not completed turn i muons(n)%path_1 = float(i-2) + co_muon_end%s/circumference if(allocated(sumBetaCrossE))then muons(n)%betaCrossE = muons(n)%betaCrossE + sumBetaCrossE(branch%n_ele_track)%vec muons(n)%betaCrossE_time = muons(n)%betaCrossE_time + sumBetaCrossE(branch%n_ele_track)%dt muons(n)%betaDotB(1:3) = muons(n)%betaDotB(1:3) + sumBetaCrossE(branch%n_ele_track)%vec_pitch(1:3) muons(n)%path = muons(n)%path + sumBetaCrossE(branch%n_ele_track)%path muons(n)%pathx = muons(n)%pathx + sumBetaCrossE(branch%n_ele_track)%pathx muons(n)%B_field = muons(n)%B_field + sumBetaCrossE(branch%n_ele_track)%B_field muons(n)%BetaCrossB = muons(n)%BetaCrossB + sumBetaCrossE(branch%n_ele_track)%BetaCrossB endif distance = float(i-1)*circumference + co_muon_end%s if(.not. spin_tracking_on)call decay_info(distance,n,co_electron(n),co_muon_end, muons = muons(n)) if(spin_tracking_on)call create_and_propagate_electron(lat,nbranch, co_muon_end, co_electron(n),n, i, nmuons, make_movie) exit ! the muon decayed in element j so now exit loop over elements end if ! if decay in element j end do ! element search j=1,branch%n_ele_track cycle !go to next muon end if ! if decay at the end if (co(branch%n_ele_track)%t >= co(0)%path_len) then ! if(co(0)%t < muon_lifetime)then ! muons(n)%coord%state = lost$ ! muons(n)%lost = .true. ! endif endif !muon_decay_on !********************************************** muons(n)%coord%vec(1:6) = co(branch%n_ele_track)%vec(1:6) muons(n)%coord%t = co(branch%n_ele_track)%t muons(n)%coord = co(branch%n_ele_track) muons(n)%turn = i muons(n)%path_1 = float(i-1) if(allocated(sumBetaCrossE))then muons(n)%betaCrossE = muons(n)%betaCrossE + sumBetaCrossE(branch%n_ele_track)%vec muons(n)%betaCrossE_time = muons(n)%betaCrossE_time + sumBetaCrossE(branch%n_ele_track)%dt muons(n)%betaDotB(1:3) = muons(n)%betaDotB(1:3) + sumBetaCrossE(branch%n_ele_track)%vec_pitch(1:3) muons(n)%path = muons(n)%path + sumBetaCrossE(branch%n_ele_track)%path muons(n)%pathx = muons(n)%pathx + sumBetaCrossE(branch%n_ele_track)%pathx muons(n)%B_field = muons(n)%B_field + sumBetaCrossE(branch%n_ele_track)%B_field muons(n)%BetaCrossB = muons(n)%BetaCrossB + sumBetaCrossE(branch%n_ele_track)%BetaCrossB if(nmuons == 1)call write_BetaCrossE_turn(i, co(branch%n_ele_track), sumBetaCrossE(branch%n_ele_track), muons(n), BetaCrossE_start_time) !%betaCrossE(1:3), muons(n)%betaCrossE_time, muons(n)%betaDotB(1:3), muons(n)%path, muons(n)%pathx) endif ENDDO ! n muons if(i == 1)then !first turn - compute moments and spins, save phase space, through injection line do j=1,lat%branch(0)%n_ele_track distance = float(i-1)-(lat%branch(0)%ele(lat%branch(0)%n_ele_track)%s-lat%branch(0)%ele(j)%s)/circumference if(j < lat%branch(0)%n_ele_track .and. start_tracking_at_inflector_exit)cycle !if tracking starts at inf exit, there is data only for the last element in the inj channel forall(n=1:nmuons)muons_temp(n) = muons_ele_inj(n,j) call compute_moments(muons_temp,distance,moment_ele, nmuons,nmu, 1, lat%branch(0)%ele(j)%name, j) if(write_phase_space_file)call write_phase_space(nmuons,muons_temp, lat%branch(0)%ele(j)%name, tot, turn) !write phase space after each element in injection line if(write_phase_space_every_n_turns >0 .and. j == lat%branch(0)%n_ele_track) & call write_phase_space(nmuons,muons_temp, trim(lat%branch(0)%ele(j)%name)//'000', tot, turn) !write phase space at end of injection line if (spin_tracking_on) call compute_spins(muons_temp(1:nmuons),distance,nmuons,1, lat%branch(0)%ele(j)%name, j,lat%branch(0)%ele(j)%s) !compute and write average spins for element call moment_range(moment_ele) end do endif do j=1,branch%n_ele_track distance = float(i-1)+branch%ele(j)%s/circumference if(allocated(muons_ele))then forall(n=1:nmuons)muons_temp(n) = muons_ele(n,j) else if(j /= branch%n_ele_track)cycle forall(n=1:nmuons)muons_temp(n) = muons(n) endif call compute_moments(muons_temp,distance,moment_ele, nmuons,nmu, 1, branch%ele(j)%name, j) if(i == 1)then if(write_phase_space_file .or. index(branch%ele(j)%name,'END')/=0) & call write_phase_space(nmuons,muons_temp, branch%ele(j)%name, tot, i) !write phase space after each element in first turn endif if(write_phase_space_every_n_turns > 0 .and. j == branch%n_ele_track)then ndum = write_phase_space_every_n_turns char_turn='000' if(i < 10)write(char_turn(3:3),'(i1)')i if(i >= 10 .and. i < 100)write(char_turn(2:3),'(i2)')i if(i >= 100 .and. i < 1000)write(char_turn(1:3),'(i3)')i if( (i/ndum) * ndum == i) call write_phase_space(nmuons,muons_temp, trim(branch%ele(j)%name)//char_turn, tot, i) !write phase space after each element in first turn endif if (spin_tracking_on) call compute_spins(muons_temp(1:nmuons),distance,nmuons,1, branch%ele(j)%name, j,branch%ele(j)%s) !compute and write average spins for elements if(nturns - i < vparam_cbo_turns)call moment_range(moment_ele) end do distance = float(i) call compute_moments(muons,distance,moment, nmuons,nmu, 0, branch%ele(branch%n_ele_track)%name, j-1) !j somehow becomes n_ele_track+1 at the end of the loop if(time_bin_width > 0)call compute_moments_vs_time(nturns, i,nmuons, muons, branch%ele(branch%n_ele_track)%name) if (spin_tracking_on) call compute_spins(muons(1:nmuons),distance,nmuons,0,branch%ele(branch%n_ele_track)%name, j-1,branch%ele(j-1)%s) !compute and write average spins if(nturns - i < vparam_cbo_turns)then call moment_range(moment) !compute range of moments only for last 100 turns if(100-(nturns-i)==1)write(154,'(a12,a10,3a12)')'vparam','turn','avex','max_avex','min_avex' write(154,'(es12.4,1x,i10,1x,3es12.4)')vparam,i,moment%ave(1),moment%max_ave(1),moment%min_ave(1) endif if(i > 10)call collect_times(muons, nturns, NN, bins) !start to collect times of particles that are likely to survive if(nmuons == 1)then if(spin_tracking_on) Q = muons(1)%coord%spin call write_single_particle_tbt(i,muons,Q, nturns, branch%ele(branch%n_ele_track)%ref_time) endif ENDDO ! nturns if(bin_width > 0 .and. pbin_width >0 .and. trim(time_binning) == 'momentum_slice')call time_bins(co, -1, nmuons, time_binning_start, ele_time_bins) if(muon_decay_on .and. spin_tracking_on)call electron_info_s(electrons_s,n_row,esteps*(1 + tot_track),n,.false.) !sample one electon's trajectory if(time_bin_width > 0)then call compute_moments_vs_time(nturns, -1,nmuons,muons,branch%ele(branch%n_ele_track)%name, coef_save=coef,chisq_save=chisq,ndf_save=ndf) call compute_moments_vs_time_calo(nturns, -1,co(branch%n_ele_track), branch%ele(branch%n_ele_track)%name) endif ! lun = lunget() ! open(unit = lun, file = 'muons_vs_time.dat') ! print '(2a)',' write binned time data to ','muons_vs_time.dat' ! print *,' bins = ', bins ! write(lun,'(a10,a12)')'Bin','NN(bin)' ! write(lun,'(i10,es12.4)')(bin, NN(bin), bin=1, bins) ! close(lun) call write_phase_space(nmuons,muons,'EndOfTracking', tot,nturns, BetaCrossE_start_time) print '(/,5a,/)',' Use plotting script ',"'",'sigma.gnu',"'",' to plot turn-by-turn average position and size' if (loop)then call write_loop (vparam_id, vparam_label(vparam_id), vparam_column_index(vparam_id), nmuons,kicker_params, moment, tot_inf_scatter, tot, twiss, initial_offsets, xmax, ymax, inflector_angle, lat) if(time_bin_width > 0) call write_rfquad_loop(vparam_id, vparam_label(vparam_id), vparam_column_index(vparam_id), nmuons, kicker_params, tot, coef,chisq,ndf) if(tot_max < tot)then tot_max = tot vparam_max_tot = vparam + vparam_min endif if(cbo_min >( moment%max_ave(1)-moment%min_ave(1))*1.e3)then cbo_min = (moment%max_ave(1)-moment%min_ave(1))*1.e3 vparam_cbo_min = vparam + vparam_min endif call create_gnu_input(tot_max,vparam_max_tot, cbo_min,vparam_cbo_min) endif !if looping for vparam deallocate(muons) if(allocated(muons_ele))deallocate(muons_ele) deallocate(muons_temp) deallocate(muons_ele_inj) deallocate(electrons_s) vparam = vparam + delta_vparam if(.not. loop .or. vparam_id <= 0) exit end do lun = lunclose() print '(a,i10,a)',' Close ',lun,' units' print '(a,a)',' output written to directory = ', trim(directory) call cpu_time(finish) print '(a,f10.2,a)',' CPU time = ',(finish-start)/60.,' minutes' end