program envelope ! !=================================================================== ! Implemented tune setting ! 20 Sep 2006 ! !=================================================================== ! Major modifications with goal of adding long-range beam-beam ! interaction simulation. ! 18 May 2004 !=================================================================== ! Changed ring%param%E_TOT to ring%E_TOT ! 14 May 2004 !=================================================================== ! Removed all uses of zero_real_vec, which didn't work for double ! precision vectors! ! JAC 25 July 2003 !=================================================================== ! Removed use of n_ele_maxx as per DCS e-mail of 7/9/03 ! Added dynamic allocation of coord_struct arrays in main program ! Removed calls to TRACK_LONG and the matrices mats627 ! JAC 22 July 2003 !=================================================================== ! This program calculates aperture envelopes for CESR injection ! in routine ENVELOPE_NT ! Adapted from INJTRACK.F90. ! ! J.A. Crittenden 4/14/03 !=================================================================== ! This is the CESR injection simulation code developed by S.Henderson. ! It is described in the note of 6/1/01 "REU Injection Simulation Project" ! which was worked on during the summer of 2001 by Morris. ! Further references to CESR injection procedures can be found in ! "A Description of CESR Injection" of 10/27/95 by SDH, and in the ! Ph.D. thesis of John Seeman of 1979. ! ! J.A.Crittenden 10/25/01 !=================================================================== use bmad use bmadz_mod use bmadz_interface use z_tune_mod use injlib use envelope_beambeam_setup_mod use injparam_mod implicit none type (lat_struct) ring, ring_in, ring_out, ring_test type (ele_struct), pointer :: ele type (coord_struct), allocatable :: inj_orbit(:) type (coord_struct), allocatable :: closed_orbit(:) type (coord_struct), allocatable, save :: co(:) type (coord_struct), allocatable :: diff_co(:) type (coord_struct), allocatable :: prz_orbit(:) real(rdef), allocatable :: dk1(:) real(rdef) phi_a, phi_b integer int_Q_x, int_Q_y integer nquad logical ok type (inj_struct) inj type (synch_beam_struct) synch_beam type (pretzel_struct) pretzel type (qtune_struct) qtune type (xqune_struct) xqune type (tunes_struct) tunes type (coord_struct) delta_ip integer i, j, n, ix_start, tans integer imax integer indprint integer norbwrt1,norbwrt2 integer anal, species(2) integer lun, nturn, n_part integer n_trains_strong, n_cars_strong, train_weak, car_weak integer species_strong integer coord integer lun5 integer pick, try, state integer analysis_code integer set_lrbbi real qvt, qht, qv(2), qh(2) real injang, angmin, angmax, angstep, angmid real eff, qhmin, qhmax, qhstep, qvmin, qvmax, qvstep real x0inj, x0pinj, nsig, nsigsave, weight, totweight, injweight real i_bunch real injpos, posmin, posmax, posstep character this_lat*100, latname*40, qtname*100 character lat_file*80 character fname*40 character injection_file*40 character lrbbi_file*40 character tunescan_file*40 character phasespace_file*40 integer set_ccu /0/ real(rdef) wt_lrbbi /1./ logical on character datfname*40 integer newid character*4 apname(4)/'entr','exit','both','none'/ character date*10, time*10 integer version, track_state integer strong_particle, i_train, j_car, n_trains_tot, n_cars, slices real(rdef) bbiwt(9,10) /90*0./ logical rec_taylor, beambeam_ip, custom_pat real(rdef) current, coupling_sb logical lostit, use_lrbbi, first_loss logical use_coupling, error real(rdef) qx,qy,qz logical constant_tunes /.false./ ! !======================================= ! PAW storage real hm(30000000) common/pawc/hm !======================================= namelist / lrbbi_input / strong_particle, i_train, j_car, n_trains_tot, n_cars, current, bbiwt, rec_taylor, & beambeam_ip, coupling_sb, slices, custom_pat ! namelist / envelope_setup / this_lat, analysis_code, set_lrbbi, use_coupling, & qx, qy, qz, constant_tunes, & injection_file, lrbbi_file, tunescan_file, phasespace_file, & set_ccu, wt_lrbbi ! call date_and_time(date, time) write(6,11)date,time 11 format(' ==> Initializing Program ENVELOPE on ',a8,' at ',a10) !======================== ! Initialize HBOOK storage ! call hlimit(30000000) ! !====================================================================== ! Program setup parameters !====================================================================== ! Read in ENVELOPE program setup info lun = lunget() fname='envelope_setup.in' open(lun,file=fname, status= 'old') read(lun, nml=envelope_setup) close(lun) ! write(6,1000)fname, this_lat, analysis_code, set_lrbbi, use_coupling, & qx, qy, qz, constant_tunes, & injection_file, lrbbi_file, tunescan_file, set_ccu, wt_lrbbi 1000 format(' ==> Program ENVELOPE initialized with file ',a40/, & ' The chosen lattice is: ',a40/, & ' The analysis code is:'i4/, & ' The LRBBI setup choice is:',i4/, & ' The coupling setup choice is:',l/, & ' The chosen tune values are: ',3f7.3/, & ' CONSTANT_TUNES is set to: ',l1/, & ' The injection setup file name is: ',a40/, & ' The lrbbi setup file name is: ',a40/, & ' The tunescan setup file name is: ',a40/, & ' SET_CCU is set to:',i4/, & ' The lrbbi weight is set to:',f7.3) ! ! get lattice. Setting 5th argument to lattice name prevents interactive query ! call choose_cesr_lattice( latname, lat_file, this_lat, ring, this_lat ) ! qtname='/cdat/cesr38/disk1/critten/cesr/lat/bmad_'//this_lat ! qtname='BMAD_LAT:bmad_'//this_lat ! qtname='BMAD_LAT:'//this_lat ! qtname= '/home/sw565/chess_undulator/'//this_lat qtname = this_lat print *,'Lattice: ',qtname ! print *,' Loading lattice from file ',qtname call bmad_parser(qtname, ring) ! print *,' After loading lattice:' ! write(6,1111)(ring%ele(i)%name,ring%ele(i)%s,ring%ele(i)%a%beta,i=1,ring%n_ele_track) 1111 format(4(a10,2x,2f8.2,5x)) ele => ring%ele(720) ! if (.not. associated(ele%taylor(1)%term)) then ! print *,' Taylor term is not associated after bmad_parser' ! else ! print *,' Taylor term is associated after bmad_parser' ! endif ! Turn off compact undulators if SET_CCU=-1 ! Turn on compact undulators if SET_CCU=1 ! Otherwise do nothing if( set_ccu.eq.1 ) then on = .true. call turn_ccu_onoff(ring, on) elseif( set_ccu.eq.-1 ) then on = .false. call turn_ccu_onoff(ring, on) endif ! Allocate storage call reallocate_coord(co,ring%n_ele_max) call reallocate_coord(diff_co,ring%n_ele_max) ! Read in QTUNE knobs (RAW_QTUNE_5 and RAW_QTUNE_6) ! qtname='[cesr.bmad.lat]BMAD_QTUNES_K9A18A000.FD96E_4S_15KG_INJ' ! qtname='BMAD_LAT:BMAD_QTUNES_K9A18A000.FD96E_4S_15KG_INJ.LAT' ! qtname='/cdat/cesr38/disk1/critten/cesr/lat/bmad_qtunes_k9a18a000.fd96e_4s_15kg_inj' ! call bmad_parser2(qtname, ring) ! print *,' After loading qtunes:',(i,ring%ele(i)%name,i=1,ring%n_ele_track) ! ! write (6,1550)(i,ring%ele(i)%name, & ! i=1,ring%n_ele_track) 1550 format(' After loading lattice file:'/ & 500(8(i5,1x,a8,4x)/)) ! Read in aperture limit definitions ! qtname='/cdat/cesr38/disk1/critten/cesr/layout/cesr_c1_aperture_bends_nowigs.bmad' ! qtname='/cdat/cesr38/disk1/critten/cesr/layout/cesr_c1_aperture_bends_1wig.bmad' ! qtname='/cdat/cesr38/disk1/critten/cesr/layout/cesr_c1_aperture_bends_12wig.bmad' ! Include apertures for CHESS wigglers ! qtname='/cdat/cesr38/disk1/critten/cesr/layout/cesr_c1_aperture_bends_14wig.bmad' ! West CHESS wiggler vertical aperture changed from 2.5 to 1.9 cm ! This file in use until 27 June 2007 13:45 ! qtname='/cdat/cesr38/disk1/critten/cesr/layout/cesr_c1_aperture_bends_14wiga.bmad' ! Reduce Q28E aperture from 45 to 37 mm ! Began useing this file 27 Juen 2007 13:45 ! qtname='/cdat/cesr38/disk1/critten/cesr/layout/cesr_c1_aperture_bends_14wiga_28e.bmad' ! qtname='BMAD_LAYOUT:cesr_c1_aperture_bends_14wiga_28e.bmad' ! qtname='/cdat/cesr38/disk1/critten/cesr/layout/cesr_c1_aperture_bends_10wig.bmad' ! ! 4 Sep 09 ! Reduce Q28E aperture from 45 to 37 mm ! qtname='/cdat/cesr38/disk1/critten/cesr/layout/cesrta_aperture_28e.bmad' ! qtname='/nfs/acc/user/critten/cesr38/cesr/layout/cesrta_aperture_28e.bmad' qtname='/nfs/acc/user/critten/cesr38/cesr/layout/chess_aperture_28e.bmad' ! qtname='/nfs/acc/user/critten/cesr38/cesr/layout/chess_aperture.bmad' print *,' Aperture file is:',qtname call bmad_parser2(qtname, ring) call calc_z_tune (ring) ! write (6,1600)(ring%ele(i)%name, ring%ele(i)%s, & ! ring%ele(i)%value(x1_limit$),ring%ele(i)%value(y1_limit$), apname(ring%ele(i)%aperture_at), & ! i=1,ring%n_ele_track) 1600 format(' After loading aperture limits:'/ & 500(3(a8,' s: ',f6.2,' X: ',f6.3,' Y: ',f6.3,1x,a4,1x)/)) ! Initialize long-range beam-beam interaction simulation if( set_lrbbi == 1 ) then use_lrbbi = .true. lun = lunget() open(lun,file=lrbbi_file, status= 'old') read(lun, nml=lrbbi_input) close(lun) write(6,1050)lrbbi_file, strong_particle, i_train, j_car, n_trains_tot, n_cars, bbiwt, current, rec_taylor, & beambeam_ip, coupling_sb, slices,custom_pat 1050 format(' ==> LRBBI initialized with file ',a40/, & ' Strong beam particle type (-1:e- +1:e+): ',i4/, & ' Test train:',i4/, & ' Test car:'i4/, & ' Total nr of trains:',i4/, & ' Nr of cars/train: ',i4/, & ' BBI weights: ',9f10.2/,9(18x,9f10.2/), & ' Current/bunch: ',f4.1/, & ' Taylor tracking: ',l/, & ' Include BBI at IP: ',l/, & ' Coupling: ',f6.3/, & ' Nr of slices: ',i4/, & ' Custom Bunch pattern: ',l) ! If BBIWT not present, set it to be evenly weighted ! over the number of cars and trains if ( all (bbiwt .eq. 0) ) then print *, ' ==> LRBBI initialization found all BBIWT=0. Setting uniform bunch currents.' do i=1,9 do j=1,5 bbiwt(i,j)=1. enddo enddo endif ! For purposes of setting up the LRBBI and IP BBI, turn off ! radiation fluctuations and damping and turn off the RF bmad_com%radiation_damping_on = .false. bmad_com%radiation_fluctuations_on = .false. call set_on_off (rfcavity$, ring, off$) ! Calculate and print tunes for electrons as bbi are added ring%param%particle = electron$ ! call twiss_at_start(ring) ! co(0)%vec = 0. ! call closed_orbit_at_start(ring, co(0), 4, .true.) ! call closed_orbit_calc(ring, co, 4) ! call track_all (ring, co) ! call lat_make_mat6(ring,-1,co) ! call twiss_at_start(ring) ! call twiss_propagate_all(ring) ! Do not stop program if error in tracking_update global_com%exit_on_error = .false. bmad_com%max_aperture_limit=1. call tracking_update(ring, co, 4, diff_co) !************************************************************ ! Set desired tunes with no BBI ! Calculate phi_a and phi_b qht = ring%a%tune / twopi qvt = ring%b%tune / twopi if ( qx .ne. 0. ) qht = qx if ( qy .ne. 0. ) qvt = qy int_Q_x = int(ring%ele(ring%n_ele_track)%a%phi / twopi) int_Q_y = int(ring%ele(ring%n_ele_track)%b%phi / twopi) phi_a = ( int_Q_x + qht ) * twopi phi_b = ( int_Q_y + qvt ) * twopi ! ! If not using design tunes, set them ! allocate(dk1(ring%n_ele_max),stat=j) call re_allocate(dk1,ring%n_ele_max) if ( qx .ne. 0. .or. qy .ne. 0. ) then call choose_quads( ring, dk1 ) nquad = 0 ! print *, 'Rtn CHOOSE_QUADS chose the following quads:' do i=1,ring%n_ele_track if ( dk1(i) .ne. 0 ) then nquad = nquad + 1 ! write ( 6, 5000 ) nquad, ring%ele( i )%name, dk1(i) !5000 format(i5,2x,a15,3x,'dk = ',f6.1) endif enddo print *,' Call custom_set_tune with phi_a, phi_b:',phi_a,phi_b call custom_set_tune (phi_a, phi_b, dk1, ring, co, ok) deallocate(dk1) if (.not. ok ) then print *,' ERROR from CUSTOM_SET_TUNE! qht,qvt=',qht,qvt stop endif endif if (qz .ne. 0) call set_z_tune(ring, qz*twopi) type * type *,' Electron tunes after setting tunes with no BBI:' type *,' Qx = ',ring%a%tune/twopi,' Qy = ',ring%b%tune/twopi,' Qz = ',ring%z%tune/twopi type '(a15,6e12.4)',' Closed orbit ', co(0)%vec(1:6) ! End of tune setting !************************************************************ type * type *,' Electron tunes before adding BBI:' type *,' Qx = ',qht,' Qy = ',qvt type '(a15,4e12.4)',' Closed orbit ', co(0)%vec(1:4) ! LRBBI_SETUP adds beam-beam elements to the ring. It splits elements if necessary, including wigglers.) It returns the ring to be used for tracking electons. print *,' LRBBI_SETUP called for strong-beam particle type',strong_particle ring_in = ring call lrbbi_setup ( ring_in, ring_out, strong_particle, i_train, j_car, n_trains_tot, n_cars, current, rec_taylor, bbiwt, custom_pat) ring = ring_out ! Reallocate storage because ring size changed call reallocate_coord(co,ring%n_ele_max) call reallocate_coord(diff_co,ring%n_ele_max) ring%param%particle = electron$ call tracking_update(ring, co, 4, diff_co) type * type '(a20,f4.1,a8)', & ' After LRBBI_SETUP ', current,'mA/bunch' write(6,2070)ring%a%tune/twopi,ring%b%tune/twopi,ring%z%tune/twopi 2070 format('Electron tunes: Qx = ',f6.3,' Qy = ',f6.3, ' Qz = ',f6.3) type '(a15,4e12.4)',' Closed orbit ', co(0)%vec(1:4) ! Initialize beam-beam interaction at interaction point and print out tunes if(beambeam_ip)then ! ! ENVELOPE_BEAMBEAM_SETUP needs the sync tune (energy spread for beam size at IP) ! ring%z%tune = -0.089 * twopi ! Set up the ring for positrons as expected by ENVELOPE_BEAMBEAM_SETUP ! ring%param%particle = positron$ ! call tracking_update(ring, co, 4, diff_co) ! ring%param%n_part = 0. ! call twiss_at_start(ring) ! co(0)%vec = 0. ! call closed_orbit_at_start(ring, co(0), 4, .true.) ! call closed_orbit_calc(ring, co, 4) ! call track_all (ring, co) ! call lat_make_mat6(ring,-1, co) ! call twiss_at_start(ring) ! call twiss_propagate_all(ring) print *,'part,curr,coup,sl=',strong_particle, current, coupling_sb, slices call envelope_beambeam_setup(ring, strong_particle, current, 1.0_rp, coupling_sb, slices) print *,' ring%param%n_part=', ring%param%n_part ! Reallocate storage because ring size changed call reallocate_coord(co,ring%n_ele_max) call reallocate_coord(diff_co,ring%n_ele_max) ring%param%particle = electron$ call tracking_update(ring, co, 4, diff_co) type '(a52,4e12.4)',' Electron closed orbit at IP after creating IP BBI ', co(0)%vec(1:4) ! Type out tunes write(6,2100)ring%a%tune/twopi,ring%b%tune/twopi,ring%z%tune/twopi 2100 format('Electron tunes after adding IP BBI: Qx = ',f6.3,' Qy = ',f6.3,' Qz = ',f6.3) ! TEST_BBI used to make sure IP BBI and LRBBI are working. ! ring_test = ring ! call test_bbi(ring_test, current) endif ! Finished adding IP BBI else ! SET_LRBBI =0. Weight any pre-existing beam-beam element strengths in the input lattice. if (wt_lrbbi .eq. 0.) then use_lrbbi = .false. call weight_bbi(ring, wt_lrbbi) elseif (wt_lrbbi .ne. 1.) then call weight_bbi(ring, wt_lrbbi) use_lrbbi = .true. else use_lrbbi = .true. endif endif ! write (6,1700)(i,ring%ele(i)%name, ring%ele(i)%s, & 100*ring%ele(i)%value(x1_limit$), & 100*ring%ele(i)%value(y1_limit$), & apname(ring%ele(i)%aperture_at), & i=1,ring%n_ele_max) 1700 format(' Ring after creating or weighting BBI:'/ & 500(1x,3(i4,1x,a15,' s: ',f6.2,' X: ',f4.1,' Y: ',f4.1,1x,a4,x)/)) !************************************************************ ! Set desired tunes if ( constant_tunes ) then print *,' Call custom_set_tune with phi_a, phi_b:',phi_a,phi_b ! allocate(dk1(ring%n_ele_max),stat=j) call re_allocate(dk1,ring%n_ele_max) call choose_quads( ring, dk1 ) call custom_set_tune (phi_a, phi_b, dk1, ring, co, ok) deallocate(dk1) if (.not. ok ) then write (6,2900) qht,qvt,current 2900 format(' Error from CUSTOM_SET_TUNE for qht,qvt =',2f6.3,& ' current=',f5.2,'-- skip this point.') goto 9999 endif type * type *,' Electron tunes after setting tunes:' type *,' Qx = ',ring%a%tune/twopi,' Qy = ',ring%b%tune/twopi,' Qz = ',ring%z%tune/twopi type '(a15,6e12.4)',' Closed orbit ', co(0)%vec(1:6) endif ! End of tune setting !************************************************************ ! ! Set bumps call set_bumps( ring, 4 ) call tracking_update(ring, co, 4, diff_co) ! Ensure beams are colliding call reallocate_coord(prz_orbit,ring%n_ele_max) call init_pretzel(ring,pretzel,prz_orbit) ! Do not collide CHESS beams 10 sep 09 jac ! call collide ( ring, 4 ) call tracking_update(ring, co, 4, diff_co) write(6,1850)ring%a%tune/twopi,ring%b%tune/twopi,ring%z%tune/twopi 2001 format('e- tunes prior to calling INIT_INJECTION: Qx = ',f6.3, ' Qy = ',f6.3, ' Qz = ',f6.3) !************************************************************ call init_injection(ring, pretzel, inj, synch_beam, & injection_file, use_lrbbi ) ring%param%particle = electron$ call tracking_update(ring, co, 4, diff_co) !************************************************************ ! Set desired tunes if ( constant_tunes ) then print *,' Call custom_set_tune with phi_a, phi_b:',phi_a,phi_b ! allocate(dk1(ring%n_ele_max),stat=j) call re_allocate(dk1,ring%n_ele_max) call choose_quads( ring, dk1 ) call custom_set_tune (phi_a, phi_b, dk1, ring, co, ok) deallocate(dk1) if (.not. ok ) then write (6,2900) qht,qvt,current goto 9999 endif type * type *,' Electron tunes after setting tunes:' type *,' Qx = ',ring%a%tune/twopi,' Qy = ',ring%b%tune/twopi,' Qz = ',ring%z%tune/twopi type '(a15,6e12.4)',' Closed orbit ', co(0)%vec(1:6) endif ! End of tune setting !************************************************************ call reallocate_coord(closed_orbit,ring%n_ele_max) call tracking_update(ring, closed_orbit, 4, diff_co) print *,' Horiz Beta and alpha functions after INIT_INJECTION:' write(6,2120)(i,ring%ele(i)%name,ring%ele(i)%s,ring%ele(i)%a%beta,ring%ele(i)%a%alpha,i=1,ring%n_ele_max) 2120 format(4(i4,1x,a15,1x,3f5.1,2x)) write(6,2150)ring%a%tune/twopi,ring%b%tune/twopi 2150 format('Tunes after initializing injection: Qx = ',f6.3, ' Qy = ',f6.3) ! ! ! ------------------------------------------------------------------------- ! Book, fill and write out ntuple containing injection envelope information ! ------------------------------------------------------------------------- call envelope_nt(ring, pretzel, inj, synch_beam, this_lat, use_coupling, n_trains_tot, n_cars, set_lrbbi, use_lrbbi, custom_pat) ! For the following analyses, ! set up six-dimensional tracking including longitudinal coordinates ! Include damping effects, but exclude fluctuations bmad_com%radiation_damping_on = .true. bmad_com%radiation_fluctuations_on = .false. ! Turn on RF cavities print *,'Before turning on RF ring%z%tune is',ring%z%tune call set_on_off (rfcavity$, ring, on$) print *,'After turning on RF ring%z%tune is',ring%z%tune ! Remove ad hoc tune setting ! 18 Feb 2015 jac ! ring%z%tune = -0.089 * twopi ! call set_z_tune(ring) !call twiss_at_start(ring) !co(0)%vec=0. !call closed_orbit_calc(ring, co, 6) ! call closed_orbit_at_start(ring, co(0), 6, .true.) ! call track_all(ring, co) !call lat_make_mat6(ring, -1, co) !call twiss_at_start(ring) ! New simplified routine added 7/2/04 by DTK call tracking_update(ring, co, 6, diff_co, error) ! Chance to quit out before all hell breaks loose if(error) then print *,'There were errors encountered during tracking and execution of the program will be halted.' goto 9999 endif write(6,1850)ring%a%tune/twopi,ring%b%tune/twopi,ring%z%tune/twopi 1850 format('e- tunes after turning on RF: Qx = ',f6.3, ' Qy = ',f6.3, ' Qz = ',f6.3) ! Calculate beam-beam separations at all elements (delta_ip) ! and print out separation at IP ! Be sure particle type is electron (weak beam) when calling beambeam_separation ring%param%particle = electron$ call beambeam_separation ( ring, delta_ip, 6 ) ! ============================================ ! Begin branching over chosen analysis type ! ============================================ ! !=========================================================== ! The analysis code is ! the sum of 2 raised to the powers N-1, where N are ! the desired analysis numbers. ! For example, an analysis code value of 3 will result ! in the analyses of types 1 and 2 to be carried out. ! 1 --> Track the stored e- beam through pulsed bump' ! 2 --> Track the injected e- beam through 10 turns' ! 3 --> Test tune setting' ! 4 --> Scan injection angle (obsolete as of 15 October 2003)' ! 5 --> Injection Tune scan' ! 6 --> Track a phase space distribution' ! 7 --> Test random phase space' ! 8 --> Find efficiency for given injection coordinates' ! 9 --> Find injection efficiency for pinger kicks' ! (Requires interactive input) ! 10 --> Scan injection positions and angles for' ! something that survives NTESTTURNS turns' !=========================================================== anal=analysis_code if( and(anal,1) .ne. 0 ) then ! track stored e- beam about pulsed bump write(6,1100) 1100 format(/' -------- Begin Analysis Type 1 ----------') ! Write file of stored electron orbits without bump or pinger lun = lunget() write(datfname,1198) 1198 format('dat/stored_ele_orbits_nopb.dat') open(lun, file=datfname, status = 'unknown') call reallocate_coord(inj_orbit,ring%n_ele_max) inj_orbit(0)%vec(:) = 0 ! call closed_orbit_at_start (ring, inj_orbit(0), 4, .true.) call closed_orbit_calc (ring, inj_orbit, 4 ) ! do n=-5,5 n=-5 track_state = moving_forward$ inj_orbit(0)%vec(2) = -inj_orbit(0)%vec(2) inj_orbit(0)%vec(4) = -inj_orbit(0)%vec(4) inj_orbit(0)%vec(5) = -inj_orbit(0)%vec(5) do while ( n .le. 15 .and. track_state == moving_forward$) ! call pulsed_bump(ring, n, inj ) ! bmad_com%aperture_limit_on=.true. ring%param%particle = antiparticle(ring%param%particle) call track_many (ring, inj_orbit, 0, 0, -1, track_state = track_state) ring%param%particle = antiparticle(ring%param%particle) if (track_state /= moving_forward$) then imax = track_state else imax = ring%n_ele_track endif do i= imax,0,-1 write(lun,100) n, ring%ele(i)%s, & 1000.*inj_orbit(i)%vec(1),1000.*inj_orbit(i)%vec(2), & 1000.*inj_orbit(i)%vec(3),1000.*inj_orbit(i)%vec(4), & 1000.*inj_orbit(i)%vec(5),1000.*inj_orbit(i)%vec(6) 100 format(2x, i6, f10.4, 2(1x,f8.3),2(1x,f10.5),1x,f8.3,1x,f11.6) enddo n = n + 1 enddo close(lun) ! Write file of stored electron orbits including bump, pinger and pretzel lun = lunget() write(datfname,1200) 1200 format('dat/stored_ele_orbits.dat') open(lun, file=datfname, status = 'unknown') inj_orbit(0)%vec(:) = 0. ! call closed_orbit_at_start (ring, inj_orbit(0), 4, .true.) call closed_orbit_calc (ring, inj_orbit, 4) ! We are tracking electrons backwards, so change sign of angles and time ! 16 Dec 2014 jac & dcs inj_orbit(0)%vec(2) = -inj_orbit(0)%vec(2) ! x-prime inj_orbit(0)%vec(4) = -inj_orbit(0)%vec(4) ! y-prime inj_orbit(0)%vec(5) = -inj_orbit(0)%vec(5) ! time ! do n=-5,5 n=-5 track_state = moving_forward$ do while ( n .le. 15 .and. track_state == moving_forward$) call pulsed_bump(ring, n, inj ) ! if(n.lt.6) then if(n.eq.0) then print *,' Stored-electron orbit calculation calling PRINT_BUMPING for turn ',n call print_bumping(ring, n, inj ) endif ! bmad_com%aperture_limit_on=.true. ring%param%particle = antiparticle(ring%param%particle) call track_many (ring, inj_orbit, 0, 0, -1, track_state = track_state) ring%param%particle = antiparticle(ring%param%particle) ! print *,' At IP on turn',n,': ',inj_orbit(0)%vec(1), inj_orbit(0)%vec(2) if (track_state /= moving_forward$) then imax = track_state else imax = ring%n_ele_track endif do i = imax, 0, -1 write(lun,100) n, ring%ele(i)%s, & 1000.*inj_orbit(i)%vec(1),1000.*inj_orbit(i)%vec(2), & 1000.*inj_orbit(i)%vec(3),1000.*inj_orbit(i)%vec(4), & 1000.*inj_orbit(i)%vec(5),1000.*inj_orbit(i)%vec(6) enddo n = n + 1 enddo close(lun) endif if ( and(anal,2) .ne. 0 ) then ! Track injected e- beam for 10 turns write(6, '(/, a)') '-------- Begin Analysis Type 2 ----------' lun = lunget() datfname = 'dat/inj_ele_orbits.dat' open(lun, file=datfname, status = 'unknown') call inject_particle( ring, .false., 0, 0., weight, inj, synch_beam, inj_orbit ) ! Ensure TRACK_MANY will check aperture limits ! It stops tracking if particle lost bmad_com%aperture_limit_on = .true. first_loss=.true. ! ring%ele(inj%ix_injpt)%value(x1_limit$)=0.10 ! ring%ele(inj%ix_injsex)%value(x1_limit$)=0.10 ! Be sure to include overlords between n_ele_use and n_ele_max, ! because they control the slaves which are created by the bbi elements do i=1,ring%n_ele_max if(ring%ele( i )%name(1:7).eq.'SEX_34E'.or.ring%ele( i )%name(1:4).eq.'Q34E'.or.ring%ele( i )%name(1:4).eq.'B35E')then ring%ele(i)%value(x1_limit$)=0.10 ring%ele(i)%value(x2_limit$)=0.10 ! type *,' Expanding aperture limit to 10 cm for element ',i,ring%ele(i)%name endif enddo ! Take the beam to the IP ix_start = inj%ix_injpt norbwrt1 = 0 norbwrt2 = inj%ix_injpt ! do n=0,10 ! ! Track from injection point (ix_start) to IP (ix=0) on turn 0 ! 16 April 2003 Since the last version of MANY_TRACK, the Z coordinate it expects to receive ! has changed sign. This means that xp,yp, and the time z need to have their sign changed ! before this first call. ! 10/7/2003 Moved this sign change from inside the loop over turns just before ! call to TRACK_MANY to outside the loop ! inj_orbit(ix_start)%vec(2) = - inj_orbit(ix_start)%vec(2) ! inj_orbit(ix_start)%vec(4) = - inj_orbit(ix_start)%vec(4) ! inj_orbit(ix_start)%vec(5) = - inj_orbit(ix_start)%vec(5) ! n=-1 track_state = moving_forward$ do while (n.le.15 .and. track_state == moving_forward$) n=n+1 call pulsed_bump(ring, n, inj ) call init_coord (inj_orbit(ix_start), inj_orbit(ix_start)%vec, ring%ele(ix_start), & downstream_end$, ring%param%particle, -1) ring%param%particle = antiparticle(ring%param%particle) call track_many (ring, inj_orbit, ix_start, 0, -1, track_state = track_state) ring%param%particle = antiparticle(ring%param%particle) ! Print out recent tracking history if lost ! Always print out tracking near injection point if (track_state /= moving_forward$ .or. n .eq. 0 ) then if (track_state /= moving_forward$) then norbwrt1 = track_state write(6,5000)n,ring%ele(track_state)%name 5000 format('Lost injected electron on turn ',i2,' at element ',a12) if (first_loss) then first_loss=.false. indprint = track_state print *, ' Some tracking history follows:' print *, & 'Element s kick k1 k2 orb-x orb-xp x-limit' else ! don't print anything if not first loss indprint=-1000 endif else ! Always print injection turn tracking around injection point print *, 'Electron survived injection turn' indprint=505 print *, ' Some tracking history follows:' print *, & 'Element s kick k1 k2 orb-x orb-xp x-limit' endif do i=indprint+190,indprint-5,-1 if(i.gt.0.and.i.le.ring%n_ele_track)then write(6,5100)i, & ring%ele( i )%name, & ring%ele( i )%s, & ring%ele( i )%value(angle$), & ring%ele( i )%value(k1$), & ring%ele( i )%value(k2$), & 100*inj_orbit( i )%vec(1), & 100*inj_orbit( i )%vec(2), & 100*ring%ele( i )%value(x1_limit$), & ring%ele(i)%aperture_at, & apname(ring%ele(i)%aperture_at) 5100 format(1x,i4,1x,a10,f8.2,3f8.3,f8.2,f8.3,f8.2,1x,i4,1x,a4) endif enddo endif ! print *,' At IP: ',inj_orbit(0)%vec(1), inj_orbit(0)%vec(2) do i=norbwrt1,norbwrt2 write(lun,101) n, ring%ele(i)%s, & 1000.*inj_orbit(i)%vec(1),1000.*inj_orbit(i)%vec(2), & 1000.*inj_orbit(i)%vec(3),1000.*inj_orbit(i)%vec(4), & 1000.*inj_orbit(i)%vec(5),1000.*inj_orbit(i)%vec(6) 101 format(2x, i6, f10.4, 2(1x,f8.3),2(1x,f10.5),1x,f8.3,1x,f11.6) enddo ! After injection turn, track from IP and write orbit ! to file for whole ring ix_start = 0 norbwrt2 = ring%n_ele_track ! After injection turn, set aperture near injection point back to ! beampipe wall position ! ring%ele(inj%ix_injpt)%value(x1_limit$)=0.045 ! ring%ele(inj%ix_injsex)%value(x1_limit$)=0.045 ! Be sure to include overlords between n_ele_use and n_ele_max, ! because they control the slaves which are created by the bbi elements do j=1,ring%n_ele_max if(ring%ele( j )%name(1:7).eq.'SEX_34E'.or.ring%ele( j )%name(1:4).eq.'Q34E'.or.ring%ele( j )%name(1:4).eq.'B35E')then ring%ele(j)%value(x1_limit$)=0.045 ring%ele(j)%value(x2_limit$)=0.045 ! type *,' Resetting aperture limit to 4.5 cm for element ',j,ring%ele(j)%name endif enddo enddo close(lun) endif if ( and(anal,4) .ne. 0 ) then ! test set_tunes routine write(6,3300) 3300 format(/' -------- Begin Analysis Type 3 ----------') species(1) = electron$ species(2) = positron$ ! print *,' Qtune with quads for one beam (1) or tune for both beams (2)' ! read *, tans tans=2 qh(2) = 0.54 qv(2) = 0.55 do qht = 0.52, 0.60, 0.01 do qvt = 0.52, 0.60, 0.01 qh(1) = qht qv(1) = qvt if( tans == 1 ) then call set_tunes(ring, pretzel, qtune, xqune, 1, & .true., .false., species, qh, qv, 0.001 ) else call set_tunes(ring, pretzel, qtune, xqune, 2, & .true., .true., species, qh, qv, 0.001 ) endif enddo enddo endif if( and(anal,8) .ne. 0 ) then ! Scan injection angle ! write(6,4400) !4400 format(/' -------- Begin Analysis Type 4 ----------') ! type *,'INIT_INJECTION always does a scan over injection angles so this analysis type 4 option is obsolete' ! 15 October 2003 JAC endif if( and(anal,16) .ne. 0 ) then ! Do injected beam tune scan write(6,5500) 5500 format(/' -------- Begin Analysis Type 5 ----------') call inj_tunescan( ring, inj, synch_beam, pretzel, & closed_orbit, & use_lrbbi, tunescan_file) ! ! Close tunescan ntuple file call tunescan_nt(1,0.,0.,0.,0,0.,0.) ! endif if( and(anal,32) .ne. 0 ) then ! Inject many particles, look at phase space write(6,6600) 6600 format(/' -------- Begin Analysis Type 6 ----------') call phase_space_track( ring, inj, synch_beam, pretzel, & use_lrbbi, phasespace_file ) endif if( and(anal,64) .ne. 0 ) then ! Test random phase space write(6,7700) 7700 format(/' -------- Begin Analysis Type 7 ----------') nsig = 3.0 totweight = 0.0 ! type *,' Enter number of particles' ! read *, n_part n_part=100 do i = 1,n_part ! call inject_particle with style=2 (for weighted uniform phase-sp) call inject_particle( ring, .true., 2, nsig, weight, inj, & synch_beam, inj_orbit ) totweight = totweight + weight write(10,400) weight, (inj_orbit(inj%ix_injpt)%vec(j), j=1,6) 400 format( 1x, 7(1x, e12.5) ) enddo type *,' Total weight = ',totweight endif if( and(anal,128) .ne. 0 ) then ! Efficiency of Injection Coordinates write(6,8800) 8800 format(/' -------- Begin Analysis Type 8 ----------') print *,'Find injection efficiency for which phase space coordinate?' print *,'1=x, 2=xprime, 3=y, 4=yprime' read *,coord call inj_coord_eff(ring, inj, synch_beam, closed_orbit, coord) endif if( and(anal,256) .ne. 0 ) then !Pinger Injection Efficiency write(6,9900) 9900 format(/' -------- Begin Analysis Type 9 ----------') call inj_pinger_eff(ring, inj, synch_beam) endif if( and(anal,512) .ne. 0) then !Scan for something that works! write(6,10100) 10100 format(/' -------- Begin Analysis Type 10 ----------') ! lun = lunget() ! open(lun, file='dat/scan_ele_orbits.dat', status = 'unknown') lun5 = lunget() datfname = 'dat/scan_output.dat' open(lun5, file=datfname, status = 'unknown') posmin = -0.10 posmax = -0.02 posstep = 0.004 angmin = -0.010 angmax = 0.0 angstep = 0.0004 try = 0 do injpos=posmin,posmax,posstep do injang=angmin,angmax,angstep print *, ' Injpos, injang=', injpos,injang try=try+1 if (mod(try,1000)==0) print *,'Running try: ',try call inject_particle( ring, .false., 0, 0., weight, inj, & synch_beam, inj_orbit ) inj_orbit(inj%ix_injpt)%vec(2)=injang inj_orbit(inj%ix_injpt)%vec(1)=injpos ! ring%ele(inj%ix_injpt)%value(x1_limit$)=0.10 ! ring%ele(inj%ix_injsex)%value(x1_limit$)=0.10 ! Be sure to include overlords between n_ele_use and n_ele_max, ! because they control the slaves which are created by the bbi elements do i=1,ring%n_ele_max if(ring%ele( i )%name(1:7).eq.'SEX_34E'.or.ring%ele( i )%name(1:4).eq.'Q34E'.or.ring%ele( i )%name(1:4).eq.'B35E')then ring%ele(i)%value(x1_limit$)=0.10 ring%ele(i)%value(x2_limit$)=0.10 ! type *,' Expanding aperture limit to 10 cm for element ',i,ring%ele(i)%name endif enddo ! Track from injection point (ix_start) to IP (ix=0) on turn 0 ! 16 April 2003 Since the last version of MANY_TRACK, the Z coordinate it expects to receive ! has changed sign. This means that xp,yp, and the time z need to have their sign changed ! before the first call. inj_orbit(ix_start)%vec(2) = - inj_orbit(ix_start)%vec(2) inj_orbit(ix_start)%vec(4) = - inj_orbit(ix_start)%vec(4) inj_orbit(ix_start)%vec(5) = - inj_orbit(ix_start)%vec(5) ! Take the beam to the IP ix_start = inj%ix_injpt norbwrt2 = inj%ix_injpt n=0 lostit=.false. do while ( n < 21 .and. .not. lostit ) call pulsed_bump(ring, n, inj ) ring%param%particle = antiparticle(ring%param%particle) call track_many (ring, inj_orbit, ix_start, 0, -1, track_state = track_state) ring%param%particle = antiparticle(ring%param%particle) ! if (lostit) print *,'Lostit!', try, n ! if (lostit==.false.) then ! write (lun5,*)'injpos: ',injpos,' injang: ',injang ! write (lun5,*)'On turn:',n,' At IP: ',inj_orbit(0)%vec(1), & ! inj_orbit(0)%vec(2) ! ! do i=0, norbwrt2 ! write(lun,501) try, n, ring%ele(i)%s, 1000.*inj_orbit(i)%vec(1), & ! 1000.*inj_orbit(i)%vec(2) ! 501 format(2x, i6,2x,i6, f10.4, 2x, f8.3,2x,f8.3) ! enddo ! endif ix_start = 0 norbwrt2 = ring%n_ele_track ! ring%ele(inj%ix_injpt)%value(x1_limit$)=0.045 ! ring%ele(inj%ix_injsex)%value(x1_limit$)=0.045 ! Be sure to include overlords between n_ele_use and n_ele_max, ! because they control the slaves which are created by the bbi elements do i=1,ring%n_ele_max if(ring%ele( i )%name(1:7).eq.'SEX_34E'.or.ring%ele( i )%name(1:4).eq.'Q34E'.or.ring%ele( i )%name(1:4).eq.'B35E')then ring%ele(i)%value(x1_limit$)=0.045 ring%ele(i)%value(x2_limit$)=0.045 ! type *,' Resetting aperture limit to 4.5 cm for element ',i,ring%ele(i)%name endif enddo n=n+1 enddo if (track_state == moving_forward$ ) then print *,'Found position and angle with successful injection:' print *,' Injpos: ',injpos,' injang: ',injang write (lun5,*)'Injpos - ',injpos,' injang - ',injang else print *,'Lost it. Injpos: ',injpos,' injang: ',injang endif enddo enddo ! close(lun) close(lun5) endif ! ============================================ ! End branching over chosen analysis type ! ============================================ ! 9999 continue ! deallocate (inj_orbit) ! deallocate (closed_orbit) ! call date_and_time(date, time) write(6,9100)date,time 9100 format(' ==> Exiting Program ENVELOPE on ',a8,' at ',a10) end program envelope !################################################################# subroutine phase_space_track( ring, inj, synch_beam, pretzel, & use_lrbbi, fname ) use injlib implicit none type (lat_struct) ring type (coord_struct), allocatable, save :: closed_orbit(:) type (coord_struct), allocatable, save :: co(:) type (coord_struct) inj_orbit(0:ring%n_ele_max) type (inj_struct) inj type (synch_beam_struct) synch_beam type (pretzel_struct) pretzel type (normal_modes_struct) mode type (coord_struct), allocatable :: diff_co(:) type (coord_struct), allocatable :: orbit(:) integer i_dim real(rdef), allocatable :: dk1(:) real(rdef) phi_a, phi_b integer int_Q_x integer int_Q_y logical ok integer i,j,species(2), ix_start, n integer lun real totweight, weight, injweight, qh(2), qv(2) real qht, qvt, eff, epsx, epsy, dele character fname*40 logical lostit integer track_state character datfname*40 ! Namelist variables integer nturns, n_part, nwrite, style real starting_tunes(3), nsig logical use_default_tunes, random, use_lrbbi logical use_coupling, aplimenable namelist / ps_track_input / random, & use_default_tunes, starting_tunes, use_coupling, & nturns, nsig, n_part, nwrite, style, aplimenable ! allocate storage call reallocate_coord (closed_orbit, ring%n_ele_max) call reallocate_coord(co,ring%n_ele_max) call reallocate_coord (orbit, ring%n_ele_max) call reallocate_coord(diff_co,ring%n_ele_max) ! call reallocate_coord (inj_orbit, ring%n_ele_max) ! Read namelist file ! print *,' Enter the phase space tracking input filename' ! read *, fname ! Default value for backward compatibility aplimenable = .false. print *,' Routine PHASE_SPACE_TRACK opening input file ',fname lun = lunget() open(lun,file=fname, status= 'old') read(lun, nml=ps_track_input) close(lun) ! Print phasespace tracking input parameters ! write(6,1000)fname, random, & use_default_tunes, starting_tunes, use_coupling, & nturns, nsig, n_part, nwrite, style, aplimenable 1000 format(' ==> Routine PHASE_SPACE_TRACK initialized with file ',a40/, & ' The RANDOM choice is: ',l/, & ' USE_DEFAULT_TUNES is:'l/, & ' The custom tune choices are:',3f7.3/, & ' The coupling setup choice is:',l/, & ' The number of turns is: ',i6/, & ' The injection phasespace sigma range is: ',f6.1/, & ' The number of tracked particles is: ',i6/, & ' The writing frequence NWRITE is: ',i5/, & ' The injection randomness style is: ',i3/, & ' Aperture limit enable is:',l) ! Open text file for turn-by-turn phase space at injection point print *,' Routine PHASE_SPACE_TRACK opening output file ',fname lun = lunget() write(datfname,1100) 1100 format('dat/ps_track.dat') open(lun, file=datfname, status = 'unknown') ! open(lun, file='dat/ps_track.dat', status = 'unknown') ! Choose whether to stop tracking if particle lost bmad_com%aperture_limit_on = aplimenable ! keep tracking no matter what ! species(1) = electron$ ! species(2) = positron$ ! Calculate closed orbit for electrons ! call twiss_and_track( ring, closed_orbit ) call twiss_at_start(ring) closed_orbit(0)%vec = 0. call closed_orbit_calc(ring, closed_orbit, 4 ) ok = .true. if( use_default_tunes ) then qht = ring%a%tune/twopi qvt = ring%b%tune/twopi else qht = starting_tunes(1) qvt = starting_tunes(2) ! Now set the tunes ! call set_tunes(ring, pretzel, qtune, xqune, 1, & ! .true., .false., species, qh, qv, 0.001 ) int_Q_x = int(ring%ele(ring%n_ele_track)%a%phi / twopi) int_Q_y = int(ring%ele(ring%n_ele_track)%b%phi / twopi) phi_a = ( int_Q_x + qht ) * twopi phi_b = ( int_Q_y + qvt ) * twopi ! allocate(dk1(ring%n_ele_max)) call re_allocate(dk1,ring%n_ele_max) call choose_quads(ring, dk1) call custom_set_tune (phi_a, phi_b, dk1, ring, closed_orbit, ok) deallocate(dk1) if ( .not. ok ) then print *,' Rtn PHASE_SPACE_TRACK: ERROR from CUSTOM_SET_TUNE, quitting.' goto 999 endif endif ! Now set the longitudinal tune ! Set up six-dimensional tracking including longitudinal coordinates ! Include damping effects, but exclude fluctuations bmad_com%radiation_damping_on = .true. bmad_com%radiation_fluctuations_on = .false. ! Turn on RF cavities call set_on_off (rfcavity$, ring, on$) ! print *,'ring%z%tune is',ring%z%tune if ( .not. use_default_tunes ) ring%z%tune = starting_tunes(3) * twopi call set_z_tune(ring) !------------------------------------------------------------------------ ! Update obsolete tracking 30 nov 2009 jac ! call twiss_at_start(ring) ! co(0)%vec=0. ! call closed_orbit_calc(ring, co, 6) ! call lat_make_mat6(ring, -1, co) ! call twiss_at_start(ring) ring%param%particle = electron$ i_dim = 4 call tracking_update (ring, orbit, i_dim, diff_co) !------------------------------------------------------------------------ write(6,1850)ring%a%tune/twopi,ring%b%tune/twopi,ring%z%tune/twopi 1850 format('Rtn PHASE_SPACE_TRACK: e- tunes after turning on RF: Qx = ',f6.3, ' Qy = ',f6.3, ' Qz = ',f6.3) ! turn on LRBBI and do emittance calc if( use_lrbbi ) then ! call emitt_calc( ring, all$, mode ) call radiation_integrals( ring, closed_orbit, mode ) epsx = mode%a%emittance dele = mode%sige_e if( use_coupling ) then epsx = epsx/2. epsy = epsx else epsy = 0.01*epsx endif ! call lrbb_strong_init (ring, 1, use_coupling, n_lr, ix_lr ) ! Turn on both IP and parasitic crossing BBI call lrbbi_ele_setup (ring, .true., .true. ) print *,' Rtn PHASE_SPACE_TRACK: LRBBI with ',ring%param%n_part,' particles/bunch' ! do i=1,n_lr ! call type_ele( ring.ele_(ix_lr(i)), .false., 0, .false., .false., ring) ! type *,' s = ',ring.ele_(ix_lr(i)).s,' e- = ', ! + closed_orbit(ix_lr(i)).vec(1),' e+ = ', ! + ring.ele_(ix_lr(i)).value(x_offset$) ! enddo endif nsig = abs(nsig) totweight = 0.0 injweight = 0.0 do i=1,n_part ! Tunes are where we want them, now inject ! style = 0 gives sampling of the 2-d phase space uniform in (x,xp), etc within a square of +-nsig. Gauss weight returned. ! style = 1 gives sampling of the 2-d phase space uniform in (x,xp), etc within a radius of nsig. Gauss weight returned. ! style = 2 gives sampling of the 2-d phase space uniform in r=sqrt(x**2+xp**2), etc. Gauss weight returned. ! style = 3 gives gaussian-distributed 2-d phase space. Limit to +-nsig. Returned weight = 1.0 call inject_particle( ring, random, style, nsig, weight, inj, & synch_beam, inj_orbit ) totweight = totweight + weight call track_injected_particle( ring, inj, electron$, inj_orbit, nturns, & .true., nwrite, lun, i, weight, lostit,track_state) if ( .not. lostit ) injweight = injweight + weight ! if( mod(i,10) == 0 ) then ! print *,' Finished tracking for particle ',i ! endif enddo 999 continue close( lun ) end subroutine phase_space_track !################################################################# subroutine inj_tunescan( ring, inj, synch_beam, pretzel, & closed_orbit, use_lrbbi, fname) use injlib implicit none type (lat_struct) ring type (coord_struct) closed_orbit(:) type (coord_struct), allocatable :: diff_co(:) type (coord_struct) delta_ip type (coord_struct), allocatable :: inj_orbit(:) type (coord_struct), allocatable, save :: co(:) type (inj_struct) inj type (synch_beam_struct) synch_beam type (pretzel_struct) pretzel type (normal_modes_struct) mode real(rdef), allocatable :: dk1(:) logical ok integer i,j, species(2), ix_start, n integer lun, lun2 real lost_array(0:ring%n_ele_track) real totweight, weight, injweight, qh(2), qv(2) real qht, qvt, eff, epsx, epsy, dele real(rdef) phi_a, phi_b character fname*40 logical lostit, err_flag integer i_dim, track_state integer ih,iv,nhstep,nvstep real qhstart,qvstart integer int_Q_x, int_Q_y integer nquad ! character date*10, time*10 character datfname*40 ! Namelist variables integer single_or_two_beam_scan, nturns, n_part, style real starting_ele_tunes(2), starting_pos_tunes(2) real qhmin, qhmax, qhstep, qvmin, qvmax, qvstep, nsig, dq_coupling real qz logical use_default_tunes, random, use_lrbbi logical use_coupling namelist / scan_input / random, & use_coupling, dq_coupling, & qz, & qhmin, qhmax, qhstep, qvmin, qvmax, qvstep, nturns, nsig, n_part, style ! allocate storage call reallocate_coord (inj_orbit, ring%n_ele_max) call reallocate_coord (co, ring%n_ele_max) call reallocate_coord(diff_co,ring%n_ele_max) ! do i=0,ring%n_ele_track lost_array(i) = 0. enddo ! Read namelist file lun = lunget() open(lun,file=fname, status= 'old') read(lun, nml=scan_input) close(lun) write(6,1000)fname, random, & use_coupling, dq_coupling, & qz, & qhmin, qhmax, qhstep, qvmin, qvmax, qvstep, nturns, nsig, n_part, style 1000 format(' ==> Program INJ_TUNESCAN initialized with file ',a40/, & ' RANDOM is: ',l/, & ' USE_COUPLING is: 'l/, & ' DQ_COUPLING is: ',f6.3/, & ' Q_z is: ',f6.3/, & ' QHMIN is: ',f6.3/, & ' QHMAX is: ',f6.3/, & ' QHSTEP is: ',f6.3/, & ' QVMIN is: ',f6.3/, & ' QVMAX is: ',f6.3/, & ' QVSTEP is: ',f6.3/, & ' NTURNS is: ',i6/, & ' NSIG is: ',f6.3/, & ' N_PART is: ',i6/, & ' STYLE is: ',i6) lun = lunget() datfname = 'dat/inj_tunescan.dat' open(lun, file=datfname, status = 'unknown') lun2 = lunget() datfname = 'dat/inj_lost_array.dat' open(lun2, file=datfname, status = 'unknown') if( use_lrbbi ) then ! We will calculate the tunes to be set and the injection ! efficiencies including the effects of both the LRBBI ! at the parasitic crossings and the BBI at the interaction ! point. Normally the latter should be horizontally defocussing ! and vertically focusing, just like the former, ! because PRZ13 should have been set to keep the beams separated ! horizontally. call lrbbi_ele_setup (ring, .true., .true. ) endif ! ! Print out initial tunes ring%param%particle = electron$ ! call twiss_at_start(ring) ! co(0)%vec = 0. ! call closed_orbit_calc(ring, co, 6) !! call closed_orbit_at_start(ring, co(0), 6, .true.) !! call track_all (ring, co) ! call lat_make_mat6(ring,-1,co) ! call twiss_at_start(ring) ! call twiss_propagate_all(ring) call tracking_update ( ring, co, 6 , diff_co) type * type *,' Rtn INJ_TUNESCAN: Initial e- tunes and closed orbit ' write(6,3000)ring%a%tune/twopi,ring%b%tune/twopi,ring%z%tune/twopi 3000 format(' Qx = ',f6.3,' Qy = ',f6.3,' Qz = ',f6.3) type '(a15,4e12.4)',' Closed orbit ', co(0)%vec(1:4) ! Set synchrotron tune ring%z%tune = qz * twopi call set_z_tune(ring) ! Use 6-dimensional tracking i_dim = 6 ! call twiss_at_start(ring) ! co(0)%vec = 0. ! call closed_orbit_calc(ring, co, 6) !! call closed_orbit_at_start(ring, co(0), i_dim, .true.) !! call track_all (ring, co) ! call lat_make_mat6(ring,-1,co) ! call twiss_at_start(ring) ! call twiss_propagate_all(ring) call tracking_update ( ring, co, i_dim , diff_co) type * type *,' Rtn INJ_TUNESCAN: e- tunes and closed orbit after setting Q_z' write(6,3000)ring%a%tune/twopi,ring%b%tune/twopi,ring%z%tune/twopi type '(a15,4e12.4)',' Closed orbit ', co(0)%vec(1:4) qhstart = ring%a%tune/twopi qvstart = ring%b%tune/twopi ! Enable aperture limits for tracking bmad_com%aperture_limit_on = .true. ! Turn on Coupling if( use_coupling ) then endif ! Set range over which injection phase space will be generated nsig = abs(nsig) ! Avoid exiting in closed_orbit_at_start global_com%exit_on_error = .false. ! If step size .le. zero then use only design tune if ( qhstep .gt. 0 )then nhstep = int((qhmax-qhmin+qhstep/2.)/qhstep)+1 else nhstep = 0 endif if ( qvstep .gt. 0 )then nvstep = int((qvmax-qvmin+qvstep/2.)/qvstep)+1 else nvstep = 0 endif ! do qht = qhmin, qhmax, qhstep ! do qvt = qvmin, qvmax, qvstep ! ! Always calculate efficiency for the design tune do ih = 0,nhstep do iv = 0,nvstep ! call date_and_time(date, time) write(6,2100)date,time 2100 format(' ==> Rtn INJ_TUNESCAN tune loop on ',a8,' at ',a10) if ( ih .eq. 0 .and. iv .eq. 0 ) then qht = qhstart qvt = qvstart else if ( ih .eq. 0 ) then qht = qhstart else qht = qhmin + ( ih - 1 ) * qhstep endif if ( iv .eq. 0 ) then qvt = qvstart else qvt = qvmin + ( iv - 1 ) * qvstep endif ! Qtune int_Q_x = int(ring%ele(ring%n_ele_track)%a%phi / twopi) int_Q_y = int(ring%ele(ring%n_ele_track)%b%phi / twopi) phi_a = ( int_Q_x + qht ) * twopi phi_b = ( int_Q_y + qvt ) * twopi ! allocate(dk1(ring%n_ele_max)) call re_allocate(dk1,ring%n_ele_max) call choose_quads(ring, dk1) nquad = 0 print *, 'Rtn CHOOSE_QUADS chose the following quads:' do i=1,ring%n_ele_track if ( dk1(i) .ne. 0 ) then nquad = nquad + 1 ! write ( 6, 5000 ) nquad, ring%ele( i )%name, dk1(i) 5000 format(i5,2x,a15,3x,'dk = ',f6.1) endif enddo call custom_set_tune (phi_a, phi_b, dk1, ring, co, ok) deallocate(dk1) if (.not. ok ) then print *,' Rtn INJ_TUNESCAN: ERROR from CUSTOM_SET_TUNE! qht,qvt=',qht,qvt exit endif endif ! Calculate beam-beam separations at all elements (delta_ip) ! and print out separation at IP ! Be sure particle type is electron (weak beam) when calling beambeam_separation ring%param%particle = electron$ call beambeam_separation ( ring, delta_ip, 6 ) if ( abs(delta_ip%vec(1) ) .lt. 0.5e-3 )then write ( 6, 4000) 1000*delta_ip%vec(1) 4000 format(' Rtn INJ_TUNESCAN Warning! Horizontal beam separation is only ', & f8.3,' mm !!') endif ! ! call twiss_at_start(ring) ! co(0)%vec = 0. ! call closed_orbit_calc(ring, co, i_dim ) !! call closed_orbit_at_start(ring, co(0), i_dim, .true.) !! call track_all (ring, co) ! call lat_make_mat6(ring,-1,co) ! call twiss_at_start(ring) ! call twiss_propagate_all(ring) call tracking_update ( ring, co, i_dim , diff_co, err_flag) ! type * ! type *,' Rtn INJ_TUNESCAN: e- tunes and closed orbit after setting Q_x, Q_y' ! write(6,3000)ring%a%tune/twopi,ring%b%tune/twopi,ring%z%tune/twopi ! type '(a15,4e12.4)',' Closed orbit ', co(0)%vec(1:4) ! Check to see if bmad status okay. If not, skip this tune value if ( err_flag ) then write ( 6, 6000 ) qht, qvt 6000 format ( ' Rtn INJ_TUNESCAN: Could not calculate tunes for set values QHT, QVT=',2f8.3,' Skipping to next set values.' ) exit endif totweight = 0.0 injweight = 0.0 do i=1,n_part ! Tunes are where we want them, now inject (with style=2 for weighted) call inject_particle( ring, random, style, nsig, weight, inj, & synch_beam, inj_orbit ) totweight = totweight + weight call track_injected_particle( ring, inj, electron$, inj_orbit, nturns, & .false., 20, 1, i, weight, lostit, track_state) if( lostit ) then lost_array(track_state) = lost_array(track_state) + weight else injweight = injweight + weight endif enddo eff = injweight/totweight ! Write out efficiency scan data write(lun,300) ring%a%tune/twopi,ring%b%tune/twopi, eff write(6,3100) ring%a%tune/twopi,ring%b%tune/twopi,ring%z%tune/twopi, & totweight, injweight, 100*eff 3100 format(' Rtn INJ_TUNESCAN: Qx= ',f6.3,' Qy= ',f6.3,' Qz= ',f6.3, & ' Totwt= ',e12.5,' Injwt= ',e12.5,' Efficiency: ',f6.2,'%') ! print *,' H tune, V tune, injection efficiency=',qht, qvt, eff ! print *,' injw = ',injweight,' totw = ',totweight 300 format(1x,f7.4,1x,f7.4,1x,f7.4) ! Write out the lost array here do i=0,ring%n_ele_track lost_array(i) = lost_array(i)/totweight write(lun2, 301) i, ring%ele(i)%s, lost_array(i) 301 format(1x,i5, 1x, f10.4, 1x, f10.6 ) ! Fill tune scan ntuple call tunescan_nt(0, qht, qvt, eff, i, real(ring%ele(i)%s), lost_array(i)) ! enddo enddo enddo end subroutine inj_tunescan !####################################### subroutine print_bumping( ring, n, inj ) use bmad use injparam_mod implicit none type (lat_struct) ring type (inj_struct) inj integer n, i real arg do i= 1,3 arg = inj%omega*( real(n) *ring%param%total_length/c_light + & inj%t0(i) ) write(6,1000)ring%ele(inj%ix_bumper(i))%name,10**9*inj%t0(i),180*arg/pi,1000*ring%ele(inj%ix_bumper(i))%value(hkick$) 1000 format(' *** Bumper ',a10,' T0=',f8.3,' ns', & ' Phase of half sine wave=',f6.1,' degs', & ' Kick=',f7.3,' mrad') enddo ! Now do the pinger write(6,2000)1000*ring%ele(inj%ix_pinger)%value(hkick$) 2000 format(' *** Pinger kick=',f7.3' mrad') end subroutine print_bumping !########################################################### subroutine tunescan_nt(iflag, qht, qvt, eff, iele, s, lostwt) ! ! Book and fill ntuple with tune scan injection efficiency data ! implicit none ! Input arguments real qht !horizontal tune value real qvt !vertical tune value real eff !injection efficiency integer iele !ordinal number of position in ring real s !longitudinal trajectory position in ring real lostwt !fraction of total inefficiency lost at this position integer iflag !zero if data entry, 1 if job end ! integer, save :: ntlun integer istat,icycle character ntname*40 ! integer i,k,id0,nvar character*12 htitle character*8 tags(6) real vec(6) ! data nvar/6/ data htitle/' Tune Scan '/ data tags/' qht ',' qvt ',' eff ','ielement', & ' s ',' lostwt '/ ! integer ienter, lunget data ienter/0/ ! !----------------------------------------------------------------------- ! if(iflag.eq.1)goto 900 ! ienter=ienter+1 if(ienter.eq.1)then ! Open file for ntuple ntlun=lunget() ntname='ntuples/tunescan.ntu' print *,' Opening ntuple file ',ntname call hrOPEN(ntlun,'NT',ntname,'N',1024,istat) ! Book ntuple id0=10 call hbookn(id0,htitle,nvar,'//NT',10000,tags) endif ! vec(1)=qht vec(2)=qvt vec(3)=eff vec(4)=iele vec(5)=s vec(6)=lostwt ! print *,vec call hfn(id0,vec) goto 950 ! 900 continue !================================== ! Write out histograms and ntuples print *,' Writing out ntuples to file ', ntname call hcdir('//NT',' ') CALL HROUT(0,ICYCLE,' ') CALL HREND('NT') close(ntlun) ! 950 continue ! end subroutine tunescan_nt !######################################################################## subroutine envelope_nt(ring, pretzel, inj, synch_beam, latname, coupling, n_trains_tot, n_cars, set_lrbbi, use_lrbbi, custom_pat) ! ! Book and fill ntuple with orbit envelope data ! use injlib use bunchcross_mod use constraints_mod ! implicit none type (lat_struct) ring type (pretzel_struct) pretzel type (inj_struct) inj type (coord_struct), allocatable, save :: inj_orbit(:), & elorbit(:), posorbit(:), & noprz_orbit(:),pb_orbit(:), diff_co(:) type (normal_modes_struct) mode type (synch_beam_struct) synch_beam type (pc_struct) pc type (constraint_struct) con ! !----------------------------- ! Input arguments character(*) latname logical coupling integer n_trains_tot, n_cars integer set_lrbbi logical use_lrbbi, custom_pat !----------------------------- character*50 name integer prefixi character*4 strjob ! character*80 infoname integer infolun real beaminfo(50) ! integer ntlun integer istat,icycle character ntname*80 ! integer n real pbx(0:ring%n_ele_track) real pby(0:ring%n_ele_track) real injorbx(0:ring%n_ele_track) real pbetax(0:ring%n_ele_track) real palphax(0:ring%n_ele_track) real petax(0:ring%n_ele_track) real ebetax(0:ring%n_ele_track) real ealphax(0:ring%n_ele_track) real eetax(0:ring%n_ele_track) real sign real xsave1,xsave2 ! integer ipcf(0:ring%n_ele_track) ! real dx,dxmin integer nearele(0:ring%n_ele_track) integer nbuntot,ip,ie,iside real xclear(0:ring%n_ele_track) real yclear(0:ring%n_ele_track) ! integer jj,nn ! real eps_ws,delta_ws,eps_wpi,delta_wi ! real eps_wi real epsx,betax,etax,xinj,x_lim real epsy,betay,etay,yinj,y_lim real xpretzel,ypretzel,xp34 ! integer i,ic,k integer id01,nvar1 character*12 htitle1 character*8 tags1(36) real vec(36) ! data nvar1/36/ data htitle1/' Inj Envelope '/ data tags1/' ind ',' key ',' s ', & 'pbetax ','palphax ','petax ', & 'ebetax ','ealphax ','eetax ', & ' xlim ', & 'storede ','storedp ','przdisp ','pbdisp ',' injorb ', & ' sigmax ', & 'xmaxinj ','xstrclrw','xinjclrw','xinjclrp', & ' sigmay ', & 'ymaxinj ','ystrclrw','yinjclrw','yinjclrp', & ' pcflag ', & 'xinjwbet','xinjweta', & 'strdex ','strdexp ', 'strdey ','strdeyp ', & 'strdpx ','strdpxp ', 'strdpy ','strdpyp '/ ! integer id02,nvar2 character*12 htitle2 character*8 tags2(14) ! data nvar2/14/ data htitle2/'Paracrossings '/ data tags2/' pbunch',' ebunch',' ipcross ',' side ', & ' nearel',' s ',' sdiff ', & ' xepsep ',' yepsep ', & ' xclear',' xsigp ',' xsige ',' yclear',' ysige '/ interface subroutine bunchcross(ring, con, pc) use bmad use bunchcross_mod use constraints_mod implicit none type (pc_struct) pc type (lat_struct) ring type (constraint_struct) con end subroutine end interface ! ! allocate storage call reallocate_coord (inj_orbit, ring%n_ele_max) call reallocate_coord (elorbit, ring%n_ele_max) call reallocate_coord (posorbit, ring%n_ele_max) call reallocate_coord (noprz_orbit, ring%n_ele_max) call reallocate_coord (pb_orbit, ring%n_ele_max) call reallocate_coord(diff_co, ring%n_ele_max) ! ! Define prefix for output files prefixi=index(latname,'.')-1 name=latname(1:prefixi) write (6,1000)name 1000 format(/'==> Entering routine ENVELOPE_NT with file prefix ',a30) ! !=================================== ! Open file for ntuple !=================================== ntlun=lunget() ntname='ntuples/job.ntu' write(6,3000)ntname 3000 format('Routine ENVELOPE_NT opening ntuple file ',a30) call hrOPEN(ntlun,'NT',ntname,'N',1024,istat) ! Book ntuples id01=10 call hbookn(id01,htitle1,nvar1,'//NT',10000,tags1) if( set_lrbbi.eq.1 .and. use_lrbbi )then id02=20 call hbookn(id02,htitle2,nvar2,'//NT',10000,tags2) endif ! ! Get pretzel-off orbit. Species choice should not matter. ! write(6,3333)ring%ele(inj%ix_injpt)%name,ring%ele(inj%ix_injpt)%a%beta 3333 format('Rtn ENVELOPE_NT before pretzel off. Beta-x(q34e):',a10,f10.2) !! print *,' Before pretzel off:',(i,ring%ele(i)%name,ring%ele(i)%a%beta,i=1,ring%n_ele_track) ! Turn Pretzel off print *,' Routine ENVELOPE_NT turning off pretzel' call set_pretzel(0, ring, pretzel) noprz_orbit(0)%vec(:) = 0. ! call closed_orbit_at_start (ring, noprz_orbit(0), 4, .true.) call closed_orbit_calc (ring, noprz_orbit, 4 ) ! We are tracking electrons backwards, so change sign of angles and time ! 19 Dec 2014 jac noprz_orbit(0)%vec(2) = -noprz_orbit(0)%vec(2) ! x-prime noprz_orbit(0)%vec(4) = -noprz_orbit(0)%vec(4) ! y-prime noprz_orbit(0)%vec(5) = -noprz_orbit(0)%vec(5) ! time ring%param%particle = antiparticle(ring%param%particle) call track_many (ring, noprz_orbit, 0, 0, -1) ring%param%particle = antiparticle(ring%param%particle) print *,' No-pretzel x position at Q34E:', noprz_orbit(inj%ix_injpt)%vec(1) write(6,4444)ring%ele(inj%ix_injpt)%name,ring%ele(inj%ix_injpt)%a%beta 4444 format('Rtn ENVELOPE_NT after pretzel off. Beta-x(q34e):',a10,f10.2) !! print *,' After pretzel off:',(i,ring%ele(i)%name,ring%ele(i)%a%beta,i=1,ring%n_ele_track) ! Turn pretzel on print *,' Routine ENVELOPE_NT turning pretzel on' ! Positrons ring%param%particle = positron$ ! Turn BBI elements off for positrons call lrbbi_ele_setup ( ring, .false., .false. ) ! SET_PRETZEL recalculates the transfer matrices and twiss functions by calling tracking_update (6/1/07) call set_pretzel(1, ring, pretzel) ! Store beta, alpha, eta do i=1,ring%n_ele_track pbetax(i)=ring%ele(i)%a%beta palphax(i)=ring%ele(i)%a%alpha petax(i)=ring%ele(i)%a%eta enddo call tracking_update (ring, posorbit, 4, diff_co) print *,' Pretzel e+ x position at Q34E:', posorbit(inj%ix_injpt)%vec(1) write(6,5555)ring%ele(inj%ix_injpt)%name,ring%ele(inj%ix_injpt)%a%beta 5555 format('Rtn ENVELOPE_NT after pretzel back on. e+ Beta-x(q34e):',a10,f10.2) !! print *,' After pretzel back on:',(i,ring%ele(i)%name,ring%ele(i)%a%beta,i=1,ring%n_ele_track) ! Electrons ring%param%particle = electron$ ! Turn BBI elements on for electrons call lrbbi_ele_setup ( ring, .true., .true. ) ! SET_PRETZEL recalculates the transfer matrices and twiss functions call set_pretzel(1, ring, pretzel) ! Store beta, alpha, eta do i=1,ring%n_ele_track ebetax(i)=ring%ele(i)%a%beta ealphax(i)=ring%ele(i)%a%alpha eetax(i)=ring%ele(i)%a%eta enddo call tracking_update (ring, elorbit, 4, diff_co) print *,' Pretzel e- x position at Q34E:', elorbit(inj%ix_injpt)%vec(1) write(6,6666)ring%ele(inj%ix_injpt)%name,ring%ele(inj%ix_injpt)%a%beta 6666 format('Rtn ENVELOPE_NT after pretzel back on. e- Beta-x(q34e):',a10,f10.2) ! Print out positions and angles of stored e+ and e- orbits at IP write(6,6700)(posorbit(1)%vec(i), i=1,4) 6700 format('Positron x,xprime,y,yprime at IP:',4e12.3) write(6,6710)(elorbit(1)%vec(i), i=1,4) 6710 format('Electron x,xprime,y,yprime at IP:',4e12.3) ! Now track stored beam with pretzel and bump pb_orbit(0)%vec(:) = 0. ! call closed_orbit_at_start (ring, pb_orbit(0), 4, .true.) call closed_orbit_calc (ring, pb_orbit, 4 ) ! We are tracking electrons backwards, so change sign of angles and time ! 19 Dec 2014 jac pb_orbit(0)%vec(2) = -pb_orbit(0)%vec(2) ! x-prime pb_orbit(0)%vec(4) = -pb_orbit(0)%vec(4) ! y-prime pb_orbit(0)%vec(5) = -pb_orbit(0)%vec(5) ! time do n=-5,50 call pulsed_bump(ring, n, inj ) ring%param%particle = antiparticle(ring%param%particle) pb_orbit(:)%species = antiparticle(ring%param%particle) call track_many (ring, pb_orbit, 0, 0, -1) ring%param%particle = antiparticle(ring%param%particle) if ( n .eq. 0 ) then do i=1,ring%n_ele_track pbx(i) = pb_orbit(i)%vec(1)- elorbit(i)%vec(1) ! 2 August 2005 Fixed bug (2)-->(3) ! pby(i) = pb_orbit(i)%vec(2)- elorbit(i)%vec(2) pby(i) = pb_orbit(i)%vec(3)- elorbit(i)%vec(3) enddo endif enddo ! Now track injected beam with pretzel and bump do i=1,ring%n_ele_track inj_orbit(i)%vec(:) = 0. enddo ! call closed_orbit_at_start (ring, inj_orbit(0), 4, .true.) call closed_orbit_calc (ring, inj_orbit, 4 ) ! Track injected electrons backward from injection point do i=1,6 sign=1. if (i.eq.2.or.i.eq.4.or.i.eq.5)sign=-1. inj_orbit(inj%ix_injpt)%vec(i) = sign*synch_beam%centroid(i) enddo do n=0,10 ! Raise aperture limits at injection point for first turn if(n == 0)then do i=1,ring%n_ele_track if(ring%ele( i )%name(1:7).eq.'SEX_34E'.or.ring%ele( i )%name(1:4).eq.'Q34E'.or.ring%ele( i )%name(1:4).eq.'B35E')then ring%ele(i)%value(x1_limit$)=0.10 ring%ele(i)%value(x2_limit$)=0.10 ! type *,' Expanding aperture limit to 10 cm for element ',i,ring%ele(i)%name endif enddo elseif(n==1)then do i=1,ring%n_ele_track if(ring%ele( i )%name(1:7).eq.'SEX_34E'.or.ring%ele( i )%name(1:4).eq.'Q34E'.or.ring%ele( i )%name(1:4).eq.'B35E')then ring%ele(i)%value(x1_limit$)=0.045 ring%ele(i)%value(x2_limit$)=0.045 ! type *,' Resetting aperture limit to 4.5 cm for element ',i,ring%ele(i)%name endif enddo ! ring%ele(inj%ix_injpt)%value(x1_limit$)=xsave1 ! ring%ele(inj%ix_injsex)%value(x1_limit$)=xsave2 endif call pulsed_bump(ring, n, inj ) ring%param%particle = antiparticle(ring%param%particle) call track_many (ring, inj_orbit, inj%ix_injpt, 0, -1) ring%param%particle = antiparticle(ring%param%particle) if ( n .eq. 0 ) then do i=1,ring%n_ele_track injorbx(i) = inj_orbit(i)%vec(1) enddo endif enddo ! ! Calculate squared normalized wake amplitudes here !--------------------------------------------------------------------- ! Squared normalized wake amplitude for the stored beam from the pulsed bump and pinger [m-rad] eps_ws = inj%xrms_ws**2 / ebetax(inj%ix_injpt) print *,' The effective emmittance contr. to the stored beam from pb and ping is ',eps_ws !--------------------------------------------------------------------- ! Squared normalized amplitude for the pinger's kick as given to the injected beam [m-rad] eps_wpi = inj%xrms_wpi**2 / ebetax(inj%ix_injpt) + ebetax(inj%ix_injpt)*inj%xprms_wpi**2 ! ! Calculate squared normalized amplitude of the injection oscillation xp34 = elorbit(inj%ix_injpt)%vec(1) - noprz_orbit(inj%ix_injpt)%vec(1) delta_wi = abs( synch_beam%centroid(6) ) eps_wi = abs( injorbx(inj%ix_injpt) ) - abs( xp34 ) & - abs( delta_wi * ring%ele(inj%ix_injpt)%a%eta ) & - sqrt( eps_wpi * ring%ele(inj%ix_injpt)%a%beta ) & - abs ( pbx(inj%ix_injpt) ) eps_wi = eps_wi**2 / ring%ele(inj%ix_injpt)%a%beta print *,' The effective emittance of the injection oscillation is',eps_wi ! ! Get parasitic interaction points ! 9x5 for full set ! 9x4 for 6 wiggler lattice ! 8x5 for chess lattice if ( set_lrbbi.eq.1 .and. use_lrbbi ) then con%n_trains = n_trains_tot con%n_cars = n_cars con%n_14ns_space = 1 nbuntot = con%n_trains*con%n_cars if (custom_pat) then con%BunchPattern='bunchcross.in' else if (n_trains_tot.eq.10)then con%BunchPattern='tenfour.sh_five' else con%BunchPattern='NULL' endif endif call bunchcross ( ring, con, pc ) ! print *,' Rtn BUNCHCROSS finds nr parasitic crossings: ',pc%total_pc ! ! Print output of Routine BUNCHCROSS write(6,4050)pc%total_pc 4050 format(/' --- Routine BUNCHCROSS for ',i3,' parasitic crossings ---') do jj=1,pc%total_pc write(6,4100) jj, pc%cross(jj)%ele%s 4100 format(' Crossing',i4,' at ',f5.1) write(6,4200) (pc%beam_beam(nn,jj), nn=1,nbuntot) 4200 format(54l2) enddo endif ! !======================================================================== ! Calculate ancillary information !======================================================================== beaminfo(1) = ring%ele(0)%value(E_TOT$)/1.e6 beaminfo(2) = ring%a%tune/twopi beaminfo(3) = ring%b%tune/twopi ! call emitt_calc( ring, all$, mode ) ring%param%particle = electron$ elorbit(:)%species = electron$ call radiation_integrals( ring, elorbit, mode ) ! Emittance for electrons beaminfo(4) = mode%a%emittance beaminfo(5) = mode%b%emittance beaminfo(6) = mode%sige_e beaminfo(7) = mode%sig_z ring%param%particle = positron$ posorbit(:)%species = positron$ call radiation_integrals( ring, posorbit, mode ) ! Emittance for positrons beaminfo(8) = mode%a%emittance beaminfo(9) = mode%b%emittance beaminfo(10) = mode%sige_e beaminfo(11) = mode%sig_z !!! For fully coupled case, share the hor and vert emmitances if ( coupling ) then beaminfo(4) = beaminfo(4)/2 beaminfo(5) = beaminfo(4) beaminfo(8) = beaminfo(8)/2 beaminfo(9) = beaminfo(8) endif ! Stored beam orbit position at pulsed bump maximum at Q34E beaminfo(12) = inj%xpb_stored(1) ! Stored beam orbit angle at pulsed bump maximum at Q34E beaminfo(13) = inj%xpb_stored(2) ! Injected beam position at Q34E beaminfo(14) = synch_beam%centroid(1) ! Injected beam angle at Q34E beaminfo(15) = synch_beam%centroid(2) ! Synch beam energy offset beaminfo(16) = synch_beam%centroid(6) ! Squared normalized wake amplitude for the stored beam from the pulsed bump and pinger [m-rad] beaminfo(17) = eps_ws ! Energy oscillation amplitude of the wake for the stored beam beaminfo(18) = inj%delta_ws ! Squared normalized amplitude for the pinger's kick as given to the injected beam [m-rad] beaminfo(19) = eps_wpi ! Energy oscillation of the injected beam beaminfo(20) = delta_wi ! Stored orbit without bump or pinger at injection point beaminfo(21) = elorbit(inj%ix_injpt)%vec(1) beaminfo(22) = elorbit(inj%ix_injpt)%vec(2) beaminfo(23) = elorbit(inj%ix_injpt)%vec(3) beaminfo(24) = elorbit(inj%ix_injpt)%vec(4) ! Electron orbit twiss at injection point beaminfo(25) = ebetax(inj%ix_injpt) beaminfo(26) = ealphax(inj%ix_injpt) beaminfo(27) = eetax(inj%ix_injpt) beaminfo(28) = inj%ix_injpt ! !=================================================== ! Load ntuple !=================================================== ! ! print *,' In envel:',(i,ring%ele(i)%name,ring%ele(i)%a%beta,i=1,ring%n_ele_track) ! ! Initialize p.c. flag preventing more than one ntuple entry per crossing do i=1,ring%n_ele_track ipcf(i) = 0 enddo do i=1,ring%n_ele_track vec(1) = i vec(2) = ring%ele(i)%key vec(3) = ring%ele(i)%s vec(4) = pbetax(i) vec(5) = palphax(i) vec(6) = petax(i) vec(7) = ebetax(i) vec(8) = ealphax(i) vec(9) = eetax(i) ! Horizontal aperture vec(10) = ring%ele(i)%value(x1_limit$) ! Stored electron orbit x vec(11) = elorbit(i)%vec(1) ! Stored positron orbit x vec(12) = posorbit(i)%vec(1) ! Pretzel displacement ! xpretzel = elorbit(i)%vec(1) - noprz_orbit(i)%vec(1) xpretzel = ( elorbit(i)%vec(1) - posorbit(i)%vec(1) ) / 2 ! Bug found Aug 1, 2005. vec(2)-->vec(3) ! ypretzel = ( elorbit(i)%vec(2) - posorbit(i)%vec(2) ) / 2 ypretzel = ( elorbit(i)%vec(3) - posorbit(i)%vec(3) ) / 2 vec(13) = xpretzel ! Pulsed bump displacement vec(14) = pbx(i) ! Injected orbit vec(15) = injorbx(i) ! Horizontal beam size for stored electrons epsx = beaminfo(4) betax = ring%ele(i)%a%beta etax = ring%ele(i)%a%eta if (epsx.gt.0.)then vec(16) = sqrt ( epsx*betax + (beaminfo(6) * etax)**2 ) else print *,' Warning in Rtn ENVELOPE_NT: EPSX <= 0' vec(16) = 0. endif ! Maximum horiz displacement during the pulsed bump or after the pinger xinj = max ( abs(pbx(i)) , sqrt ( eps_ws * betax ) ) vec(17) = xinj x_lim = ring%ele(i)%value(x1_limit$) ! Horiz Clearance of stored beam's centroid from the wall vec(18) = x_lim - abs( xpretzel ) - xinj - abs ( etax * inj%delta_ws ) ! HorizClearance of the injected beam's centroid from the wall vec(19) = x_lim - abs( xpretzel ) - sqrt( eps_wi * betax ) - abs( delta_wi * etax ) ! Horiz Clearance of the injected beam's centroid from the positron beam's centroid vec(20) = 2 * abs( xpretzel ) - sqrt( eps_wi * betax) - abs( delta_wi * etax ) & - sqrt( eps_ws * betax ) - abs( inj%delta_ws * etax ) xclear(i) = vec(20) ! Vertical beam size for injected electrons epsy = beaminfo(5) betay = ring%ele(i)%b%beta etay = ring%ele(i)%b%eta if(epsy.gt.0)then vec(21) = sqrt ( epsy*betay + (beaminfo(6) * etay)**2 ) else print *,' Warning in Rtn ENVELOPE_NT: EPSY <= 0' vec(21) = 0. endif ! Maximum vert displacement during the pulsed bump or after the pinger yinj = max ( abs(pby(i)) , sqrt ( eps_ws * betay ) ) vec(22) = yinj y_lim = ring%ele(i)%value(y1_limit$) ! Vertical Clearance of stored beam's centroid from the wall vec(23) = y_lim - abs( ypretzel ) - yinj - abs ( etay * inj%delta_ws ) ! Vertical Clearance of the injected beam's centroid from the wall ! 10 Aug 2004 Removed eps_wi contribution, because eps_wi is determined for horiz only ! vec(24) = y_lim - abs( ypretzel ) - sqrt( eps_wi * betay ) - abs( delta_wi * etay ) vec(24) = y_lim - abs( ypretzel ) - abs( delta_wi * etay ) ! Vertical Clearance of the injected beam's centroid from the positron beam's centroid ! 10 Aug 2004 Removed eps_wi and eps_ws contributions, because they were determined for horiz only ! vec(25) = 2 * abs( ypretzel ) - sqrt( eps_wi * betay) - abs( delta_wi * etay ) & ! - sqrt( eps_ws * betay ) - abs( inj%delta_ws * etay ) vec(25) = 2 * abs( ypretzel ) - abs( delta_wi * etay ) - abs( inj%delta_ws * etay ) yclear(i) = vec(25) ! vec(25) = 2 * abs( ypretzel ) ! ! Flag parasitic interaction points. Find close elements. vec(26) = 0. do ic = 1, pc%total_pc dx = abs ( ring%ele(i)%s - pc%cross(ic)%ele%s ) ! print *,' pcnr, pc_s, dx',ic,pc%cross(ic)%ele%s,dx if ( dx .lt. 1.5 ) then ! Load only if not already loaded if ( ipcf(ic).eq.0 ) then vec(26) = ic ipcf(ic) = 1 endif endif enddo ! ! Individual contributions to the injection electron beam's wall ! clearance from beta and eta vec(27) = sqrt( eps_wi * betax ) vec(28) = abs( delta_wi * etax ) ! ! Stored electron orbit x,xp,y,yp vec(29) = elorbit(i)%vec(1) vec(30) = elorbit(i)%vec(2) vec(31) = elorbit(i)%vec(3) vec(32) = elorbit(i)%vec(4) ! Stored positron orbit x,xp,y,yp vec(33) = posorbit(i)%vec(1) vec(34) = posorbit(i)%vec(2) vec(35) = posorbit(i)%vec(3) vec(36) = posorbit(i)%vec(4) !! call hfn(id01,vec) ! if ( ring%ele(i)%name .eq. 'Q34E' .or. & ring%ele(i)%name .eq. 'Q34W' ) write (6,5000) ring%ele(i)%name, vec 5000 format(' Routine ENVELOPE_NT makes ntuple entry for element ',a30,': ',5(10f10.3/)) ! enddo if( set_lrbbi.eq.1 .and. use_lrbbi )then! ! Find element closest to each parasitic crossing do ic = 1, pc%total_pc dxmin = 800. do i=1,ring%n_ele_track dx = abs ( ring%ele(i)%s - pc%cross(ic)%ele%s ) if ( dx .lt. dxmin ) then nearele(ic)=i dxmin=dx endif enddo enddo ! ! Now load parasitic crossing ntuple do ip=1,nbuntot do ie=1,1 ! electron only 1 bunch for now do iside=1,2 vec(1) = ip vec(2) = ie ic = pc%posbunch_elebunch(ip,ie,iside) vec(3) = ic vec(4) = iside vec(5) = nearele(ic) vec(6) = pc%cross(ic)%ele%s vec(7) = ring%ele(nearele(ic))%s - pc%cross(ic)%ele%s vec(8) = elorbit(nearele(ic))%vec(1) - posorbit(nearele(ic))%vec(1) vec(9) = elorbit(nearele(ic))%vec(3) - posorbit(nearele(ic))%vec(3) ! Clearance of injected electrons to positron beam centroid vec(10) = xclear ( nearele(ic) ) ! Horizontal beam size for injected electrons epsx = beaminfo(4) betax = ring%ele(nearele(ic))%a%beta etax = ring%ele(nearele(ic))%a%eta vec(11) = sqrt ( epsx*betax + (beaminfo(10) * etax)**2 ) vec(12) = sqrt ( epsx*betax + (beaminfo(6) * etax)**2 ) ! Clearance of injected electrons to positron beam centroid vec(13) = yclear ( nearele(ic) ) ! Vertical beam size for injected electrons epsy = beaminfo(5) betay = ring%ele(nearele(ic))%b%beta etay = ring%ele(nearele(ic))%b%eta vec(14) = sqrt ( epsy*betay + (beaminfo(6) * etay)**2 ) call hfn(id02,vec) enddo enddo enddo endif ! ! End of loading ntuples ! 900 continue !================================== ! Write out ntuples !================================== write( 6, 9000) ntname 9000 format(' Routine ENVELOPE_NT Writing out ntuples to file ',a20/ & '==> Exiting Routine ENVELOPE_NT'/) call hcdir('//NT',' ') CALL HROUT(0,ICYCLE,' ') CALL HREND('NT') close(ntlun) ! !=================================================== ! Write file with ancillary global information !=================================================== ! print *,'Stored beam orbit position at pulsed bump maximum at Q34E', inj%xpb_stored(1) print *,'Stored beam orbit angle at pulsed bump maximum at Q34E', inj%xpb_stored(2) print *,'Injected beam position at Q34E', synch_beam%centroid(1) print *,'Injected beam angle at Q34E', synch_beam%centroid(2) print *,'Synch beam energy offset', synch_beam%centroid(6) print *,'Displacement wake amplitude for the stored beam from the pulsed bump and pinger ', & inj%xrms_ws print *,'Squared normalized wake amplitude for the stored beam from the pulsed bump and pinger ', & eps_ws print *,'Energy oscillation amplitude of the wake for the stored beam', inj%delta_ws print *,'Displacement wake amplitude for the injected beam from the pulsed bump and pinger ', & inj%xrms_wpi print *,'Angular wake amplitude for the injected beam from the pulsed bump and pinger ', & inj%xprms_wpi print *,'Squared normalized amplitude for the pinger kick as given to the injected beam ', eps_wpi print *,'Energy oscillation of the injected beam ',delta_wi ! infolun = lunget() ! infoname = latname(1:prefixi)//'_'//strjob//'_info'//'.txt' ! infoname = 'job_'//strjob//'_info'//'.txt' infoname = 'job_info.txt' open(infolun, file=infoname, status = 'unknown') write(infolun,1100)beaminfo 1100 format(G22.7) close(infolun) ! end subroutine envelope_nt !######################################################################## ! Subroutine to eliminate differential displacement of strong and weak beam at IP ! by adjustment of PRETZING 13 ! ! Input: ! RING -- lat_struct: Ring ! FINAL_POS_IN -- Coord_struct: optional, position to which closed orbit should converge ! FINAL_POS_OUT -- Coord_struct: optional, actual position to which closed orbit has converged ! Output: ! RING -- lat_struct: Ring with separators adjusted to bring beams into collision ! J.A.Crittenden 18 April 2007 Adapted from BEAMBEAM_SCAN in the BSIM library subroutine collide (ring, i_dim) use injlib implicit none type ( lat_struct ) ring type(coord_struct) final_pos_in, final_pos_out type (coord_struct), allocatable, save :: co(:) type (coord_struct), allocatable, save :: diff_co(:) type (ele_struct) beambeam_ele integer i, i_dim integer ix_ip real(rdef) Qx, Qy, Qx_in , Qy_in real(rdef) rgamma, den, xi_v, xi_h real(rdef) r_0/2.81794092e-15/ ! print *,' Rtn COLLIDE called with i_dim=',i_dim call reallocate_coord(co,ring%n_ele_max) call reallocate_coord(diff_co,ring%n_ele_max) ! ! save initial tunes Qx = ring%a%tune/twopi Qx_in = Qx Qy = ring%b%tune/twopi Qy_in = Qy final_pos_in%vec(:) = 0. final_pos_out%vec(:) = 0 call close_pretzel (ring, i_dim, final_pos_in, final_pos_out) ! call twiss_at_start(ring) ! call closed_orbit_at_start(ring, co(0), i_dim, .true.) ! call track_all(ring, co) ! call lat_make_mat6(ring, -1, co) ! call twiss_at_start(ring) call tracking_update (ring, co, i_dim, diff_co) type*,'Rtn COLLIDE: after close pretzel but before close vertical: ' type '(1x,3(a9,f12.4))',' Qx = ',ring%a%tune/twopi,' Qy = ',ring%b%tune/twopi,' Qz = ',ring%z%tune/twopi call close_vertical(ring,i_dim,final_pos_in,final_pos_out) ! call twiss_at_start(ring) ! forall( i=0:ring%n_ele_track) co(i)%vec = 0. ! call closed_orbit_at_start(ring, co(0), i_dim, .true.) ! call track_all(ring, co) ! call lat_make_mat6(ring, -1, co) ! call twiss_at_start(ring) ! call twiss_propagate_all(ring) call tracking_update (ring, co, i_dim, diff_co) type * type *,'Rtn COLLIDE: After CLOSE VERTICAL ' type '(1x,3(a9,f12.4))',' Qx = ',ring%a%tune/twopi,' Qy = ',ring%b%tune/twopi,' Qz = ', ring%z%tune/twopi type '(a15,a8,a2,4e12.4)',' Closed orbit for ',species_name(ring%param%particle),'S:', co(0)%vec(1:4) write(6,'(a36,f7.4,a12,f7.4)')' Beam beam tune shifts: Delta Qx = ', ring%a%tune/twopi - Qx, & ' Delta Qy =',ring%b%tune/twopi-Qy type '(a47,f8.3,3x,f8.3)','Horizontal and Vertical Separations at IP (mm):', & 1000*diff_co(0)%vec(1),1000*diff_co(0)%vec(3) ! If vertical closing has ruined horizontal closure, call close_pretzel again if ( abs(diff_co(0)%vec(1)) > 1.e-5) then print *,' Rtn COLLIDE calling CLOSE_PRETZEL a second time' call close_pretzel (ring, i_dim, final_pos_in, final_pos_out) call tracking_update (ring, co, i_dim, diff_co) type '(a78,f8.3,3x,f8.3)', & ' Horizontal and Vertical Separations at IP (mm) after second and final attempt:', & 1000*diff_co(0)%vec(1),1000*diff_co(0)%vec(3) endif ! set rgamma, den, xi_v, and xi_h do i = 1, ring%n_ele_track if ( ring%ele(i)%name(1:5) .eq. 'IP_CO' ) then beambeam_ele = ring%ele(i) rgamma = ring%ele(0)%value(E_TOT$) den = twopi*rgamma*beambeam_ele%value(sig_x$)+beambeam_ele%value(sig_y$) xi_v = r_0*ring%param%n_part*ring%ele(0)%b%beta/beambeam_ele%value(sig_y$)/den xi_h = r_0*ring%param%n_part*ring%ele(0)%a%beta/beambeam_ele%value(sig_x$)/den write(6,'(2(a10,e12.4))')' xi_v = ', xi_v,' xi_h = ',xi_h endif enddo ! print *,' Exiting Rtn COLLIDE' ! end subroutine collide !######################################################################## ! Set bumps according to input file ! ! 22 May 2007 JAC subroutine set_bumps ( ring, i_dim ) use injlib implicit none ! type (lat_struct) ring type (coord_struct),allocatable :: orbit(:) integer i_dim type (coord_struct),allocatable :: diff_co(:) ! integer lun,lun1,istat character*40 fname integer nbmax data nbmax /200/ character*40 orbname character*40 bname(200) character*40 bn integer is integer iset(200) integer nbump integer i,ib,isin,found ! call reallocate_coord(orbit,ring%n_ele_max) call reallocate_coord(diff_co,ring%n_ele_max) ! ! Get present orbits call tracking_update (ring, orbit, i_dim, diff_co) ! ! Write initial orbit file lun1=lunget() write(orbname,1000) 1000 format('dat/bump_1.dat') open(lun1, file=orbname, status = 'unknown' ) print *,' Rtn SET_BUMPS writing file ',orbname do i=0, ring%n_ele_track write(lun1,1100)ring%ele(i)%s, & ring%ele(i)%name(1:10), & ring%ele(i)%key, & ring%ele(i)%a%phi, & ring%ele(i)%a%beta,ring%ele(i)%a%eta,& ring%ele(i)%b%phi, & ring%ele(i)%b%beta,ring%ele(i)%b%eta,& ring%ele(i)%mat6(2,1),ring%ele(i)%mat6(4,3), & diff_co(i)%vec(1),diff_co(i)%vec(2), & diff_co(i)%vec(3),diff_co(i)%vec(4), & diff_co(i)%vec(5),diff_co(i)%vec(6), & orbit(i)%vec(1),orbit(i)%vec(2), & orbit(i)%vec(3),orbit(i)%vec(4), & orbit(i)%vec(5),orbit(i)%vec(6), & ring%ele(i)%value(sig_x$),ring%ele(i)%value(sig_y$) 1100 format(f11.6,2x,a10,i6,22f15.10) enddo close(lun1) ! lun=lunget() fname='bump_setup.in' open(lun,file=fname,err=100,status='unknown',iostat=istat) goto 150 100 continue write(6,1200)istat,fname 1200 format('Rtn SET_BUMPS finds error',i5,' opening file ',a20, & ' ---- exiting') goto 990 150 continue ! nbump=0 200 continue read(lun,*,end=500,err=220)bn,is 2000 format(a20,i10) nbump=nbump+1 bname(nbump)=bn iset(nbump)=is if ( nbump == nbmax ) then print *,' Rtn SET_BUMPS reached limit. nbump=',nbump goto 500 endif goto 200 220 continue goto 200 500 continue close(lun) ! if ( nbump == 0) then print *,' Rtn SET_BUMPS found no bump entries in file ',fname, & ' -- exiting' goto 990 else print *,' Rtn SET_BUMPS found ',nbump,' bump entries in file ',fname endif ! ! do ib=1,nbump found=0 do i=1,ring%n_ele_max ! if (ring%ele(i)%name(1:15) == bname(ib)(1:15) ) then if (ring%ele(i)%name == bname(ib) ) then found=1 isin = ring%ele(i)%control%var(1)%value ring%ele(i)%control%var(1)%value = iset(ib) write(6,2900)bname(ib)(1:20),isin,iset(ib) 2900 format(' Rtn SET_BUMPS changes command for ',a20, & ' from ',i6,' to ',i6) ! print *,' Rtn SET_BUMPS changes command for ', bname(ib)(1:20), & ! ' from ',isin,' to ', iset(ib) endif enddo if ( found == 0 ) print *,' Rtn SET_BUMPS found no lattice element for ', & bname(ib) enddo ! ! Write orbit file with bumps ! call tracking_update (ring, orbit, i_dim, diff_co) ! lun1=lunget() write(orbname,3000) 3000 format('dat/bump_2.dat') open(lun1, file=orbname, status = 'unknown' ) print *,' Rtn SET_BUMPS writing file ',orbname do i=0, ring%n_ele_track write(lun1,3100)ring%ele(i)%s, & ring%ele(i)%name(1:10), & ring%ele(i)%key, & ring%ele(i)%a%phi, & ring%ele(i)%a%beta,ring%ele(i)%a%eta,& ring%ele(i)%b%phi, & ring%ele(i)%b%beta,ring%ele(i)%b%eta,& ring%ele(i)%mat6(2,1),ring%ele(i)%mat6(4,3), & diff_co(i)%vec(1),diff_co(i)%vec(2), & diff_co(i)%vec(3),diff_co(i)%vec(4), & diff_co(i)%vec(5),diff_co(i)%vec(6), & orbit(i)%vec(1),orbit(i)%vec(2), & orbit(i)%vec(3),orbit(i)%vec(4), & orbit(i)%vec(5),orbit(i)%vec(6), & ring%ele(i)%value(sig_x$),ring%ele(i)%value(sig_y$) 3100 format(f11.6,2x,a10,i6,22f15.10) enddo close(lun1) ! 990 continue ! end subroutine set_bumps !######################################################################## ! Turn off wigglers/undulators ! ! 23 December 2014 JAC subroutine turn_ccu_onoff ( ring, on) use bmad implicit none ! type (lat_struct) ring logical on integer i logical is_on1,is_on2 do i=0, ring%n_ele_track if( ring%ele(i)%key .eq. wiggler$ .and. ring%ele(i)%name(1:3).eq.'CCU' )then is_on1 = ring%ele(i)%is_on is_on2 = on ring%ele(i)%is_on = is_on2 write(6,'(a,x,i5,a,a,a,l,x,l)'),'Rtn TURN_CCU_ONOFF sets element ', & i,': ',trim(ring%ele(i)%name),' IS_ON before/after = ',is_on1,is_on2 endif enddo ! end subroutine turn_ccu_onoff !######################################################################## ! Weight beam-beam element strengths ! ! 24 December 2014 JAC subroutine weight_bbi ( ring, wt_lrbbi ) use bmad implicit none ! type (lat_struct) ring integer i real(rdef) wt_lrbbi real wt1, wt2 do i=0, ring%n_ele_track if( ring%ele(i)%key .eq. beambeam$ )then wt1 = ring%ele(i)%value(charge$) wt2 = wt_lrbbi * wt1 ring%ele(i)%value(charge$) = wt2 write(6,'(a,x,i5,a,a,f7.3,x,f7.3)'),'Rtn WEIGHT_BBI changes beam-beam element strength ', & i,ring%ele(i)%name,' CHARGE before/after = ',wt1,wt2 endif enddo ! end subroutine weight_bbi