program tracking_master use bmad use muon_mod use muon_interface use parameters_bmad use runge_kutta_mod use materials_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 (coord_struct), allocatable :: from_orbit(:) type (coord_struct), allocatable, save :: co(:),co2(:) !co2 for closed orbit calculation type (coord_struct), allocatable :: coSteps(:),co_electron(:), co_elec(:) type (coord_struct), allocatable :: low_e_orb(:), high_e_orb(:) type (coord_struct) to_orbit type (coord_struct) start_orb, orb_at_s type (coord_struct) opt_orb type (coord_struct) co_muon_end type (coord_struct) co_elec_end type (all_pointer_struct) a_ptr type (ele_pointer_struct), allocatable :: eles(:) type (em_field_struct) field type (track_struct) track type (ele_struct) ele integer pgopen, istat1, istat2 integer i, j, k,l,steps/1000/,l1 !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 lun,lun32 integer nargs, iargc, ios integer end, nturns, nmuons integer nmu integer bin,bins/0/ integer n integer track_state integer unit integer turn 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/ integer ix_end integer seed/1/ integer i2,i3 integer lunclose 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) offset_m !misalignment of the quadrupole logical first_offset/.true./ logical err integer n_loc integer tot_max/0/ character*10 offset_cm !misalignment of the quadrupole character*120 line, lat_file, lat_file_name/''/,new_file/' '/, muon_file/' '/ character*16 ele_name character*20 vparam_label(18)/'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]' / integer vparam_column_index(18)/17,20,18,19,21,24,25,22,23,26,27,28,29,30,1,2,3,31/ logical err_flag, inf_end_us/.false./, inf_end_ds/.false./, quad_plate/.false./ logical first/.true./, twiss_file/.false./,first_file/.true./,first_file2/.true./ 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 record_fiber_data/.true./ logical use_lattice_twiss/.false./ logical write_phase_space_file/.false./ logical save_element_by_element_info/.true./ logical cl_exist 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 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 character*16 epsdistr, tdistr, pzdistr, inf_aperture, twiss_ref, ring_twiss ! Input namelist structure ("/g-2/test/input.dat") namelist /input/ lat_file_name, nmuons, nturns, & ! GENERAL 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, enerloss, inflector_angle, inflector_field, & ! INFLECTOR ring_theta, & ! ELEMENT POSITIONING kickerPlates, kickerCircuit, 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, & 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, & record_fiber_data, & ! record fiber harp data 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 ! 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 kicker_params%dtRise = 0. kicker_params%dtFall = 0. quad_params%short_quad_field_index = -99. quad_params%long_quad_field_index = -99. rf_quad(1:4)%amp_h=0 rf_quad(1:4)%amp_v=0 OPEN (UNIT=5, FILE='input.dat', STATUS='old', IOSTAT=ios) READ (5, NML=input, IOSTAT=ios) WRITE(6,NML=input) print *, 'ios=', ios rewind(unit=5) ! READ (5, NML=input) CLOSE(5) 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 tilt = inflector_angle ! save initial kicks kicker_params_0%kicker_field(1:3) = kicker_params%kicker_field(1:3) ! 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) 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(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 = cesr_iargc() if (nargs == 1) then call cesr_getarg(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 ! 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) 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. ! Scattering in inflector etc. is off while parser computes transfer matrices call initializeMaterials() call bmad_parser (lat_file, lat) do l1 =1, lat%branch(1)%n_ele_track if (lat%branch(1)%ele(l1)%name == 'QUAD1_SHORT') then print *, "x_offset-------------------",lat%branch(1)%ele(l1)%value(x_offset$) exit end if end do do l1 =1, lat%branch(1)%n_ele_track if (lat%branch(1)%ele(l1)%name == 'QUAD1_SHORT') then lat%branch(1)%ele(l1)%value(x_offset$) = 0.00 print *, "x_offset-------------------",lat%branch(1)%ele(l1)%value(x_offset$) exit end if end do call write_log(lat) ! lat%ele(1:lat%n_ele_track)%alias= 'No' ! nbranch = ubound(lat%branch,1) !refers to the ring branch nbranch = 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 eloss = .false. us_scatter = .false. ds_scatter = .false. quad_p = quad_plate quad_plate = .false. ! if(index(inflector_ape,'989')/=0)inflector_width=0.018 lun=lunget() open(unit=lun,file='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) 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 ring ! These values are set in the input file branch%ele(0)%a%beta = twiss%betax branch%ele(0)%b%beta = 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_quad_params(lat, nbranch) ! 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 branch%ele(i)%value(kicker_field$) = 0. endif end do call sleep(2) !set quadruple x_offset and loop over it. We use vparam_id > 100 to distinguish it from other vparam_id. Also there is a different do loop. vparam = 0 do while (vparam <= (vparam_max-vparam_min)+0.1*delta_vparam) if(loop .and. vparam_id == 100)then offset_m = vparam_min + vparam ! if(abs(offset_m) < 10.**(-10)) offset_m = 0. ! print *, (abs(offset_m) < 10.**(-10)) ! print *, "abs_offset1", abs(offset_m) call lat_ele_locator('QUAD1_SHORT',lat,eles,n_loc,err) if(n_loc == 0 .or. err)then print *,' Cannot find element ','QUAD1_SHORT', ' Try again' endif call pointer_to_attribute(eles(1)%ele, 'X_OFFSET', .false., a_ptr, err_flag, .false.) ! if(ix_attr .eq. 0)then if(err_flag)then call type_ele (eles(1)%ele, .true., 0, .false.,0, .false.) print *,' Attributes of ',ele%name,' shown above. Try again' endif a_ptr%r = offset_m print *, 'a_ptr__________________________________________________',a_ptr%r call set_flags_for_changed_attribute (eles(1)%ele, eles(1)%ele%value(x_offset$)) call lattice_bookkeeper(lat, err_flag) call lat_make_mat6(lat, -1) do l1 =1,branch%n_ele_track if (branch%ele(l1)%name == 'QUAD1_SHORT') then print*, 'x_offset__________________', branch%ele(l1)%value(x_offset$) ! call lattice_bookkeeper(lat, err_flag) exit end if end do end if !if (first_offset) then ! this if statement is introduced so that the following codes are run only once inside the loop of different element offsets !first_offset = .false. !end if !if first_offset !print *, "vparam_min", vparam_min !print *, "vparam_max", vparam_max !print *, "delta_vparam",delta_vparam print *, 'offset', offset_m do i=1,branch%n_ele_track if (branch%ele(i)%name == 'QUAD1_SHORT') then write (offset_cm,'(f5.2)') (vparam_min+vparam)*100 exit end if end do co = co2 !initialization before closed orbit calculation ! Calculate the ring's Twiss parameters, coupling, and dispersion parameters call closed_orbit_calc(lat, co, i_dim=4, ix_branch=nbranch) ! if (first_file2) then lun = lunget() open(unit=lun,file= 'closed_orbit_shift'//trim(adjustl(offset_cm))//'cm.dat', status = 'replace') close(lun) ! first_file2 = .false. ! end if do i=1,branch%n_ele_track lun = lunget() open(unit=lun,file= 'closed_orbit_shift'//trim(adjustl(offset_cm))//'cm.dat', access = 'append') write(lun,'(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), branch%ele(i)%value(field_index$), branch%ele(i)%slave_status close(lun) if (branch%ele(i)%name == 'QUAD1_SHORT') then ! inquire(file='close_orbit_vs_shift.dat',exist= cl_exist) if (.not.first_file ) then lun = lunget() open(unit=lun,file= 'closed_orbit_vs_shift.dat',status = 'old',access = 'append') write(lun,'(i10,1x,a16,1x,3es12.4,1x,es12.4,1x,i10,es12.4)'),i,branch%ele(i)%name,branch%ele(i)%s,co(i)%vec(1),co(i)%vec(3), branch%ele(i)%value(field_index$), branch%ele(i)%slave_status, vparam_min + vparam close(lun) ! first_file = .false. ! end if else lun = lunget() open(unit=lun,file= 'closed_orbit_vs_shift.dat',status = 'replace') write(lun,'(i10,1x,a16,1x,3es12.4,1x,es12.4,1x,i10,es12.4)'),i,branch%ele(i)%name,branch%ele(i)%s,co(i)%vec(1),co(i)%vec(3), branch%ele(i)%value(field_index$), branch%ele(i)%slave_status, vparam_min + vparam close(lun) first_file = .false. end if end if print '(i10,1x,a16,1x,3es12.4,1x,es12.4,1x,i10,a10)',i,branch%ele(i)%name,branch%ele(i)%s,co(i)%vec(1),co(i)%vec(3), branch%ele(i)%value(field_index$), branch%ele(i)%slave_status ! call lat_make_mat6(lat,ix_ele = i,ref_orb = co, ix_branch=nbranch) end do if (vparam_id >= 100) vparam = vparam + delta_vparam if(.not. loop .or. vparam_id < 100 ) exit end do if (vparam_id >=100) stop('fin') 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 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 if(any(rf_quad_0(1:4)%amp_h /= 0) .or. any(rf_quad_0(1:4)%amp_v /= 0)) & print '(a13,es12.4,a13,es12.4,a,es12.4)',' rf_quad_h = ',rf_quad(1)%freq_h,' rf_quad_v = ', rf_quad(1)%freq_v,' rf_quad_amp_h = ',rf_quad(1)%amp_h 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 ! 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 call MatchToRing(matchlat, initial_offsets, 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 call write_rfquad if (vparam_id <100) 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:3) = kicker_params_0%kicker_field(1:3) * (vparam+vparam_min) ! all three kickers scaled by vparam if(vparam_id == 5) kicker_params%kick_width(1:3) = 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:3) = vparam +vparam_min if(vparam_id == 9) kicker_params%dtFall(1:3) = vparam +vparam_min if(vparam_id == 10) inflector_field = vparam +vparam_min 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 quad_params%short_quad_field_index(1:4) = quad_params_0%short_quad_field_index(1:4) * (vparam+vparam_min) !all short quads quad_params%long_quad_field_index(1:4) = quad_params_0%long_quad_field_index(1:4) * (vparam+vparam_min) !all short quads 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 branch%ele(i)%value(kicker_field$) = 0. 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), branch%ele(i)%value(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 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 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. circumference = branch%ele(branch%n_ele_track)%s 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 (coSteps, steps) call reallocate_coord (co_elec, branch%n_ele_max) call init_coord(co(0)) call init_coord(coSteps(0)) start_orb%vec = 0 start_orb%t = initial_offsets%tmean start_orb%s = 0 start_orb%vec(1) = initial_offsets%x_mean start_orb%vec(2) = initial_offsets%pxmean start_orb%vec(3) = initial_offsets%y_mean start_orb%vec(4) = initial_offsets%pymean start_orb%vec(5) = initial_offsets%tmean start_orb%vec(6) = initial_offsets%pzmean co(0) = start_orb print '(a,6es12.4)','start_orb',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, opt_orb) !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) ! Multiply matrices for injection line together then compute inverse to pass to CREATE_PHASE_SPACE call prop_phase_space(mat_inv,twiss, twiss_start_inj_line) if(.not. use_lattice_twiss .and. twiss%betax /= 0 )then !use twiss parameters propagated backwards from inflector 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 field print '(a)',' BETA ZERO at START. INITIALIZE in LATTICE file or in INPUT.DAT' ! stop endif endif call create_phase_space(nmuons, create_new_distribution, new_file, muon_file, & epsdistr, epsx, epsy, tdistr, tlength, tsigma, pzdistr, pz, pzsigma, twiss,muons, twiss_ref, mat, mat_inv) !mat_inv propagates from end of inflector back to start of injection line ! call compute_beam_params(muons, 1, 'PropagatedBackToStart') do j = 1,nmuons if (muons(j)%coord%state /= alive$) cycle 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 quad_plate = quad_p runge_kutta_com%check_limits = .true. runge_kutta_com%check_wall_aperture = quad_plate !if true then scatter in plates print '(a,1x,l)', 'runge_kutta_com%check_wall_aperture = ', runge_kutta_com%check_wall_aperture runge_kutta_com%hit_when = wall_transition$ 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 i=0 if(nmuons == 1)then call write_single_particle_tbt(i,muons) 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)) endif branch%ele(1:branch%n_ele_track)%alias= 'trajectory' lat%branch(0)%ele(1:lat%branch(0)%n_ele_track)%alias = 'trajectory' DO i=1, nturns IF (i>1) runge_kutta_com%check_wall_aperture = .false. 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$ IF (muons(n)%coord%state /= alive$) CYCLE if(start_tracking_at_inflector_exit)then do j=1,lat%branch(0)%n_ele_track from_orbit(j)%vec = muons(n)%coord%vec from_orbit(j)%t = muons(n)%coord%t end do else from_orbit(0)%vec = muons(n)%coord%vec from_orbit(0)%t = muons(n)%coord%t endif ! call track_all(lat,from_orbit, ix_branch = 0, track_state = track_state, err_flag = err_flag) do j=1,lat%branch(0)%n_ele_track if(start_tracking_at_inflector_exit .and. j < lat%branch(0)%n_ele_track-1 ) cycle 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)) if(j < lat%branch(0)%n_ele_track)then if(lat%branch(0)%ele(j)%name == 'BACKLEG_START' )then from_orbit(j)%vec = from_orbit(j)%vec + start_orb%vec endif endif lat%branch(0)%ele(j)%tracking_method = ele_save%tracking_method muons_ele_inj(n,j)%coord = from_orbit(j) end do IF (nmuons==1 .and. .not. start_tracking_at_inflector_exit )then !if there is only one muon then compute its trajectory and twiss parameters along injection line call Bfield_along_axis do j=1,lat%branch(0)%n_ele_track call write_single_particle_by_element(i,lat%branch(0)%ele(j)%name,from_orbit(j)%vec, lat%branch(0)%ele(j)%s/circumference, & ((i-1)*lat%branch(0)%ele(lat%branch(0)%n_ele_track)%s +lat%branch(0)%ele(j)%s)/circumference, Q) call lat_make_mat6 (lat, j,from_orbit, ix_branch=0, err_flag=err_flag) end do call twiss_propagate_all(lat, ix_branch=0) do j=1,lat%branch(0)%n_ele_track if(make_movie)call step_floor_around_ring(lat,j,0,from_orbit,lat%branch(0)%ele(j-1)%floor%r) end do ! call optimize_twiss(lat, twiss, twiss_opt) !twiss_opt are the twiss parameters at the upstream end of the injection line ! that yield the desired twiss parameters in the inflector (namely twiss) ! superseded by prop_phase_space s_end=0.00001 orb_at_s%vec=from_orbit(0)%vec ix_end=1 call write_injection_line_trajectory(s_end, orb_at_s%vec, from_orbit(1)%vec,lat%ele(1)%name, ix_end,lat%ele(1)%floor%r) do while(s_end < lat%ele(lat%branch(0)%n_ele_track)%s +0.05) !the 0.05 is to make sure we overlap the very end if(s_end > lat%ele(lat%branch(0)%n_ele_track)%s) s_end = lat%ele(lat%branch(0)%n_ele_track)%s !track to the very end call twiss_and_track_at_s(lat, s_end, ele_at_s, from_orbit, orb_at_s, ix_branch=0,compute_floor_coords=.true.) ix_end = element_at_s (lat, s_end, .true.) s_rel = s_end - lat%ele(ix_end-1)%s call em_field_custom(lat%ele(ix_end), lat%param, s_rel,t_rel, orb_at_s, local_ref, field,calcd, err_flag) write(73,'(es12.4,1x,6es12.4,3es12.4)')s_end, orb_at_s%vec, field%B(1:3) call write_injection_line_trajectory(s_end, orb_at_s%vec, from_orbit(ix_end)%vec,lat%ele(ix_end)%name, ix_end, ele_at_s%floor%r) !end_orb%vec call write_injection_line_twiss(ele_at_s) s_end = s_end + 0.1 end do endif !nmuons = 1 to_orbit = from_orbit(ie_from) if(start_tracking_at_inflector_exit)to_orbit%vec = to_orbit%vec + start_orb%vec muons(n)%coord%vec = to_orbit%vec muons(n)%coord%t = to_orbit%t muons(n)%coord = to_orbit ! generate the same spin for every muon if (spin_tracking_on) muons(n)%coord%spin = (/(0.,0.),(1.,0.)/) 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) endif !for first turn only do n=1, nmuons if(allocated(muons_ele))muons_ele(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 call track_all(lat, co, nbranch, track_state, err_flag) if(nmuons == 1)then do j=1,branch%n_ele_track if(make_movie)call step_floor_around_ring(lat,j,nbranch,co ,branch%ele(j-1)%floor%r) end do endif if (i == 1 .and. n == 6) then do j=0, branch%n_ele_track write(107,'(es12.4)') co(j)%s !sampling end do end if 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)call compute_spin_single_muon(co(j), Q) call write_single_particle_by_element(i,branch%ele(j)%name,co(j)%vec, branch%ele(j)%s/circumference, & ((i-1)*branch%ele(branch%n_ele_track)%s +branch%ele(j)%s)/circumference,Q) 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)%vec, 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$ .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 if (co(branch%n_ele_track)%t >= co(0)%path_len .and. co(0)%t < co(0)%path_len) then do j = 1, branch%n_ele_track if (co(j)%t >= co(0)%path_len) 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.) ! print '(i10,1x,a16,1x,i10,1x,a9,es12.4,a15,es12.4)',j,branch%ele(j)%name, track%n_pt, 'co(j)%t =',co(j)%t, 'co(0)%path_len= ',co(0)%path_len ! do k=1,track%n_pt ! print '(a11,es12.4,a9,es12.4,a9,3es12.4)',' co(j-1)%t ',co(j-1)%t,' co(j)%t ',co(j)%t,'track_orb',track%orb(k)%s,track%orb(k)%t, co(0)%path_len ! end do if(allocated(muons_ele))muons_ele(n,j:branch%n_ele_track)%coord%state = lost$ call bracket_index(track%orb(:)%t, 0,track%n_pt,co(0)%path_len,ix) !bracket index where muon decayed ! print '(a,i10,a,i10)',' track%n_pt = ',track%n_pt,' ix = ',ix co_muon_end = track%orb(ix) if(ix < track%n_pt)then !interpolate dt=track%orb(ix+1)%t- track%orb(ix)%t dt1 = co(0)%path_len - track%orb(ix)%t dt2 = dt-dt1 co_muon_end%s = (dt2*track%orb(ix)%s+dt1*track%orb(ix+1)%s)/dt co_muon_end%t = (dt2*track%orb(ix)%t+dt1*track%orb(ix+1)%t)/dt co_muon_end%vec = (dt2*track%orb(ix)%vec+dt1*track%orb(ix+1)%vec)/dt endif call decayElectron(co_muon_end,co_electron(n),n) distance = float(i-1) + co_muon_end%s/circumference if(i==1 .and. co_muon_end%s == 0.00001) distance = 0. ! write(105,'(4es12.4,i10,3es12.4)') distance,co_muon_end%s, & ! circumference,circumference2,n,ie_from, from_orbit(ie_from)%s,co(0)%s,branch%ele(0)%s call decay_info(distance,n,co_electron(n),co_muon_end) !record decay information bmad_com%aperture_limit_on = .false. bmad_com%spin_tracking_on = .false. ! runge_kutta_com%check_limits = .false. do k=1,branch%n_ele_track !set aperture limits to zero. (Equivalent to setting to maximum) ele_limit(k,1) = branch%ele(k)%value(x1_limit$) !save aperture limits ele_limit(k,2) = branch%ele(k)%value(x2_limit$) ele_limit(k,3) = branch%ele(k)%value(y1_limit$) ele_limit(k,4) = branch%ele(k)%value(y2_limit$) branch%ele(k)%value(x1_limit$)=0. !set aperture limits to max branch%ele(k)%value(x2_limit$)=0. branch%ele(k)%value(y1_limit$)=0. branch%ele(k)%value(y2_limit$)=0. end do ! print '(a,i10,a,6es12.4)','n =',n,'co_electron(n) = ',co_electron(n)%vec ix_end = branch%n_ele_track call track_from_s_to_s(lat, co_muon_end%s,branch%ele(ix_end)%s,co_electron(n),co_elec_end,all_orb=co_elec,ix_branch=nbranch, track_state = track_state) !track electron to end of lattice ! do k=1,ix_end ! print '(a10,1x,i10,1x,a16,1x,11es12.4)',' electron ',k,branch%ele(k)%name,& ! co_elec_end%s,branch%ele(ix_end)%s,co_elec_end%vec,co_elec(k)%s,co_elec(k)%vec(1), co_elec(k)%p0c ! end do co_elec(ix_end) = co_elec_end do while(track_state == moving_forward$) !keep going around until the positron is lost co_elec(0)=co_elec(ix_end) call track_all(lat,co_elec, ix_branch=nbranch,track_state=track_state, err_flag = err_flag) !track electron around ring until it is lost ! do k=1,ix_end ! print '(a10,1x,i10,1x,a16,1x,7es12.4)',' electron ',k,branch%ele(k)%name,co_elec(k)%vec, co_elec(k)%p0c ! end do end do ! print '(a,a16,a,a,a15,1x,i10,1x,a)',' Positron lost at element ', branch%ele(track_state)%name, ' with coord state = ', & ! coord_state_name(co_elec(track_state)%state),' track_state = ',co_elec_end%state, & ! coord_state_name(co_elec_end%state) ! bmad_com%aperture_limit_on = .false. bmad_com%aperture_limit_on = .true. bmad_com%spin_tracking_on = .true. do k=1,branch%n_ele_track !restore aperture limits branch%ele(k)%value(x1_limit$)=ele_limit(k,1) branch%ele(k)%value(x2_limit$)=ele_limit(k,2) branch%ele(k)%value(y1_limit$)=ele_limit(k,3) branch%ele(k)%value(y2_limit$)=ele_limit(k,4) end do ! runge_kutta_com%check_limits = .true. 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 muons(n)%coord%state = lost$ muons(n)%lost = .true. 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 < co(0)%path_len)then muons(n)%coord%state = lost$ muons(n)%lost = .true. endif endif !muon_decay_on !********************************************** muons(n)%coord%vec = co(branch%n_ele_track)%vec muons(n)%coord%t = co(branch%n_ele_track)%t muons(n)%coord = co(branch%n_ele_track) 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) !write phase space after each element in 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) !collect data for 180 degree fiber monitor if (distance-floor(distance) == 0.5 .and. index(branch%ele(j)%name,'TRACK')/=0 ) call compute_moments(muons_temp,distance,moment_ele, nmuons,nmu, 2, branch%ele(j)%name, j) if (distance-floor(distance) == 0.5 .and. index(branch%ele(j)%name,'TRACK')/=0 ) call compute_moments(muons_temp,distance,moment_ele, nmuons,nmu, 2, branch%ele(j)%name, -1) !get distributions if (index(branch%ele(j)%name,'TRACKER2')/=0 ) call compute_beam_params(muons_temp,-i) ! calculate beam parameters if (index(branch%ele(j)%name,'TRACKER3')/=0 ) call compute_beam_params(muons_temp,-i-1000000) ! calculate beam parameters ! if (distance-floor(distance) == 0.5 .and. index(branch%ele(j)%name,'TRACK')/=0 ) call compute_beam_params(muons_temp,-i) ! calculate emittance if(i == 1)then if(write_phase_space_file)call write_phase_space(nmuons,muons_temp, branch%ele(j)%name, tot) !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 < 100)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 (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 < 100)call moment_range(moment) !compute range of moments only for last 100 turns if(i > 10)call collect_times(muons, nturns, NN, bins) !start to collect times of particles that are likely to survive if(nmuons == 1) call write_single_particle_tbt(i,muons) ENDDO ! nturns if(muon_decay_on)call electron_info_s(electrons_s,n_row,esteps*(1 + tot_track),n,.false.) !sample one electon's trajectory 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) print '(/,5a,/)',' Use plotting script ',"'",'sigma.gnu',"'",' to plot turn-by-turn average position and size' if (loop .and. vparam_id < 100)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) 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) if (vparam_id <100) vparam = vparam + delta_vparam if(.not. loop .or. vparam_id <= 0 .or. vparam_id > 100) exit end do lun = lunclose() print '(a,i10,a)',' Close ',lun,' units' if(record_fiber_data) then call fiber_records !process the fiber monitor data !call fiber_records !process the fiber monitor data call fft_analysis("fort.82",128,85,87,1,1) !analyze x data from 180 degree fiber monitor !call sleep(20) call fft_analysis("fort.95",128,96,97,2,2) !data from general simulation at 180 degree call fft_analysis("svd_results.dat",128,108,109,2,10) !call sleep(20) call fft_analysis("fort.83",128,86,88,1,3) !x data from 270 degree fiber monitor !call fft_analysis("fort.84",128,89,90,1,6) !y data from 180 degree fiber monitor !call fft_analysis("fort.79",128,77,78,1,7) !y data from 270 degree fiber monitor !call fft_analysis("params_in_ring_180.dat",128,98,99,3,4) !fit eta function for 180 degree fiber monitor !call fft_analysis("params_in_ring_180.dat",128,100,101,4,5) !fit beta function for 180 degree fiber monitor !call fft_analysis("params_in_ring_270.dat",128,102,103,3,8) !fit eta function for 270 degree fiber monitor !call fft_analysis("params_in_ring_270.dat",128,104,105,4,9) !fit beta function for 270 degree fiber monitor call svd_analysis(130) end if end