! Program to ! - create phase space distribution at point of injection (from linac) into the synchrotron. Average energy is 150 MeV ! - Track 'beam' from injection point to extraction point (thin septum). ! - Continue to track each turn to extraction point while ramping synchrotron magnets from 150 MeV to flat top (branch 0) ! - Ramp 3 component bump at turn and with rate specified in ! - Ramp quads that shift tunes onto resonance as specified in ! - Ramp sextupole that shift tunes onto resonance as specified in ! Particles that clear the first (thin) septum continue in the extraction line branch (branch 1) ! ! Particles that clear each of the septa in the extraction branch receive the appropiate kick ! ! (Note: There are two branches. branch 0 - synchrotron. branch 1 - elements immediately after thin septum all the way to target. ! this will include several synchrotron magnets. But the trajectory is in the fringe field of these magnets. All of the ! information about the fringe can be in branch 1.) ! ! Input parameters are in ! If the number of particles to track is unity, use rather than program ramp use bmad use beam_mod use multibunch_mod use multibunch_interface use bmadz_other_code_interface use ramp_parameters implicit none type (lat_struct) lat type (branch_struct), pointer :: branch ! ring type (coord_struct), allocatable::co(:), co_fft(:), coa(:), co_many(:), co_2000(:) type (coord_struct) co_start type (ele_struct), pointer :: fork_ele type (beam_init_struct) beam_init type (beam_struct) beam, ext_beam, beam_on_target type (moment_struct), allocatable :: bunch_moments(:) type (ele_pointer_struct), allocatable :: eles(:) type (bunch_struct) bunch_dum type (all_pointer_struct) a_ptr integer ix, nargs integer i, j, k integer lun integer lun_trajectory, lun_beam_trajectory integer lun_ext_phasespace integer lun_phasespace integer lun_target, lun_sep1_beam, lun_sep2_beam, lun_sep0_beam, lun_sep1_phasespace, lun_sep2_phasespace, lun_sep0_phasespace, lun_target_moments integer lunx integer lun_tune integer nturns integer ios integer n_fft/256/ integer nparticles integer n_clear, n_lost, n_hit integer ix_sext integer n_extract_tot(0:2)/3*0/, n_lost_tot(0:2)/3*0/, n_sep_tot(0:2)/3*0/ integer Nbin(-200:200), nn integer n_loc integer track_state integer ix_inj, ix_essep, ix_l3, ie_from, ix_magsept1(0:1), ix_magsept2(0:1)/0,0/, ix_extraction_fork integer, parameter:: end_of_line$=-1 integer status/0/, system integer datetime_values(8) integer ix_custom1, ix_custom2 real(rp) start, finish real(rp) delta real(rp) e_init, omega, rf_phase_twopi real(rp) initial_offset(6) real(rp) Q(3), phase(3) real(rp) qvert,qhorz, ext_time real(rp) emit_x, emit_y,sig_z, sig_e real(rp) x real(rp) width/0.0002/ real(rp) dx/0.002/ real(rp) bump_amp/1./ real(rp) hkick(2000) character*120 string character*140 lat_file/'*'/ character*120 line, lineout character*10 iteration character*16 date, time, zone logical err logical ramp_synch, phase_space/.false./ logical sext/.false./ logical aperture_limit_on/.false./ logical scale_emit/.false./ logical err_flag namelist/input/lat_file, & ! lattice e_init, & !starting energy ramp_synch, & ! ramp t0, & ! start time nturns, & ! number of turns rf_phase_twopi, & ! RF phase initial_offset, & ! coord vec n_fft, & ! number of turns in each fft nparticles, & ! number of particles qvert,qhorz, & ! strength of detuning vertical and horizontal quads sextk2, & ! strength of sextupole octk3, & ! strength of octupole ext_time, & ! extraction start time, when pulse bump is on full turn_on_time, & ! time for detuning quads to get to full strength as fractin of ext_time emit_x, emit_y, sig_z, sig_e, & ! initial beam parameters ext_param, & ! multiplier of extraction sextupole or octupole x_septum, t_septum, k_septum, & !extraction septum - center line to far edge of septum (0:2), septum thickness (0:2), septum kick (0:2) x_wall, & !center line to wall(0:2) of vacuum chamber, thickness t1, & ! time to stop increase in octupole oct_ramp, quad_ramp, bump_ramp,sext_ramp, & ! oct_ramp(0:7000), octupole function, quad_ramp(0:7000) quad function, bump_ramp(0:7000) bump ramp * 60 hz phase_space, & ! if true write phase space for each particle each turn hbump, bump_bend, & ! hbump(1:3) in radians and the names of the 3 bump bends aperture_limit_on, & ! set true to turn on aperture limits scale_emit, & !scale emittance computed in init_beam_distribution to injection energy bump_amp ! amplitude of pulse bump call cpu_time(start) bmad_com%auto_bookkeeper = .false. !initialize ramp functions oct_ramp = -1. quad_ramp = -1. bump_ramp = -1. sext_ramp = -1. !set default septum parameters x_septum(0) = -0.031 !displacement of thin septum 0 from center of chamber x_septum(1) = -0.03 !displacement of septum 1 from center of chamber x_septum(2) = 0.03 !displacement of septum 2 from center of chamber x_wall(0) = -0.045 !distance from center of chamber to outer wall of thin septum (septum 0) x_wall(1) = -0.045 !distance from center of chamber to outer wall of septum 1 x_wall(2) = 1. !distance from center of chamber to outer wall of septum 1 t_septum(0) = 0.0001 !thickness of septum 1 t_septum(1)=0.006 !thickness of septum 1 t_septum(2)=0.006 !thickness of septum 1 call date_and_time(date,time,zone,datetime_values) print '(a10,a10,a1,a10)',' date = ',date,' time = ', time ! create subdirectory for output of this run directory = trim(date)//'_'//trim(time(1:6)) status= system('mkdir ' // directory) if(status == 0)print '(a20,a)',' create directory = ', trim(directory) ! read input file lun = lunget() print *,'lunget() = lun ',lun open(unit=lun, file='input.in', STATUS ='old') read(lun, nml=input, IOSTAT=ios) rewind(lun) read(lun, nml=input, IOSTAT=ios) close(unit=lun) write(6,nml=input) OPEN (UNIT=5, FILE='input.in', STATUS='old', IOSTAT=ios) lun=lunget() open(unit=lun, file=trim(directory)//'/input.in') do while(ios == 0) read(5,'(a)', IOSTAT=ios)string if(ios == 0)write(lun,'(a)')string end do close(unit=lun) qv = qvert qh=qhorz t_ext = ext_time ! set amplitude of 3 component bump hbump(1:3) = bump_amp*hbump(1:3) !ready an input file from the command line to read lattice file if(lat_file == '*')then nargs = cesr_iargc() if(nargs == 1)then call cesr_getarg(1,lat_file) print *, 'Using ', trim(lat_file) else lat_file = 'bmad.' print '(a,$)',' Lattice file name ? (default= synch.lat_ext_branch_l3_inside) ' read(5,'(a)') line call string_trim(line, lineout, ix) lat_file = lineout ! if(ix == 0) lat_file = '/home/jdp279/Documents/Track6610/Outgoing/darkphoton/files/synch.lat' if(ix == 0) lat_file = 'synch.lat_ext_branch_l3_inside' endif 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) call bmad_parser(lat_file, lat) print '(a)',' Bmad_parser successful' call reallocate_coord(co,lat%n_ele_track) call reallocate_coord(co_many,lat%n_ele_track) call reallocate_coord(co_fft, nturns) call reallocate_coord(coa,n_fft) call lattice_bookkeeper(lat) ! find location of Injection point (POSINJSEP), electrostatic septum (ESSEP), and BPM at l3 (L3) call lat_ele_locator('POSINJSEP',lat,eles,n_loc, err) ! Element POSINJSEP is the injection point into the synchrotron ix_inj = eles(1)%ele%ix_ele print '(a14,a16,a9,i10)',' Injection at ',eles(1)%ele%name,' element ',ix_inj call lat_ele_locator('ESSEPT',lat,eles,n_loc,err) !Elemet ESSEPT is the thin extraction septum (septup 0) ix_essep = eles(1)%ele%ix_ele hkick(ix_essep) = eles(1)%ele%value(kick$) print '(a22,1x,a16,1x,a9,1x,i10,1x,a9,1x,es12.4)',' First thin septum at ',eles(1)%ele%name,' element ',ix_essep,' hkick = ', hkick(ix_essep) lat%ele(ix_essep)%value(kick$)=0. call lat_ele_locator('L3',lat,eles,n_loc,err) ix_l3 = eles(1)%ele%ix_ele print '(a14,a16,a9,i10)',' Quad at ',eles(1)%ele%name,' element ',ix_l3 do j=0,1 !both branches call lat_ele_locator('MAGSEPT1',lat,eles,n_loc,err,ix_dflt_branch=j) !Element MAGSEPT1 is (septum 1) ix_magsept1(j) = eles(1)%ele%ix_ele hkick(ix_magsept1(j)) = eles(1)%ele%value(kick$) lat%branch(j)%ele(ix_magsept1(j))%value(kick$)=0. print '(a15,1x,a16,1x,a9,1x,i10,1x,a9,es12.4)',' Second septum ',eles(1)%ele%name,' element ',ix_magsept1(j),' hkick = ', hkick(ix_magsept1(j)) call lat_ele_locator('MAGSEPT2',lat,eles,n_loc,err, ix_dflt_branch=j) !Element MAGSEPT2 is (septum 2) ix_magsept2(j) = eles(1)%ele%ix_ele hkick(ix_magsept2(j)) = eles(1)%ele%value(kick$) lat%branch(j)%ele(ix_magsept2(j))%value(kick$)=0. if(eles(1)%ele%key == sbend$)then hkick(ix_magsept2(j)) = eles(1)%ele%value(g$) lat%branch(j)%ele(ix_magsept2(j))%value(g$)=0. endif print '(a14,1x,a16,1x,a9,1x,i10,1x,a9,es12.4,a9,es12.4)',' Third septum ',eles(1)%ele%name,' element ',ix_magsept2(j),' hkick = ', hkick(ix_magsept2(j)),' g = ',hkick(ix_magsept2(j)) end do call lat_ele_locator('RES_EXT_SEXT',lat,eles,n_loc,err) !Element RES_EXT_SEXT is the sextupole that drives the resonance ix_sext = eles(1)%ele%ix_ele print '(a23,a16,a9,i10)',' Resonstance extraction sext',eles(1)%ele%name,' element ',ix_sext call lat_ele_locator('EXTRACTION_FORK',lat,eles,n_loc,err) ! Fork from ring (branch 0) to extraction line (branch 1) ix_extraction_fork = eles(1)%ele%ix_ele print '(a23,a16,a9,i10)',' To extracted beam line ',eles(1)%ele%name,' element ',ix_extraction_fork lat%param%ixx = 0 call closed_orbit_calc(lat,co,6) print '(a)',' Closed_orbit_calc successful' call track_all(lat,co) print '(a)','Lat_make_mat6 successful' call lat_make_mat6(lat,-1) ! compute transfer matrices for all of the elements in the lattice call twiss_at_start(lat) ! compute twiss parameters at the starting point print '(a)','Twiss_at_start successful' call twiss_propagate_all(lat) ! propagate twiss parameters to all elements in the lattice print '(a)','Twiss_propagate_all successful' call write_septa_info(lat,ix_inj, ix_essep, ix_l3, ie_from, ix_magsept1(0), ix_magsept2(0), ix_extraction_fork, ix_sext) ie_from = lat%branch(1)%ix_from_ele fork_ele => lat%branch(0)%ele(ie_from) call transfer_twiss(fork_ele, lat%branch(1)%ele(0)) call twiss_propagate_all(lat,ix_branch=1, err_flag =err) call writeinfo(lat) ! write initial twiss parameters to file call floor_plan(lat) ! write physical coordinates of relevant components ! initialize septum kicks lat%ele(ix_essep)%value(kick$)=hkick(ix_essep) lat%branch(0)%ele(ix_magsept1(0))%value(kick$)=hkick(ix_magsept1(0)) lat%branch(0)%ele(ix_magsept2(0))%value(kick$)=hkick(ix_magsept2(0)) lat%branch(1)%ele(ix_magsept1(1))%value(kick$)=hkick(ix_magsept1(1)) lat%branch(1)%ele(ix_magsept2(1))%value(kick$)=hkick(ix_magsept2(1)) print '(a,5es12.4)',' septa kicks = ',hkick(ix_essep),hkick(ix_magsept1(:)),hkick(ix_magsept2(:)) ! call track_sept1_end(lat,co, ix_essep, ix_magsept1, ix_magsept2, ix_extraction_fork) lun_trajectory = lunget() ! print *,'lunget() = lun2 ',lun2 open(unit=lun_trajectory, file = trim(directory)//'/trajectory.dat') !Turn by turn trajectory in single particle mode lun_beam_trajectory = lunget() ! print *,'lunget() = lun_beam_trajectory ',lun_beam_trajectory open(unit=lun_beam_trajectory, file = trim(directory)//'/beam_trajectory.dat') !Turn by turn centroid and width in multiparticle mode lun_phasespace = lunget() ! print *,'lunget() = lun_phasespace ',lun_phasespace if(phase_space)open(unit=lun_phasespace, file=trim(directory)//'/phasespace.dat') ! Turn by turn phase space of circulating beam lun_target = lunget() ! print *,'lunget() = lun_target ',lun_target open(unit=lun_target, file = trim(directory)//'/target_beam.dat') ! Particles at the target lun_sep1_beam = lunget() ! print *,'lunget() = lun_sep1_beam ',lun_sep1_beam open(unit=lun_sep1_beam, file = trim(directory)//'/sep1_beam.dat') ! Particles that clear sep1 (MAGSEPT1) lun_sep2_beam = lunget() ! print *,'lunget() = lun_sep2_beam ',lun_sep2_beam open(unit=lun_sep2_beam, file = trim(directory)//'/sep2_beam.dat') !Particles that clear sep2 (MAGSEPT2) lun_sep0_beam = lunget() ! print *,'lunget() = lun_sep0_beam ',lun_sep0_beam open(unit=lun_sep0_beam, file = trim(directory)//'/sep0_beam.dat') !Particles that clear sep0 (ESSEPT) lun_sep0_phasespace=lunget() ! print *,'lunget() = lun_sep0_phasespace ',lun_sep0_phasespace open(unit=lun_sep0_phasespace, file = trim(directory)//'/lun_sep0_phasespace.dat') !Phase space of particles that clear sep0 lun_sep1_phasespace=lunget() ! print *,'lunget() = lun_sep1_phasespace ',lun_sep1_phasespace open(unit=lun_sep1_phasespace, file = trim(directory)//'/lun_sep1_phasespace.dat') !Phase space of particles that clear sep1 lun_sep2_phasespace=lunget() ! print *,'lunget() = lun_sep2_phasespace ',lun_sep2_phasespace open(unit=lun_sep2_phasespace, file = trim(directory)//'/lun_sep2_phasespace.dat') !Phase space of particles that clear sep2 lun_ext_phasespace=lunget() open(unit=lun_ext_phasespace, file = trim(directory)//'/lun_ext_phasespace.dat') lun_target_moments = lunget() ! print *,'lunget() = lun_target_moments ',lun_target_moments open(unit=lun_target_moments, file=trim(directory)//'/target_moments.dat') !Beam moments at the target lun_tune = lunget() open(unit=lun_tune, file=trim(directory)//'/tune_matrix_turn.dat') write(lun_tune,'(a10,2a12)')'turn', 'x tune','y tune' do j=0,1 do i=1,lat%branch(j)%n_ele_max ix_custom1 = scratch1$ ix_custom2 = scratch2$ if(lat%branch(j)%ele(i)%key == rbend$ .or. lat%branch(j)%ele(i)%key == quadrupole$ .or.lat%branch(j)%ele(i)%key == sbend$)then ! save reference values for k1, k2 lat%branch(j)%ele(i)%value(ix_custom1) = lat%branch(j)%ele(i)%value(k1$) lat%branch(j)%ele(i)%value(ix_custom2) = lat%branch(j)%ele(i)%value(k2$) endif if(lat%branch(j)%ele(i)%key == multipole$)then !save reference values of multipoles allocate(lat%branch(j)%ele(i)%r(1:2,0:n_pole_maxx,0:1)) lat%branch(j)%ele(i)%r(1,0:n_pole_maxx,1) = lat%branch(j)%ele(i)%a_pole(:) lat%branch(j)%ele(i)%r(2,0:n_pole_maxx,1) = lat%branch(j)%ele(i)%b_pole(:) endif if(index(lat%branch(j)%ele(i)%name,'RES_EXT_SEXT')/=0)then lat%branch(j)%ele(i)%value(ix_custom2) = lat%branch(j)%ele(i)%value(k2$) !save reference value of sextupole print '(a11,a16,a7,a7,es12.4,a8,es12.4)',' Sextupole ',lat%branch(j)%ele(i)%name,' found ',' k2 = ',lat%branch(j)%ele(i)%value(k2$),' at s = ',lat%branch(j)%ele(i)%s sext = .true. endif if(index(lat%branch(j)%ele(i)%name,'ESSEPT')/=0)then if(k_septum(0) /= 0.)lat%branch(j)%ele(i)%value(kick$) = k_septum(0) lat%branch(j)%ele(i)%value(ix_custom2) = lat%branch(j)%ele(i)%value(kick$) print '(a11,a16,a7,1x,a,i10,1x,a,es12.4)',' Es sept ',lat%branch(j)%ele(i)%name,' found ',' ix_custom2 =', & ix_custom2, ' kick =',lat%branch(j)%ele(i)%value(ix_custom2) endif if(index(lat%branch(j)%ele(i)%name,'MAGSEPT1')/=0)then if(k_septum(1) /= 0.)lat%branch(j)%ele(i)%value(kick$) = k_septum(1) lat%branch(j)%ele(i)%value(ix_custom2) = lat%branch(j)%ele(i)%value(kick$) print '(a11,a16,a7,1x,a,i10,1x,a,es12.4)',' Mag sept1 ',lat%branch(j)%ele(i)%name,' found ',' ix_custom2 =' & ,ix_custom2, ' kick =',lat%branch(j)%ele(i)%value(ix_custom2) endif if(index(lat%branch(j)%ele(i)%name,'MAGSEPT2')/=0)then if(k_septum(2) /= 0.)lat%branch(j)%ele(i)%value(kick$) = k_septum(2) lat%branch(j)%ele(i)%value(ix_custom2) = lat%branch(j)%ele(i)%value(kick$) if(lat%branch(j)%ele(i)%key == sbend$)lat%branch(j)%ele(i)%value(ix_custom2) = lat%branch(j)%ele(i)%value(g$) print '(a11,a16,a7,1x,a,i10,1x,a,es12.4)',' Mag sept2 ',lat%branch(j)%ele(i)%name,' found ' ,' ix_custom2 =' & ,ix_custom2, ' kick =',lat%branch(j)%ele(i)%value(ix_custom2) endif if(lat%branch(j)%ele(i)%key == rfcavity$)then lat%branch(j)%ele(i)%value(phi0$) = rf_phase_twopi/twopi ! set RF phase print '(a,i10,a,es12.4)',' Element ',i, lat%branch(j)%ele(i)%name, lat%branch(j)%ele(i)%value(phi0$) endif end do end do !both branches if(.not. sext)then print *,' No resonance driving sextupole found ' stop endif omega = twopi * 60 co_start = co(ix_inj) if(ramp_synch)then lat%param%ixx = 1 !if 0 no ramp. delta = e_init/lat%ele(0)%value(E_TOT$) !fractional energy at injection energy. is the energy out of the linac t0 = acos(1.-2*delta)/omega ! fractional energy delta = (1- cos(omega*t))/2 co(ix_inj)%vec(6) = -(1-delta) co_start = co(ix_inj) call init_coord(co_start, co(ix_inj), lat%ele(ix_inj), downstream_end$, positron$) co_start%t = t0 print '(a,es12.4)',' Ramp initiated, co(0)%vec(6) = ', co(0)%vec(6) endif ! co_start%vec = co_start%vec + initial_offset print '(a20,i10,a20,6es12.4)',' nparticles = ',nparticles,' at line 346', co_start%vec if(nparticles == 1)then k=0 do while(co_start%vec(1) < 0.0012) co(ix_inj) = co_start co_many(ix_inj) = co_start print '(a,6es12.4)','co_start%vec',co_start%vec print '(a,6es12.4)','co(0)%vec',co(0)%vec write(lun_trajectory,'(/,a,a7,es12.4,/)')'#','xoffset = ', co(ix_inj)%vec(1) call track_many(lat,co_many,ix_inj,ix_essep,direction=1,ix_branch=0, track_state=track_state) !track from injection point to thin septum do i=1,nturns print *,'turn = ',i if(i==2000)call turn_2000(lat, co_many(0)) !record some details call track_many(lat,co_many,ix_essep,ix_essep,direction=1,ix_branch=0, track_state=track_state) !track full turn at thin septum co_fft(i) = co_many(ix_l3) write(lun_trajectory, '(i10,1x,es12.4,1x,6es12.4,1x,es12.4)')i,i*lat%ele(lat%n_ele_track)%s, co_many(ix_essep)%vec,co_many(ix_essep)%t if(track_state /= moving_forward$) write(lun_trajectory,'(a23,1x,i5,1x,a10,a16,a10,1x,6es12.4)')'Particle lost at index ',track_state,' element =',lat%ele(track_state)%name, & 'vec = ',co_many(track_state)%vec if(track_state /= moving_forward$)exit if((i/10)*10 == i )then call lat_make_mat6(lat,ix_ele = -1,ref_orb = co_many, ix_branch=0) ! compute transfer matrices for all of the elements in the lattice call twiss_at_start(lat) ! compute twiss parameters at the starting point, and tunes call write_full_turn(lat,i,co_many, lat%n_ele_track) write(lun_tune,'(i10,1x,2es12.4)')i,lat%a%tune/twopi, lat%b%tune/twopi !tunes from full turn matrix endif end do !nturns co_start%vec(1) = co_start%vec(1)+dx end do !loop over offsets (dx) at injection endif ! special case for single particle ! these 3 arrays count particles that clear septa, are lost at septa, or hit the septum wall n_clear = 0 n_lost = 0 n_hit = 0 if(nparticles > 1)then beam_init%n_bunch = 1 beam_init%dt_bunch = 0 beam_init%n_particle=nparticles beam_init%bunch_charge = 1 beam_init%a_emit = emit_x beam_init%b_emit = emit_y beam_init%sig_z = sig_z beam_init%sig_e = sig_e*delta beam_init%center = co_start%vec allocate(bunch_moments(0:beam_init%n_bunch)) print '(a,a16)','Injection at element ',lat%ele(ix_inj)%name print '(a10,es12.4,a11,es12.4, a10,es12.4)', 'beta x = ',lat%ele(ix_inj)%a%beta, ' alpha x = ',lat%ele(ix_inj)%a%alpha, 'eta x = ',lat%ele(ix_inj)%x%eta print '(a10,es12.4,a11,es12.4, a10,es12.4)', 'beta y = ',lat%ele(ix_inj)%b%beta, ' alpha y = ',lat%ele(ix_inj)%b%alpha, 'eta y = ',lat%ele(ix_inj)%y%eta call init_beam_distribution(lat%ele(ix_inj),lat%param, beam_init, beam) call init_beam_distribution(lat%branch(1)%ele(0),lat%param, beam_init, ext_beam) call init_beam_distribution(lat%branch(1)%ele(lat%branch(1)%n_ele_track),lat%param, beam_init, beam_on_target) beam_on_target%bunch(1)%particle(:)%state = not_set$ forall(i=1:size(ext_beam%bunch(1)%particle(:)))ext_beam%bunch(1)%particle(i)%state = lost$ ! Adjust for energy. That is, energy at start of ramp is not the energy for which distribution is created ! if(scale_emit)forall(i=1:6)beam%bunch(1)%particle(:)%vec(i) = beam%bunch(1)%particle(:)%vec(i)/sqrt(delta) ! forall(i=1:6)beam%bunch(1)%particle(:)%vec(i) = co_start%vec(i)+beam%bunch(1)%particle(:)%vec(i) + initial_offset(i) beam%bunch(1)%particle(:)%vec(2) = beam%bunch(1)%particle(:)%vec(2)*(1+beam%bunch(1)%particle(:)%vec(6)) beam%bunch(1)%particle(:)%vec(4) = beam%bunch(1)%particle(:)%vec(4)*(1+beam%bunch(1)%particle(:)%vec(6)) beam%bunch(1)%particle(:)%t = co_start%t write(lun_beam_trajectory,'(a10,2a72,1x,a10,1x,2a12,1x,3a10,1x,3a10,1x,3a10)')'turn','bunch_moments%ave(1:6)', 'bunch_moments%sigma(1:6,1:6)','n_alive','ave time','sigma_t', & 'n_ext(0)','n_lost(0)','n_sep(0)','n_ext(1)','n_lost(1)','n_sep(1)','n_ext(2)','n_lost(2)','n_sep(2)' if(phase_space)write(lun_phasespace, '(6a12)') 'x','xp','y','yp','z','zp' bmad_com%aperture_limit_on = aperture_limit_on i=0 if(phase_space)then write(lun_phasespace, '(/,a1,a7,i10,/)')'#', 'turn = ',i do j=1,size(beam%bunch(1)%particle(:)) if(beam%bunch(1)%particle(j)%state == alive$)write(lun_phasespace, '(6es12.4)')beam%bunch(1)%particle(j)%vec end do endif ! track from injection point to septum if (ix_inj < ix_essep) call track_beam(lat,beam,ele1=lat%ele(ix_inj),ele2=lat%ele(ix_essep), err=err) if(ix_inj > ix_essep)then call track_beam(lat,beam,ele1=lat%ele(ix_inj),ele2=lat%ele(lat%n_ele_track), err=err) call track_beam(lat,beam,ele1=lat%ele(0),ele2=lat%ele(ix_essep), err=err) endif do i=1,nturns ! if(abs(lat%branch(0)%ele(ix_sext)%value(k2$)) > 3) then ! print *,' sextupole is on at ', lat%branch(0)%ele(ix_sext)%value(k2$) ! endif call track_beam(lat,beam,ele1=lat%ele(ix_essep),ele2=lat%ele(lat%n_ele_track), err=err) !track from septum to end call track_beam(lat,beam,ele1=lat%ele(0),ele2=lat%ele(ix_essep), err=err) ! track from beginning to septum if(phase_space)then write(lun_phasespace, '(/,a1,a7,i10,/)')'#', 'turn = ',i do j=1,size(beam%bunch(1)%particle(:)) if(beam%bunch(1)%particle(j)%state == alive$)write(lun_phasespace, '(6es12.4)')beam%bunch(1)%particle(j)%vec end do endif if(i > ext_time/2.52131e-6)then ! begin extraction call check_septum_clearance(beam%bunch(1), n_clear,n_hit, n_lost, ext_beam%bunch(1), 0 ) !thin septum n_extract_tot(0) = n_extract_tot(0) + n_clear n_lost_tot(0) = n_lost_tot(0) + n_lost n_sep_tot(0) = n_sep_tot(0) + n_hit ! Track particles that clear septum 0 to septum 1 if(n_clear > 0)then call write_phase_space(lun_sep0_beam,i,ext_beam) !particles that clear septum 0 do j=0,ix_magsept1(1)-1 !track from start of extraction branch to magsept1 call track_beam(lat,ext_beam,ele1=lat%branch(1)%ele(j), ele2=lat%branch(1)%ele(j+1), err=err) !track from septum 1 to 2 in ext branch do k=1,size(ext_beam%bunch(1)%particle(:)) if(ext_beam%bunch(1)%particle(k)%state == alive$) & write(lun_sep0_phasespace,'(i10,1x,a12,es12.4,1x,i10,1x,6es12.4,1x,i10)')j+1,lat%branch(1)%ele(j+1)%name, & lat%branch(1)%ele(j+1)%s,k, ext_beam%bunch(1)%particle(k)%vec, ext_beam%bunch(1)%particle(k)%state end do end do call check_septum_clearance(ext_beam%bunch(1), n_clear, n_hit, n_lost, bunch_dum, 1) n_extract_tot(1) = n_extract_tot(1) + n_clear n_lost_tot(1) = n_lost_tot(1) + n_lost n_sep_tot(1) = n_sep_tot(1) + n_hit endif !septum 0-1 ! Track particles that clear septum 1 to septum 2 if(n_clear > 0)then call write_phase_space(lun_sep1_beam,i,ext_beam) !particles that clear septum 1 do j=ix_magsept1(1),ix_magsept2(1)-1 call track_beam(lat,ext_beam,ele1=lat%branch(1)%ele(j), ele2=lat%branch(1)%ele(j+1), err=err) !track to from septum 1 to 2 in ext branch do k=1,size(ext_beam%bunch(1)%particle(:)) if(ext_beam%bunch(1)%particle(k)%state == alive$) & write(lun_sep1_phasespace,'(i10,1x,a12,es12.4,1x,i10,1x,6es12.4,1x,i10)')j+1,lat%branch(1)%ele(j+1)%name,& lat%branch(1)%ele(j+1)%s,k, ext_beam%bunch(1)%particle(k)%vec, ext_beam%bunch(1)%particle(k)%state end do end do call check_septum_clearance(ext_beam%bunch(1), n_clear, n_hit, n_lost, bunch_dum, 2) n_extract_tot(2) = n_extract_tot(2) + n_clear n_lost_tot(2) = n_lost_tot(2) + n_lost n_sep_tot(2) = n_sep_tot(2) + n_hit endif !septum 1-2 ! write(21, '(a,i5,3x,i10,a,i5,a,i5,a)')'Turn ',i,n_clear,' clear septum 2', n_lost,' at septum 2', n_hit,' hit septum 2' ! track particles that clear septum 2 to start of ext branch if(n_clear >0)then call write_phase_space(lun_sep2_beam,i,ext_beam) ! call track_beam(lat,ext_beam,ele1=lat%branch(1)%ele(ix_magsept2(1)), ele2=lat%branch(1)%ele(lat%branch(1)%n_ele_track), err=err) !track from septum 2 to end of extraction line in ext branch do j=ix_magsept2(1),lat%branch(1)%n_ele_track-2 call track_beam(lat,ext_beam,ele1=lat%branch(1)%ele(j), ele2=lat%branch(1)%ele(j+1), err=err) !track to from septum 1 to 2 in ext branch do k=1,size(ext_beam%bunch(1)%particle(:)) if(ext_beam%bunch(1)%particle(k)%state == alive$) & write(lun_sep2_phasespace,'(i10,1x,a12,es12.4,1x,i10,1x,6es12.4,1x,i10)')j+1,lat%branch(1)%ele(j+1)%name,& lat%branch(1)%ele(j+1)%s,k, ext_beam%bunch(1)%particle(k)%vec, ext_beam%bunch(1)%particle(k)%state if(j+2 > lat%branch(1)%n_ele_track)then print *,' ramp_test: 538, j+2 = ', j+2, lat%branch(1)%n_ele_track endif if(lat%branch(1)%ele(j+2)%name == 'RING_TO_EXT' .and. ext_beam%bunch(1)%particle(k)%state == alive$) & write(lun_ext_phasespace,'(a12,es12.4,1x,i10,1x,6es12.4,1x,i10)')lat%branch(1)%ele(j+1)%name,& lat%branch(1)%ele(j+1)%s,k, ext_beam%bunch(1)%particle(k)%vec, ext_beam%bunch(1)%particle(k)%state end do end do call write_phase_space(lun_target, i, ext_beam) !particles that get to the target do k=1,size(ext_beam%bunch(1)%particle(:)) if(ext_beam%bunch(1)%particle(k)%state == alive$)then beam_on_target%bunch(1)%particle(k) = ext_beam%bunch(1)%particle(k) !copy particles that make it to target to beam_on_target ext_beam%bunch(1)%particle(k)%state = end_of_line$ print '(i10,1x,6es12.4)',ext_beam%bunch(1)%particle(k)%state,ext_beam%bunch(1)%particle(k)%vec endif end do endif endif !extracted particles call beam_moments(beam, bunch_moments) write(lun_beam_trajectory,'(i10,12es12.4,1x,i10,1x,2es12.4,1x,3i10,1x,3i10,1x,3i10)')i,bunch_moments(1)%ave(1:6), bunch_moments(1)%sigma(1,1), bunch_moments(1)%sigma(2,2), & bunch_moments(1)%sigma(3,3),bunch_moments(1)%sigma(4,4),bunch_moments(1)%sigma(5,5),bunch_moments(1)%sigma(6,6), & bunch_moments(1)%n_alive, bunch_moments(1)%ave_t, bunch_moments(1)%sigma_t, & n_extract_tot(0), n_lost_tot(0), n_sep_tot(0),n_extract_tot(1), n_lost_tot(1), n_sep_tot(1),n_extract_tot(2), n_lost_tot(2), n_sep_tot(2) co_fft(i)%vec = bunch_moments(1)%ave end do call beam_moments(beam_on_target, bunch_moments) write(lun_target_moments,'(i10,1x,6es12.4,1x,6es12.4,1x,2es12.4)')bunch_moments(1)%n_alive, bunch_moments(1)%ave(1:6),& bunch_moments(1)%sigma(1,1), bunch_moments(1)%sigma(2,2), bunch_moments(1)%sigma(3,3), bunch_moments(1)%sigma(4,4),& bunch_moments(1)%sigma(5,5), bunch_moments(1)%sigma(6,6), bunch_moments(1)%ave_t, bunch_moments(1)%sigma_t endif ! multiple particles ! compute fft for n_fft turns, at n_fft/2 intervals ix=1 lun=lunget() open(unit=lun,file =trim(directory)//'/tune_vs_turn.dat') write(lun,'(a10,a10,6a12)')'start turn','end turn','Qx','Qy','Qz','Phase x','Phase y','Phase z' do while(ix + n_fft <= nturns) coa(1:n_fft) = co_fft(ix:ix+n_fft-1) call fft_plane(n_fft,coa,Q,phase) write(lun,'(i10,i10, 3es12.4,3es12.4)')ix,ix+n_fft-1, Q, phase ! print '(i5,a1,i5, a,3es12.4,a,3es12.4)',ix,'-',ix+n_fft-1, ' Q = ', Q, ' phase = ', phase ix = ix+n_fft/2 end do close(lun) print '(a)','Columns in file: beam_trajectory.dat' print '(6a23)','1 = turn','2-7 = ','8-13 = ','14= n_alive','15 =