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(:) 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 (em_field_struct) field type (track_struct) track 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 lun integer lun32,lun16,lun26, lun19, lun29, lun33, lun49, lun71, lun72 integer lun_inj_branch, lun_ring_branch 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 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) 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) integer tot_max/0/ character*120 line, lat_file, lat_file_name/''/,new_file/' '/, muon_file/' '/ character*16 ele_name character*20 vparam_label(17)/'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' / integer vparam_column_index(17)/17,20,18,19,21,24,25,22,23,26,27,28,29,30,1,2,3/ logical err_flag, inf_end_us/.false./, inf_end_ds/.false./, quad_plate/.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 record_fiber_data/.true./ logical use_lattice_twiss/.false./ logical write_phase_space_file/.false./ logical save_element_by_element_info/.true./ 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 (quad_params_struct) quad_params type (field_file_struct) fringe_file, inflector_file 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 ! 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. 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) ! 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) !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 ! Scattering in inflector etc. is off while parser computes transfer matrices call initializeMaterials() call bmad_parser (lat_file, lat) lun=lunget() open(unit=lun, file="log.dat",access='append') write(lun,'(a16,1x,4a12,a12)')'Name','x1_limit','x2_limit','y1_limit','y2_limit','limit_on' do nbranch=0,1 if(nbranch==0)write(lun,'(a)')'Injection line' if(nbranch==1)write(lun,'(a)')'Ring' do j=1,lat%branch(nbranch)%n_ele_track write(lun,'(a16,1x,4es12.4,5x,l)')lat%branch(nbranch)%ele(j)%name, lat%branch(nbranch)%ele(j)%value(x1_limit$), lat%branch(nbranch)%ele(j)%value(x2_limit$), & lat%branch(nbranch)%ele(j)%value(y1_limit$), lat%branch(nbranch)%ele(j)%value(y2_limit$), bmad_com%aperture_limit_on end do end do close(unit=lun) ! 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, quad_params) ! 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 ! Calculate the ring's Twiss parameters, coupling, and dispersion parameters 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 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 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) 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 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,quad_params, kicker_params,tlength, pz, initial_offsets) bmad_com%aperture_limit_on = .true. 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' lun49=lunget() open(unit=lun49, file='track_backleg_inflector.dat') do i=1,lat%branch(0)%n_ele_track print '(7es12.4,1x,a)',lat%ele(i)%s,co(i)%vec(1:6), lat%ele(i)%name write(lun49,'(7es12.4,1x,a)')lat%ele(i)%s,co(i)%vec(1:6), lat%ele(i)%name end do close(unit=lun49) 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, 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, inf_aperture, inf_end_us, inf_end_ds,enerloss, muons, twiss_ref, 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 !print *, 'alive$',alive$,'downstream',downstream_end$ !print *, 'state,location,direction',muons(j)%coord%state,muons(j)%coord%location,muons(j)%coord%direction 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(unit,nmuons,muons,'', tot_inf_scatter) call init_moment(moment) call init_moment(moment_ele) lun16=lunget() open(unit=lun16, file = 'Beam_moments_tbt.dat') lun26 = lunget() open(unit=lun26, file = 'Beam_moments_by_element.dat') lun19 = lunget() open(unit=lun19, file= 'single_particle_tbt.dat') lun29 = lunget() open(unit=lun29, file = 'single_particle_by_element.dat') lun33 = lunget() open(unit=lun33,file='lost_muons.dat') lun71 = lunget() open(unit=lun71,file='injection_line_trajectory.dat') lun72 = lunget() open(unit=lun72,file='injection_line_twiss.dat') call compute_moments(muons,0.d0,moment, nmuons,nmu, lun16, 'BEGINNING',0) call compute_moments(muons,-lat%branch(0)%ele(lat%branch(0)%n_ele_track)%s/circumference,moment, nmuons,nmu, lun26,'BEGINNING',0) !call compute_moments(muons,-lat%branch(0)%ele(lat%branch(0)%n_ele_track)%s/circumference,moment, nmuons,nmu, 95,'BEGINNING',0) !used for FFT analysis end=lat%n_ele_track i=0 if(nmuons == 1)then write(lun19,'(a10,1x,7a12)')' turn ','x','xp','y','yp','z','pz','t' write(lun19,'(i10,1x,7es12.4)')i,muons(1)%coord%vec, muons(1)%coord%t write(lun29,'(a10,1x,9a12,2a10,1x,12x,a12,12x)')' turn ','x','xp','y','yp','z','pz','t','s','s total','state','muon','polarization' write(lun29,'(i10,1x,a16,1x,9es12.4)')i,lat%branch(0)%ele(0)%name,muons(1)%coord%vec, muons(1)%coord%t, & lat%ele(0)%s/circumference, & ((i-1)*lat%ele(lat%n_ele_track)%s +lat%ele(0)%s)/circumference write(lun71,'(a12,1x,6a12)') 's','x','xp','y','yp','z','pz' write(lun71,'(es12.4,1x,6es12.4)') lat%branch(0)%ele(0)%s,start_orb%vec(1:6) write(lun72,'(7a12,1x,a16)')'s', 'beta x','beta y','alpha x','alpha y', 'eta x','etap x', 'element' write(lun72,'(7es12.4,1x,a16)')lat%branch(0)%ele(0)%s, lat%branch(0)%ele(0)%a%beta,lat%branch(0)%ele(0)%b%beta,lat%branch(0)%ele(0)%a%alpha,lat%branch(0)%ele(0)%b%alpha, & lat%branch(0)%ele(0)%x%eta,lat%branch(0)%ele(0)%x%etap,lat%branch(0)%ele(0)%name 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 write(lun29,'(i10,1x,a16,1x,9es12.4,i10)')i,lat%branch(0)%ele(j)%name,from_orbit(j)%vec, co(j)%t, & 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, from_orbit(j)%state 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 write(lun71,'(es12.4,1x,6es12.4,1x,6es12.4,1x,a16,1x,i10,3es12.4)')s_end, orb_at_s%vec, from_orbit(1)%vec,lat%ele(1)%name, ix_end,lat%ele(1)%floor%r !end_orb%vec 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) write(lun71,'(es12.4,1x,6es12.4,1x,6es12.4,1x,a16,1x,i10,3es12.4)')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 write(lun72,'(7es12.4,1x,a16)')ele_at_s%s, ele_at_s%a%beta, & ele_at_s%b%beta, ele_at_s%a%alpha, ele_at_s%b%alpha, & ele_at_s%x%eta, ele_at_s%x%etap, ele_at_s%name s_end = s_end + 0.1 end do endif !nmuons = 1 to_orbit = from_orbit(ie_from) if(start_tracking_at_inflector_exit .and. nmuons == 1)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) ! call compute_beam_params(muons, 1, 'EndInfAfterScatter') if(write_phase_space_file)call write_phase_space(unit, 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) write(lun29,'(i10,1x,a16,1x,9es12.4,i10,i10,1x,3es12.4)')i,branch%ele(j)%name,co(j)%vec, co(j)%t, & branch%ele(j)%s/circumference, & ((i-1)*branch%ele(branch%n_ele_track)%s +branch%ele(j)%s)/circumference,co(j)%state,n, 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 if(first_loss)then write(lun33,'(a10,1x,9a12)')'muon','s','x','xp','y','yp','z','pz','element','turn' print *,' lun33 = ',lun33,' lun26= ',lun26 endif write(lun33,'(i10,1x,7es12.4,1x,a12,1x,i10)')n,branch%ele(j)%s, co(j)%vec, branch%ele(j)%name,i first_loss = .false. 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 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 = (dt1*track%orb(ix)%s+dt2*track%orb(ix+1)%s)/dt co_muon_end%t = (dt1*track%orb(ix)%t+dt2*track%orb(ix+1)%t)/dt co_muon_end%vec = (dt1*track%orb(ix)%vec+dt2*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(41,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 = .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 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 forall(n=1:nmuons)muons_temp(n) = muons_ele_inj(n,j) call compute_moments(muons_temp,distance,moment_ele, nmuons,nmu, lun26, lat%branch(0)%ele(j)%name, j) if(write_phase_space_file)call write_phase_space(lun_inj_branch,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,38, 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, lun26, 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, 95, branch%ele(j)%name, j) 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(lun_ring_branch,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,38, 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, lun16, 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,39,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) write(lun19,'(i10,1x,7es12.4)') i, muons(1)%coord%vec, muons(1)%coord%t ENDDO ! nturns if(muon_decay_on)call electron_info_s(electrons_s,n_row,esteps*(1 + tot_track),42,n,.false.,43) !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(unit,nmuons,muons,'EndOfTracking', tot) print '(/,5a,/)',' Use plotting script ',"'",'sigma.gnu',"'",' to plot turn-by-turn average position and size' if (first .and.loop)then lun32 = lunget() open(unit = lun, file = 'VparamDependence.dat') write(lun32,'(a,i4)')'vparam_id =', vparam_id write(lun32,'(a,a1,a,a1)')'set xlabel ','"',vparam_label(vparam_id),'"' write(lun32,'(a,i5)')'column1 =',vparam_column_index(vparam_id) write(lun32,'(a,i10)')'nmuons =', nmuons write(lun32,'(a)')'exit' write(lun32,'(/,a,/)')'#' write(lun32,'(16a16,9a12,a12,a12, a12, a12,a12)') & 'K1[Guass]-1', 'K2[Gauss]-2', 'K3[Gauss]-3', & 'depth-4',' depth-5','depth-6','depth-7','depth-8','depth-9', & 'depth_SIG(1,1)-10','depth_SIG(3,3)-11','depth_SIG(5,5)-12','depth_SIG(6,6)-13', & 'n muons-14',' after_inf-15','remaining-16','beta x-17','eta x-18','kicker_field-19','beta y-20',' kickWid-21','dtRise-22','dtFall-23','x init','px init', & 'inf_field','pz init','etap','alphax','Kick_tStart' first=.false. endif !open lun32 if(loop)then write(lun32,'(3es16.3,10es16.6,3i16,es12.4,es12.4,5es12.4, 2es12.4, es12.4,es12.4, es12.4, es12.4,es12.4)') & kicker_params%kicker_field(1) * 1.e4, & ! Gauss kicker_params%kicker_field(2) * 1.e4, & ! Gauss kicker_params%kicker_field(3) * 1.e4, & ! Gauss (moment%max_ave-moment%min_ave) *1.e3, & ! mm or mrad (sqrt(moment%max_sigma(1,1))-sqrt(moment%min_sigma(1,1)))*1.e3, & ! mm or mrad (sqrt(moment%max_sigma(3,3))-sqrt(moment%min_sigma(3,3)))*1.e3, & ! mm or mrad (sqrt(moment%max_sigma(5,5))-sqrt(moment%min_sigma(5,5)))*1.e3, & ! mm or mrad (sqrt(moment%max_sigma(6,6))-sqrt(moment%min_sigma(6,6)))*1.e3, & ! mm or mrad nmuons, & tot_inf_scatter, & tot, twiss%betax, twiss%etax, maxval(kicker_params%kicker_field(1:3))*10000, twiss%betay, kicker_params%kick_width(1)*1.e9,& kicker_params%dtRise(1)*1.e9, kicker_params%dtFall(1)*1.e9, & initial_offsets%x_mean*1.e3, initial_offsets%pxmean*1.e3, inflector_field, initial_offsets%pzmean*1.e2, twiss%etapx, & twiss%alphax, kicker_tStart(1)*1.e9 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) close(lun16) close(lun19) close(lun29) close(lun26) close(lun33) close(lun71) close(lun72) ! twiss%etax = twiss%etax + 0.2 ! twiss%betax = twiss%betax + 0.2 ! kicker_params%kicker_field(2) = kicker_params%kicker_field(2) + 7.5e-4 vparam = vparam + delta_vparam if(.not. loop .or. vparam_id <= 0) exit end do if(loop)close(unit=lun32) 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) !call sleep(20) call fft_analysis("fort.95",128,96,97,2,2) !call sleep(20) call fft_analysis("fort.83",128,86,88,1,3) end if end