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), target :: lat type (lat_struct) lat type (coord_struct), allocatable::co(:), co_fft(:), coa(:), co_many(:) type (coord_struct) co_start type (beam_init_struct) beam_init type (beam_struct) beam type (moment_struct), allocatable :: bunch_moments(:) type (ele_pointer_struct), allocatable :: eles(:) integer ix, nargs integer i, j, k integer lun integer lun1, lun2, lun3 integer lun4, lun5 integer nturns integer ios integer n_fft/256/ integer nparticles integer n_extract, n_lost, n_sep integer n_extract_tot/0/, n_lost_tot/0/, n_sep_tot/0/ integer Nbin(-200:200), nn integer n_loc integer track_state integer ix_inj, ix_essep, ix_l3 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./ character*140 lat_file character*120 line, lineout character*10 iteration logical err logical ramp_synch, phase_space/.false./ logical sext/.false./ logical aperture_limit_on/.false./ logical scale_emit/.false./ namelist/input/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, & !extraction septum - center line to far edge of septum (0:2), septum thickness (0:2) x_wall, & !center line to wall 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, & ! hbump(1:3) in radians 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 bmad_com%auto_bookkeeper = .false. oct_ramp = -1. quad_ramp = -1. bump_ramp = -1. sext_ramp = -1. lun = lunget() 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) qv = qvert qh=qhorz t_ext = ext_time hbump(1:3) = bump_amp*hbump(1:3) !ready an input file from the command line 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) ' 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' print *, ' lat_file = ', lat_file endif 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) 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 writeinfo(lat) ! write initial twiss parameters to file lun1 = lunget() open(unit=lun1, file = 'trajectory_ele_by_ele.dat') lun2 = lunget() open(unit=lun2, file = 'trajectory.dat') lun3 = lunget() open(unit=lun3, file = 'beam_trajectory.dat') lun4 = lunget() if(phase_space)open(unit=lun4, file = 'phasespace.dat') lun5 = lunget() open(unit=lun5,file = 'extracted_beam.dat') do i=1,lat%n_ele_max if(lat%ele(i)%key == rbend$ .or. lat%ele(i)%key == quadrupole$ .or.lat%ele(i)%key == sbend$)then ! save reference values for k1, k2 lat%ele(i)%value(custom_attribute1$) = lat%ele(i)%value(k1$) lat%ele(i)%value(custom_attribute2$) = lat%ele(i)%value(k2$) endif if(lat%ele(i)%key == multipole$)then !save reference values of multipoles allocate(lat%ele(i)%r(1:2,0:n_pole_maxx,0:1)) lat%ele(i)%r(1,0:n_pole_maxx,1) = lat%ele(i)%a_pole(:) lat%ele(i)%r(2,0:n_pole_maxx,1) = lat%ele(i)%b_pole(:) endif if(index(lat%ele(i)%name,'RES_EXT_SEXT')/=0)then lat%ele(i)%value(custom_attribute2$) = lat%ele(i)%value(k2$) print '(a11,a16,a7)',' Sextupole ',lat%ele(i)%name,' found ' sext = .true. endif if(lat%ele(i)%key == rfcavity$) lat%ele(i)%value(phi0$) = rf_phase_twopi/twopi ! set RF phase end do if(.not. sext)then print *,' No resonance driving sextupole found ' stop endif ! find location of Injection point (POSINJSEP), electrostatic septum (ESSEP), and BPM at l3 (L3) call lat_ele_locator('POSINJSEP',lat,eles,n_loc, err) 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) ix_essep = eles(1)%ele%ix_ele print '(a14,a16,a9,i10)',' Extraction at ',eles(1)%ele%name,' element ',ix_essep call lat_ele_locator('L3',lat,eles,n_loc,err) ix_l3 = eles(1)%ele%ix_ele print '(a14,a16,a9,i10)',' Injection at ',eles(1)%ele%name,' element ',ix_l3 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$) t0 = acos(1.-2*delta)/omega 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 endif print '(a,es12.4)',' Ramp initiated, co(0)%vec(6) = ', co(0)%vec(6) ! co_start%vec = co_start%vec + initial_offset 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(lun2,'(/,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 electrostatic septum do i=1,nturns write(lun1, '(/,a,a7,i10,/)')'#','turn = ', i call track_many(lat,co_many,ix_essep,ix_essep,direction=1,ix_branch=0, track_state=track_state) !track full turn at electrostatic septum co_fft(i) = co_many(ix_l3) write(lun2, '(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(lun2,'(a23,1x,i5,1x,a10,a16)')'Particle lost at index ',track_state,' element =',lat%ele(track_state)%name if(track_state /= moving_forward$)exit if((i/10)*10 == i)call write_full_turn(lat,i,co_many, lat%n_ele_track) end do co_start%vec(1) = co_start%vec(1)+dx end do endif ! single particle 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) ! 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(lun3,'(a10,2a72,1x,a10,1x,2a12,1x,2a10,1x,2a10,1x,2a10)')'turn','bunch_moments%ave(1:6)', 'bunch_moments%sigma(1:6,1:6)','n_alive','ave time','sigma_t', & 'n_extract','n_lost','n_extract_tot','n_lost_tot','n_sep','n_sep_tot' if(phase_space)write(lun4, '(6a12)') 'x','xp','y','yp','z','zp' write(lun5,'(2a10,1x,6a12)')'turn','particle','x','xp','y','yp','z','zp' bmad_com%aperture_limit_on = aperture_limit_on i=0 if(phase_space)then write(lun4, '(/,a1,a7,i10,/)')'#', 'turn = ',i do j=1,size(beam%bunch(1)%particle(:)) if(beam%bunch(1)%particle(j)%state == alive$)write(lun4, '(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 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(lun4, '(/,a1,a7,i10,/)')'#', 'turn = ',i do j=1,size(beam%bunch(1)%particle(:)) if(beam%bunch(1)%particle(j)%state == alive$)write(lun4, '(6es12.4)')beam%bunch(1)%particle(j)%vec end do endif n_extract=0 n_lost = 0 n_sep = 0 Nbin(:)=0 if(i > ext_time/2.52131e-6)then do j=1,size(beam%bunch(1)%particle(:)) if(beam%bunch(1)%particle(j)%state == alive$)then x=beam%bunch(1)%particle(j)%vec(1) if(x > x_septum(0) .and. x < x_wall)then n_extract = n_extract+1 write(lun5,'(2i10,1x,6es12.4)')i,j,beam%bunch(1)%particle(j)%vec endif if(x > x_wall) n_lost = n_lost+1 if(x > x_septum(0))beam%bunch(1)%particle(j)%state = lost$ if(x < x_septum(0) .and. x > x_septum(0) - t_septum(0))n_sep = n_sep + 1 if(x < x_septum(0) .and. x > x_septum(0) - t_septum(0))beam%bunch(1)%particle(j)%state = lost$ nn = min(x/width+sign(x,0.5_rp),0.02/width) Nbin(nn) = Nbin(nn)+1 endif end do n_extract_tot = n_extract_tot + n_extract n_lost_tot = n_lost_tot + n_lost n_sep_tot = n_sep_tot + n_sep endif ! pause call beam_moments(beam, bunch_moments) write(lun3,'(i10,12es12.4,1x,i10,1x,2es12.4,1x,2i10,1x,2i10,1x,2i10)')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, n_lost, & n_extract_tot, n_lost_tot, n_sep, n_sep_tot co_fft(i)%vec = bunch_moments(1)%ave write(27,'(/,a1,/)')'#' do nn=-200,200 write(27,'(3i10)')i,nn,Nbin(nn) end do end do endif ! multiple particles ! compute fft for n_fft turns, at n_fft/2 intervals ix=1 write(13,'(a10,a10,a36,a36)')'start turn','end turn','Q','Phase' 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(13,'(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 print '(a)','Columns in file: beam_trajectory.dat' print '(6a23)','1 = turn','2-7 = ','8-13 = ','14= n_alive','15 =