program bunch_bbi ! ! Study tune characteristics of individual bunches ! ! J.A.Crittenden 8/16/04 ! ! Rename from BUNCH to BUNCH_BBI when entering into the ! CESR SVN repository jac 9/27/2012 !=================================================================== use bmad use bmadz_mod use bmadz_interface use bb_setup_mod use injparam_mod use injlib use reverse_mod ! use io_mod use jobidmod use bunchcross_mod use constraints_mod implicit none type (lat_struct) ring, ring_in, ring_out, firstring type (ele_struct), pointer :: ele type (coord_struct), allocatable :: closed_orbit(:) type (coord_struct), allocatable :: diff_co(:) type (coord_struct), allocatable :: co(:) type (coord_struct), allocatable :: first_co(:) type (coord_struct), allocatable :: orbit(:) type (coord_struct) orb0 real(rp), allocatable :: dk1(:) integer, parameter :: ncopts = 100 real(rp) firstorb(ncopts,6), orb(ncopts,6), difforb(ncopts,6) type (ele_struct) dumele type (coord_struct) here real(rp) s,step integer nstep 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 type (normal_modes_struct) mode real(rp) sige_e real(rp) sig_z real(rp) e_loss real(rp) xemm real(rp) yemm real(rp) zemm real(rp) xalphad real(rp) yalphad real(rp) zalphad real(rp) xtau real(rp) ytau real(rp) ztau real(rp) alphaxstar real(rp) alphaystar real(rp) betaxstar real(rp) betaystar real(rp) etaxstar real(rp) etaystar logical off, error integer i, ii, j, jj, n real(rp) totpol integer nrampsteps integer imax integer lun,lundat,lundat2,lundat3,lundat4 integer set_lrbbi integer itrain, icar integer istat real(rp) xbetamax integer ixbetamax real(rp) ybetamax integer iybetamax real(rp) sxmax,symax character this_lat*100, latname*40, qtname*80 character lat_file*80 character fname*40 character lrbbi_file*40 real(rp) qx,qy,qz logical constant_tunes /.false./ real(rp) prz1, prz13 logical calc_dynap /.false./ character*40 ringmod_file/' '/ integer iq,iel character*4 qn(3) /'Q06W','Q08W','Q06E'/ real rdqcur,delcur,setk1,setqx,setqy real(rp) fact character datfname*40 integer newid logical ok,okbeta real qht, qvt real(rp) phi_a, phi_b integer int_Q_x, int_Q_y integer nquad real(rp) dxip,dyip logical fluctuations_on logical damping_on integer i_dim character date*10, time*10 integer version real(rp) current, lastcurrent, bunchwt integer strong_particle, n_trains_tot, n_cars, slices logical rec_taylor, beambeam_ip, ip real(rp) bbiwt(9,10) /90*0./ real(rp) maxcurrent, stepcurr, coupling_sb integer ic,nsteps integer itrainmax, icarmax integer ix_ip integer sdatlun character*40 string character*80 sdatname integer txlun character*80 txname logical lastentry real(rp) poscurrent real(rp) bparam,bcontribution real(rp) xsig2,xdel2 type (pc_struct) pc type (constraint_struct) con logical custom_pat ! Sync light source positions real(rp) ssource(2) /170.783,597.738/ real(rp) svec(14) !======================================= ! PAW storage real hm(30000000) common/pawc/hm !======================================= interface subroutine set_custom_pretzel (ring, co, off, prz1, prz13 ) use bmad implicit none type (lat_struct) ring type (coord_struct), allocatable :: co(:) logical off real(rp) prz1,prz13 end subroutine set_custom_pretzel end interface interface subroutine pcbetcor(ring, co, fact) use bmad implicit none type (lat_struct) ring type (coord_struct), allocatable :: co(:) real(rp) fact end subroutine pcbetcor end interface namelist / lrbbi_input / strong_particle, n_trains_tot, n_cars, bbiwt, & maxcurrent, stepcurr, rec_taylor, & beambeam_ip, coupling_sb, slices,custom_pat namelist / bunch_setup / this_lat, lrbbi_file, & qx, qy, qz, constant_tunes, & prz1, prz13, & calc_dynap, & ringmod_file ! procedure(track1_custom_def) :: track1_custom track_custom_ptr => track1_custom call date_and_time(date, time) write(6,11)date,time 11 format(' ==> Initializing Program BUNCH_BBI on ',a8,' at ',a10) !======================== ! Initialize HBOOK storage ! call hlimit(30000000) ! print *,' hm(1:10)=',hm(1:10) ! ! Get job ID number call readjobid ! Update job ID file now in go script ! newid=jobid+1 ! call writejobid(newid) !====================================================================== ! Program setup parameters !====================================================================== ! Read in setup info lun = lunget() fname='bunch_bbi_setup.in' open(lun,file=fname, status= 'old') read(lun, nml=bunch_setup) close(lun) ! write(6,1000)fname, this_lat, lrbbi_file, prz1, prz13, & qx, qy, qz, constant_tunes, & calc_dynap, ringmod_file 1000 format(' ==> Program BUNCH initialized with file ',a40/, & ' The chosen lattice is: ',a100/, & ' The lrbbi setup file name is: ',a40/, & ' The chosen PRZ1 value is: ',f7.1/, & ' The chosen PRZ13 value is: ',f7.1/, & ' The chosen tune values are: ',3f7.3/, & ' CONSTANT_TUNES is set to: ',l1/, & ' Dynamic aperture calculation: ',l/, & ' Ring modification file: ',a40) ! Read in long-range beam-beam interaction simulation setup info lun = lunget() open(lun,file=lrbbi_file, status= 'old') read(lun, nml=lrbbi_input) close(lun) ! 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 write(6,1050)lrbbi_file, strong_particle, n_trains_tot, n_cars, & bbiwt, maxcurrent, stepcurr, rec_taylor, & beambeam_ip, coupling_sb, slices 1050 format(' ==> LRBBI initialized with file ',a40/, & ' Strong beam particle type (-1:e- +1:e+): ',i4/, & ' Total nr of trains:',i4/, & ' Nr of cars/train: ',i4/, & ' BBI weights: ',9f10.2/,9(18x,9f10.2/), & ' Max e+ Current/bunch: ',f4.1/, & ' e+ current/bunch step size: ',f5.2/, & ' Taylor tracking: ',l/, & ' Include BBI at IP: ',l/, & ' Coupling: ',f7.4/, & ' Nr of slices: ',i4) ! Open output tune data file lundat=lunget() write(datfname,4500)jobid 4500 format('dat/bunch_',i4.4,'.dat') open(lundat, file=datfname, status = 'unknown') ! Open output orbit data file lundat2=lunget() write(datfname,4600)jobid 4600 format('dat/bunch_orbit_',i4.4,'.dat') open(lundat2, file=datfname, status = 'unknown') ! Open output synch light source data file lundat3=lunget() write(datfname,4700)jobid 4700 format('dat/synclight_',i4.4,'.dat') open(lundat3, file=datfname, status = 'unknown') ! Open output IP element data file lundat4=lunget() write(datfname,4800)jobid 4800 format('dat/ip_',i4.4,'.dat') open(lundat4, file=datfname, status = 'unknown') ! Initialize lattice. qtname=this_lat call bmad_parser(qtname, firstring) ! Allocate storage allocate (closed_orbit(0:firstring%n_ele_max)) call reallocate_coord(first_co,firstring%n_ele_max) call reallocate_coord(diff_co,firstring%n_ele_max) ! Read in aperture limit definitions ! Reduce Q28E aperture from 45 to 37 mm qtname='/nfs/acc/user/critten/cesr38/cesr/layout/cesrta_aperture_28e.bmad' print *,' Aperture file is:',qtname call bmad_parser2(qtname, firstring) call calc_z_tune (firstring%branch(0)) ! Set pretzel with positron orbit firstring%param%particle = positron$ call tracking_update(firstring, first_co, 4, diff_co) print *,' e+/e- separation x,angle at IP before setting pretzel:', & diff_co(0)%vec(1),diff_co(0)%vec(2) print *,' Positron x,xp,y,yp,x,dele at IP before setting pretzel:', & first_co(0)%vec(1:4) off = .false. call set_custom_pretzel(firstring, first_co, off, prz1, prz13) call tracking_update(firstring, first_co, 4, diff_co) print *,' e+/e- differences x,xp,y,yp,z,dele at IP after setting pretzel:', & diff_co(0)%vec(1:4) print *,' Positron x,xp,y,yp,x,dele at IP after setting pretzel:', & first_co(0)%vec(1:4) !!!!!!!!!!!!!! Begin loop over test bunches here print *, 'Before starting car/train loop, n_ele_max=',firstring%n_ele_max if (custom_pat) then itrainmax=1 icarmax=n_cars ! total number of weak beam bunches con%BunchPattern='bunchcross.in' call bunchcross (firstring, con, pc) else itrainmax = n_trains_tot icarmax = n_cars end if do itrain = 1,itrainmax do icar = 1,icarmax ! bunch number of the weak beam (usually electron) print *,' ====> Setting up ring for train,bunch',itrain,icar ! Re-initialize lattice ring = firstring call reallocate_coord(co,ring%n_ele_max) call reallocate_coord(diff_co,ring%n_ele_max) call reallocate_coord(orbit,ring%n_ele_max) ! Re-initialize closed orbit co = first_co ! For purposes of setting up the LRBBI and IP BBI, turn off ! radiation fluctuations and damping and turn off the RF bmad_com%radiation_fluctuations_on = .false. bmad_com%radiation_damping_on = .false. call set_on_off (rfcavity$, ring, off$) ! 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. ! Dummy nonzero current current = 0.000001 ring_in = ring if (custom_pat) then call lrbbi_setup_custom(ring_in, ring_out, strong_particle, con, pc, icar, current, rec_taylor) else if ( all (bbiwt .eq. 1.) ) then ! Do not use optional argument BBIWT for uniformly populated bunches call lrbbi_setup ( ring_in, ring_out, strong_particle, itrain, icar, n_trains_tot, n_cars, current, rec_taylor) else call lrbbi_setup ( ring_in, ring_out, strong_particle, itrain, icar, n_trains_tot, n_cars, current, rec_taylor, bbiwt) endif end if ring = ring_out bmad_com%aperture_limit_on = .false. ! Reallocate co because ring size changed call reallocate_coord(co,ring%n_ele_max) call reallocate_coord(diff_co,ring%n_ele_max) ! Restore electron tracking after lrbbi_setup 12 Dec 2005 ring%param%particle = electron$ call tracking_update(ring, co, 4, diff_co) print * print '(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) print '(a15,6e12.4)',' Closed orbit ', co(0)%vec(1:4) ! Initialize beam-beam interaction at interaction point and print out tunes if(beambeam_ip)then ! Set up the ring for positrons as expected by BB_SETUP ring%param%particle = positron$ ring%param%n_part = 0. call twiss_at_start(ring) co(0)%vec = 0. call closed_orbit_calc(ring, co, 4) call lat_make_mat6(ring,-1, co) call twiss_at_start(ring) call twiss_propagate_all(ring) ! Define BUNCHWT to weight the current of the colliding positron bunch. ! Used, for example, to turn off the IP BBI for e- train 9 all bunches ! when positrons in 8x5 configuration. ! Since BUNCHWT is used to multiply VALUE(CHARGE$) of the IP BBI element, ! which is the only element-specific parameter used to weight the strength ! of the interaction, it permits setting a strength which differs ! from the strengths of the parasitic BBI. 25 Aug 2004 jac ! ! Make robust for 10x4 ! bunchwt = bbiwt(itrain,icar) bunchwt = 1. call bb_setup(ring, strong_particle, current, bunchwt, coupling_sb, slices) ! Reallocate co because ring size changed 1 nov 04 call reallocate_coord(co,ring%n_ele_max) call tracking_update(ring, co, 4, diff_co) print '(a51,6e12.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) endif ! Finished adding IP BBI ! Write file of s positions for this configuration sdatlun = lunget() write(string,2210)n_trains_tot,n_cars,itrain,icar 2210 format(i2.2,'x',i2.2,'_',i1.1,'_',i2.2,'.txt') sdatname = 'pcross/bunch'//string open(sdatlun, file=sdatname, status = 'unknown',err=2500, iostat=istat) do i=1,ring%n_ele_track if ( ring%ele(i)%name(1:5) .eq. 'LRBBI' ) write(sdatlun,'(f11.6)')ring%ele(i)%s enddo close(sdatlun) goto 2600 2500 continue write(6,2550)sdatname,istat 2550 format(' Rtn BUNCH: Error opening ',a30/ & ' iostat=',i4,' Continuing ...') 2600 continue !================================================================= ! Reverse the ring !! 22 Sep 2004 ! Backward tracking does not work with damping, ! because closed_orbit_calc returns the closed orbit ! only for the forward-tracking ring. ! DYNAMIC_APERTURE does forward tracking, so we need ! to have it do so through the reversed ring for electrons. ! print *,' Now reversing ring ...' write(6,2221)(i,ring%ele(i)%name,ring%ele(i)%s,ring%ele(i)%a%beta,ring%ele(i)%b%beta,i=0,ring%n_ele_max) 2221 format(' Before reverse'/,4(i4.4,1x,a14,f8.2,2f7.1,2x)) ring_in = ring call lat_reverse ( ring_in, ring_out ) ring = ring_out ! Ensure beams are colliding ! 8 Sep 2009 Not for CHESS jac ! call collide ( ring, 4 ) ! write(6,2222)(i,ring%ele(i)%name,ring%ele(i)%s,ring%ele(i)%a%beta,ring%ele(i)%b%beta,i=0,ring%n_ele_max) 2222 format(' After reverse and collide'/,4(i4.4,1x,a14,f8.2,2f7.1,2x)) print *,' Before closed_orbit_calc',orbit(0)%vec(1:6) call closed_orbit_calc (ring, orbit, 4 ) print '(a50,6e12.4)',' Electron closed orbit at IP after reversing ring ', orbit(0)%vec(1:6) ! Reallocate co because ring size changed 1 nov 04 call reallocate_coord(co,ring%n_ele_max) call reallocate_coord(diff_co,ring%n_ele_max) call tracking_update(ring, co, 4, diff_co) print '(a50,6e12.4)',' Electron closed orbit at IP after reversing ring ', co(0)%vec(1:4) !================================================================= ! Check pretzel in reversed ring print *,' e+/e- differences x,xp,y,yp,z,dele at IP after reversing ring:', diff_co(0)%vec(1:4) !************************************************************ ! Set desired tunes with BBI current at zero ring%param%n_part = 0. ! 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 if ( qx .ne. 0. .or. qy .ne. 0. ) then allocate(dk1(ring%n_ele_max),stat=j) call choose_quads( ring, dk1 ) nquad = 0 do i=1,ring%n_ele_track if ( dk1(i) .ne. 0 ) then nquad = nquad + 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 print * print *,' Electron tunes after reversing ring:' print *,' Qx = ',ring%a%tune/twopi,' Qy = ',ring%b%tune/twopi,' Qz = ',ring%z%tune/twopi print '(a15,6e12.4)',' Closed orbit ', co(0)%vec(1:6) ! End of tune setting !************************************************************ ! Set up six-dimensional tracking including longitudinal coordinates ! Include damping effects, but exclude fluctuations bmad_com%radiation_fluctuations_on = .false. bmad_com%radiation_damping_on = .true. ! Turn on RF cavities call set_on_off (rfcavity$, ring, on$, co) if ( qz .ne. 0. ) then ! print *,'Before setting, ring%z%tune is',ring%z%tune ! print *,'Before calling set_z_tune, ring%z%tune is',ring%z%tune call set_z_tune(ring%branch(0), qz * twopi) i_dim = 6 else i_dim = 4 endif ! print *,'ring%z%tune is',ring%z%tune ring%param%particle = electron$ call tracking_update(ring, co, i_dim, diff_co, error) if(error) then write (6,2800) itrain,icar 2800 format(' BMAD tracking error found for train',i2, & ' car',i2,' current: 0 -- skip this bunch.') goto 8500 endif print * print *,' Electron tunes after setting Qz and turning on RF:' print *,' Qx = ',ring%a%tune/twopi,' Qy = ',ring%b%tune/twopi,' Qz = ',ring%z%tune/twopi print '(a50,6e12.4)',' Electron closed orbit at IP after turning on RF ', co(0)%vec(1:6) ! Print out ring elements and orbit x, y, dele ! print *,' Ring elements and closed orbit: index, s (m), x,y (mm) and dele (per mil):' ! write(6,2850)(i,ring%ele(i)%name,ring%ele(i)%s,1000*co(i)%vec(1:2),1000*co(i)%vec(6),i=0,ring%n_ele_track) 2850 format(4(i4.4,1x,a6,f8.2,2f7.1,f7.3,2x)) ! print *,' Ring element apertures: index, s (m), x1, x2, y1, y2 (cm)' ! write(6,2860)(i,ring%ele(i)%name,ring%ele(i)%s, & ! 100*ring%ele(i)%value(x1_limit$),100*ring%ele(i)%value(x2_limit$), & ! 100*ring%ele(i)%value(y1_limit$),100*ring%ele(i)%value(y2_limit$), & ! i=0,ring%n_ele_track) 2860 format(4(i4.4,1x,a6,f7.2,4f4.1,2x)) print * print *,' Electron tunes prior to HSEP test:' print *,' Qx = ',ring%a%tune/twopi,' Qy = ',ring%b%tune/twopi,' Qz = ',ring%z%tune/twopi print '(a50,6e12.4)',' Electron closed orbit at IP after turning on RF ', co(0)%vec(1:6) !--------------------------------------------------------------------------------------------- print *,' Turning off H seps as test !!!' off = .true. call set_custom_pretzel(ring, co, off, prz1, prz13) call tracking_update(ring, co, 6, diff_co) print *,' e+/e- differences x,xp,y,yp,z,dele at IP after turning off horiz seps:', & diff_co(0)%vec(1:6) print '(a24,6e12.4)',' Seps off closed orbit ', co(0)%vec(1:6) ! Turn seps back on off = .false. call set_custom_pretzel(ring, co, off, prz1, prz13) call tracking_update(ring, co, 6, diff_co) print *,' e+/e- differences x,xp,y,yp,z,dele at IP after turning horiz seps back on:', & diff_co(0)%vec(1:6) print '(a41,6e12.4)',' Closed orbit after turning seps back on ', co(0)%vec(1:6) !--------------------------------------------------------------------------------------------- print * print *,' Electron tunes after HSEP test:' print *,' Qx = ',ring%a%tune/twopi,' Qy = ',ring%b%tune/twopi,' Qz = ',ring%z%tune/twopi print '(a50,6e12.4)',' Electron closed orbit at IP after turning on RF ', co(0)%vec(1:6) ! Calculate global ring parameters call radiation_integrals( ring, co, mode ) ! Calculate dynamic aperture for this lattice and zero positron current if ( calc_dynap ) then current = 0. lastentry = .false. orb0 = co(0) ! Removed 23 jan 06 ! call dynap ( ring, orb0, jobid, lastentry, itrain, icar, current, mode ) endif ! Do not stop program if error in tracking_update global_com%exit_on_error = .false. ! Loop over strong-beam current/bunch !do current = stepcurr, maxcurrent, stepcurr !do current = 0, maxcurrent, stepcurr nsteps=int(maxcurrent/stepcurr+0.5) do ic=0,nsteps current = ic*stepcurr print *,' Current step ',ic,' Current: ',current ! Calculate total positron current if (custom_pat) then poscurrent = current * n_trains_tot else poscurrent = 0 do i=1,9 do j=1,5 poscurrent = poscurrent + current * bbiwt(i,j) enddo enddo endif print *,' Total positron current is ',poscurrent,' mA' ring%param%n_part = current*0.001 *(ring%param%total_length/c_light)/e_charge print *,' Number of particles per full bunch is ',ring%param%n_part ring%param%particle = electron$ call tracking_update(ring, co, i_dim, diff_co) ! Ensure beams are colliding ! 8 Sep 2009 Not for CHESS jac ! call collide ( ring, i_dim ) ! call collide ( ring, 4 ) !************************************************************ ! Set desired tunes ! Turn off IP do i = 1, ring%n_ele_track if ( ring%ele(i)%name(1:5) .eq. 'IP_CO' ) then ix_ip=i ! ring%ele(ix_ip)%is_on=.false. endif enddo if ( constant_tunes ) then allocate(dk1(ring%n_ele_max),stat=j) call choose_quads( ring, dk1 ) nquad = 0 do i=1,ring%n_ele_track if ( dk1(i) .ne. 0 ) then nquad = nquad + 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 write (6,2900) qht,qvt,current 2900 format(' Error from CUSTOM_SET_TUNE for qht,qvt =',2f6.2,& ' current=',f5.2,'-- skip this point.') goto 8000 endif print * print *,' Electron tunes after setting tunes:' print *,' Qx = ',ring%a%tune/twopi,' Qy = ',ring%b%tune/twopi,' Qz = ',ring%z%tune/twopi print '(a15,6e12.4)',' Closed orbit ', co(0)%vec(1:6) endif ! End of tune setting ! ! Turn IP back on if (beambeam_ip) ring%ele(ix_ip)%is_on=.true. !************************************************************ !***************************************************************** ! Modify quad K1 values to compensate for parasitic BBI ! Correction file is calculated for a positron current of 1 ma/bunch ! 23 November JAC/MGB if ( ringmod_file .ne. ' ') then fact = current / 1. ! call ringmod ( ring, ringmod_file, fact ) endif !***************************************************************** ! Implement PCBETING correction if available fact = current / 1. ! Assume optics optimized at 2.5 mA/bunch (fit_postop_20060116_qt_beta_seps) ! fact = ( current - 2.5 ) / 1. ! No correction at 2.0 ma as per mgb request ! fact = ( current - 2.0 ) / 1. ! call PCBETCOR( ring, co, fact ) !***************************************************************** ! Update the closed orbit with the new quad settings call tracking_update(ring, co, i_dim, diff_co, error) ! Skip quad current study go to 310 ! ! Check effect of changing quad currents at Q06W, Q08W and Q06E on tunes print *,' Quad current study:' print *,' -------------------' print *,' Train, bunch, e+ current, quad, quad curr factor, tunes, design beta x and y, calculated beta x and y=' do iq=1,3 do iel=1,ring%n_ele_max if (ring%ele(iel)%name .eq. qn(iq)) then ! Calculate beta values relative to original k1 settings setk1 = ring%ele(iel)%value(k1$) setqx = ring%a%tune/twopi setqy = ring%b%tune/twopi delcur = .001 do rdqcur=1.0,1.015,delcur ring%ele(iel)%value(k1$) = rdqcur*setk1 call tracking_update(ring, co, i_dim, diff_co, error) write(6,2955)itrain,icar,current, & qn(iq),rdqcur, & 390*ring%a%tune/twopi,390*ring%b%tune/twopi, & ring%ele(iel)%a%beta, & ring%ele(iel)%b%beta, & ! abs((2*twopi/(0.6*setk1))*(ring%a%tune/twopi-setqx)/(delcur)), & ! abs((2*twopi/(0.6*setk1))*(ring%b%tune/twopi-setqy)/(delcur)) (2*twopi/(0.6*setk1))*(ring%a%tune/twopi-setqx)/(rdqcur-1.), & (2*twopi/(0.6*setk1))*(ring%b%tune/twopi-setqy)/(rdqcur-1.) 2955 format('======',2i4,f5.2,3x,a4,f8.4,2(f10.3,' kHz'),4f8.2) enddo ring%ele(iel)%value(k1$) = setk1 call tracking_update(ring, co, i_dim, diff_co, error) endif enddo enddo 310 continue !***************************************************************** ! Write text files of ring elements: s, phase, beta, positron orbit !----------------------------------------------------------------- if (itrain.eq.1.and. icar.eq.1 .and. (ic.eq.0 .or. ic.eq.nsteps ) ) then write(string,2268) n_trains_tot,n_cars,itrain,icar,current,jobid 2268 format(i2.2,'x',i2.2,'_',i1.1,'_',i2.2,'_',f4.1,'_',i4.4,'.lat') txname = 'bmad_'//string call write_bmad_lattice_file(txname,ring) ! Write text data file for all elements for analysis by phase.kumac txlun = lunget() write(string,2270)n_trains_tot,n_cars,itrain,icar,current,jobid 2270 format(i2.2,'x',i2.2,'_',i1.1,'_',i2.2,'_',f4.1,'_',i4.4,'.txt') txname = 'txring/ele_'//string print *,' Writing file ',txname open(txlun, file=txname, status = 'unknown') do i=ring%n_ele_track,1,-1 write(txlun,2280)ring%param%total_length-ring%ele(i)%s, & ring%ele(i)%name(1:10), & ring%ele(i)%key, & ring%ele(ring%n_ele_track)%a%phi-ring%ele(i)%a%phi, & ! ring%ele(i)%a%phi, & ring%ele(i)%a%beta,ring%ele(i)%a%eta,& ring%ele(ring%n_ele_track)%b%phi-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), & co(i)%vec(1),-co(i)%vec(2), & co(i)%vec(3),-co(i)%vec(4), & -co(i)%vec(5),-co(i)%vec(6), & ring%ele(i)%value(sig_x$),ring%ele(i)%value(sig_y$) ! Changed signs of angles and dt 12 Dec 2005 enddo 2280 format(f11.6,2x,a10,i6,22f15.10) close(txlun) ! ! Write text data file for selected elements for analysis by phase.kumac txlun = lunget() write(string,2270)n_trains_tot,n_cars,itrain,icar,current,jobid txname = 'txring/ele_bbi_'//string print *,' Writing file ',txname open(txlun, file=txname, status = 'unknown') do i=ring%n_ele_track,1,-1 if ( ring%ele(i)%key .eq. beambeam$ ) & write(txlun,2280)ring%param%total_length-ring%ele(i)%s, & ring%ele(i)%name(1:10), & ring%ele(i)%key, & ring%ele(ring%n_ele_track)%a%phi-ring%ele(i)%a%phi, & ring%ele(i)%a%beta,ring%ele(i)%a%eta,& ring%ele(ring%n_ele_track)%b%phi-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), & co(i)%vec(1),-co(i)%vec(2), & co(i)%vec(3),-co(i)%vec(4), & -co(i)%vec(5),-co(i)%vec(6), & ring%ele(i)%value(sig_x$),ring%ele(i)%value(sig_y$) ! Changed signs of angles and dt 12 Dec 2005 enddo close(txlun) ! Write text data file for bpm elements for analysis by phase.kumac txlun = lunget() write(string,2270)n_trains_tot,n_cars,itrain,icar,current,jobid txname = 'txring/ele_bpm_'//string print *,' Writing file ',txname open(txlun, file=txname, status = 'unknown') do i=ring%n_ele_track,1,-1 if ( ring%ele(i)%key .eq. marker$ ) & write(txlun,2280)ring%param%total_length-ring%ele(i)%s, & ring%ele(i)%name(1:10), & ring%ele(i)%key, & ring%ele(ring%n_ele_track)%a%phi-ring%ele(i)%a%phi, & ring%ele(i)%a%beta,ring%ele(i)%a%eta,& ring%ele(ring%n_ele_track)%b%phi-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), & co(i)%vec(1),-co(i)%vec(2), & co(i)%vec(3),-co(i)%vec(4), & -co(i)%vec(5),-co(i)%vec(6), & ring%ele(i)%value(sig_x$),ring%ele(i)%value(sig_y$) enddo close(txlun) endif !***************************************************************** ! Chance to quit out before all hell breaks loose if (error) then write (6,3000) itrain,icar,current 3000 format(' BMAD tracking error found for train',i2, & ' car',i2,' current',f5.2,'-- skip this point.') goto 8000 endif ! Store first orbit as reference for bunch orbit data file if (current .eq. stepcurr ) then ! 1 nov 04 call reallocate_coord(orbit,ring%n_ele_max) orbit = co step = ring%param%total_length/ncopts do nstep = 1, ncopts s = (nstep-1)*step call twiss_and_track_at_s (ring, s, dumele, orbit, here) do j=1,6 firstorb(nstep,j) = here%vec(j) enddo enddo endif ! Print out IP BBI element do i = 1, ring%n_ele_track if ( ring%ele(i)%name(1:5) .eq. 'IP_CO' ) then call type_ele(ring%ele(i)) endif enddo ! Calculate beta values at IP do i = 1, ring%n_ele_track if ( ring%ele(i)%name(1:5) .eq. 'IP_L0' ) then alphaxstar = ring%ele(i)%a%alpha alphaystar = ring%ele(i)%b%alpha betaxstar = ring%ele(i)%a%beta betaystar = ring%ele(i)%b%beta etaxstar = ring%ele(i)%a%eta etaystar = ring%ele(i)%b%eta endif enddo ! Calculate maximum beta function value xbetamax = 0. sxmax=0. do i = 1, ring%n_ele_track if (ring%ele(i)%a%beta .ge. xbetamax ) then xbetamax = ring%ele(i)%a%beta sxmax = ring%ele(i)%s ixbetamax = i endif enddo ybetamax = 0. symax = 0. do i = 1, ring%n_ele_track if (ring%ele(i)%b%beta .ge. ybetamax ) then ybetamax = ring%ele(i)%b%beta symax = ring%ele(i)%s iybetamax = i endif enddo write(6,3245)itrain,icar,current, & xbetamax,ring%ele(ixbetamax)%name, & ybetamax,ring%ele(iybetamax)%name 3245 format(' Train ',i1,' car',i2,' current ',f4.1, & ' Max horiz beta ',f6.1,' at ',a25, & ' Max vert beta ',f6.1,' at ',a25) !=========================================================== ! Skip to next point if beta max values too large okbeta=.true. if( xbetamax .gt. 150 .or. ybetamax .gt. 150. ) then write (6,3250) itrain,icar,current 3250 format(' Beta values >150 m found for train',i2, & ' car',i2,' current',f5.2,'-- skip this point.') okbeta=.false. goto 8000 endif !=========================================================== ! Convert s values to modulo a half ring length, so that ! 0 and 768 are near each other if ( sxmax .gt. ring%param%total_length / 2. ) & sxmax = sxmax - ring%param%total_length if ( symax .gt. ring%param%total_length / 2. ) & symax = symax - ring%param%total_length ! Type out tunes ! write(6,5000)itrain,icar,current,ring%a%tune/twopi,ring%b%tune/twopi,ring%z%tune/twopi 5000 format('t',i1,'.b',i1,' e- tunes after turning on RF for ',f5.2,' e+ mA/bunch: Qx = ',f6.3, ' Qy = ',f6.3, ' Qz = ',f6.3) ! Calculate orbit difference position at IP if ( current .ne. 0. ) then dxip = co(0)%vec(1) - orbit(0)%vec(1) dyip = co(0)%vec(3) - orbit(0)%vec(3) else dxip = 0 dyip = 0 endif ! Calculate Welch/Temnykh B parameter summed over parasitic crossings ! Following ST presentation at CESR-c MiniMAC 7/23/05 ! B must be less than 1.2 +- .2 to ensure adequate beam lifetime (50 minutes) ! This calculation only works if all BBIWT values are unity, since ! the positron train and bunch number for each of the crossings is ! not available here. We assume equally populated positron bunches. bparam = 0. do i = 1, ring%n_ele_track if (ring%ele(i)%key .eq. beambeam$ .and. & ring%ele(i)%name(1:5) .ne. "IP_L0" .and. & ring%ele(i)%name(1:6) .ne. "IP_COL" .and. & abs(ring%ele(i)%s - 384.21) .gt. 0.2 ) then xsig2 = mode%a%emittance * ring%ele(i)%a%beta + & ( co(i)%vec(6) * ring%ele(i)%a%eta )**2 xdel2 = diff_co(i)%vec(1)**2 ! B is weighted with total beam current (assume equal e+/e- currents) and inverse of beam energy to account for beam rigidity ! Factor 10 is from arbitrary convention to give 1.2 as limiting value ! bcontribution = (current/10.) * & bcontribution = 0.1 * & (5.3e9/ring%ele(0)%value(e_tot$)) * & ring%ele(i)%b%beta / (xdel2/xsig2) bparam = bparam + bcontribution**2 ! write(6,5090)i,ring%ele(i)%s,ring%ele(i)%b%beta,sqrt(xdel2),sqrt(xsig2),bcontribution 5090 format('B parameter: index, s, betay, xdel, xsig, B =', & i5,f8.2,f6.2,2f7.4,e12.5) endif enddo bparam = sqrt(bparam) ! write(6,5091)ring%E_TOT,current,bparam write(6,5091)ring%ele(0)%value(e_tot$),current,bparam 5091 format(' Beam energy, e+ bunch current, global B parameter:',e12.5,f10.3,e12.5) ! !=========================================================== ! Calculate global ring parameters call radiation_integrals( ring, co, mode ) ! Emittance for electrons ! per mil sige_e = 1000 * mode%sige_e ! cm sig_z = 100 * mode%sig_z ! MeV/turn e_loss = mode%e_loss / 10**6 ! nanometers xemm = 10**9 * mode%a%emittance yemm = 10**9 * mode%b%emittance zemm = 10**9 * mode%z%emittance xalphad = mode%a%alpha_damp yalphad = mode%b%alpha_damp zalphad = mode%z%alpha_damp ! milliseconds ! Fixed this 7 July 2005 (I thought I fixed this months ago!) ! xtau = 10*3 * xalphad / (ring%param%total_length/c_light) ! ytau = 10*3 * yalphad / (ring%param%total_length/c_light) ! ztau = 10*3 * zalphad / (ring%param%total_length/c_light) xtau = 1000 * (ring%param%total_length/c_light) / xalphad ytau = 1000 * (ring%param%total_length/c_light) / yalphad ztau = 1000 * (ring%param%total_length/c_light) / zalphad ! print *, 'C, c, alphax=',ring%param%total_length,c_light,xalphad !=========================================================== if( xemm .le. 0 .or. yemm .le. 0 ) then write (6,3300) itrain,icar,current 3300 format(' Negative emittance found by RADIATION_INTEGRALS for train',i2, & ' car',i2,' current',f5.2,'-- skip this point.') goto 8000 endif !=========================================================== ! Write data to file for analysis by synclight.kumac ! NOTE: i=2 gives the ELECTRON orbit at the POSITRON sync light source point do i=1,2 call twiss_and_track_at_s (ring, ssource(i), dumele, co, here) svec(7*(i-1)+1) = ssource(i) svec(7*(i-1)+2) = 100*here%vec(1) svec(7*(i-1)+3) = 1000*here%vec(2) svec(7*(i-1)+4) = 100*here%vec(3) svec(7*(i-1)+5) = 1000*here%vec(4) svec(7*(i-1)+6) = dumele%a%beta svec(7*(i-1)+7) = dumele%b%beta enddo write(lundat3,5060)itrain,icar,current, svec 5060 format(2i5,15f12.5) ! Write IP element data to file for analysis by ip.kumac do i=1, ring%n_ele_track ip = ring%ele(i)%name == 'IP_COLLISION' if(ip) then write(lundat4,5070)itrain,icar,current, & 1e6*ring%ele(i)%value(sig_x$),1e6*ring%ele(i)%value(sig_y$),1e6*ring%ele(i)%value(sig_z$), & 1e2*ring%ele(i)%value(x_offset$),1e3*ring%ele(i)%value(x_pitch$), & 1e2*ring%ele(i)%value(y_offset$),1e3*ring%ele(i)%value(y_pitch$), & 1e3*ring%ele(i)%value(tilt$) 5070 format(2i5,15f12.5) write(6,5071)itrain,icar,current, & 1e6*ring%ele(i)%value(sig_x$),1e6*ring%ele(i)%value(sig_y$),1e6*ring%ele(i)%value(sig_z$), & 1e2*ring%ele(i)%value(x_offset$),1e3*ring%ele(i)%value(x_pitch$), & 1e2*ring%ele(i)%value(y_offset$),1e3*ring%ele(i)%value(y_pitch$), & 1e3*ring%ele(i)%value(tilt$) 5071 format(' To IP data file:',2i5,15f12.5) cycle endif enddo ! Write data to file for analysis by bunch.kumac write(lundat,5100)itrain,icar,current, & ring%a%tune/twopi, & ring%b%tune/twopi, & ring%z%tune/twopi, & ! ixbetamax, & xbetamax, & ! iybetamax, & ybetamax, & sxmax,symax, & 1000*dxip,1000*dyip, & sige_e, sig_z, & xemm, yemm, zemm, & ! xtau, ytau, ztau, & alphaxstar, alphaystar, & betaxstar, betaystar, & etaxstar, etaystar, & bparam, & diff_co(0)%vec(2),diff_co(0)%vec(4),& 1000* diff_co(0)%vec(1),1000*diff_co(0)%vec(3) !5100 format(2i5,4f12.5,i5,f12.5,i5,f12.5,2f12.5,18f12.5) 5100 format(2i5,28f12.5) write(6,5150)itrain,icar,current, & ring%a%tune/twopi, & ring%b%tune/twopi, & ring%z%tune/twopi, & ! ixbetamax, & xbetamax, & ! iybetamax, & ybetamax, & sxmax,symax, & 1000*dxip,1000*dyip, & sige_e, sig_z, & xemm, yemm, zemm, & ! xtau, ytau, ztau, & alphaxstar, alphaystar, & betaxstar, betaystar, & etaxstar, etaystar, & bparam, & diff_co(0)%vec(2),diff_co(0)%vec(4),& 1000*diff_co(0)%vec(1),1000*diff_co(0)%vec(3) !5150 format(' To BUNCH Data File:',2i5,4f12.5,i5,f12.5,i5,f12.5,2f12.5,17f12.5) 5150 format(' To BUNCH Data File:',2i5,28f12.5) lastcurrent = current ! Calculate dynamic aperture for integer positron currents if(.not. error) then if( calc_dynap .and. mod(current,1.) .eq. 0. ) then lastentry = .false. orb0 = co(0) call dynap ( ring, orb0, jobid, lastentry, itrain, icar, current, mode ) endif else write (6,3800) itrain,icar,current 3800 format(' BMAD tracking error found for train',i2, & ' car',i2,' current',f5.2,'-- skip DYNAP.') endif 8000 continue ! Write out lattice file if ( itrain .eq. 1 .and. icar .eq. 3 .and. abs(current-2.4) .le. 0.1 ) then ! print *,' Writing digested lattice file ...' ! call write_digested_bmad_file('digested_t1b3_2p4ma.lat',ring) ! call write_bmad_lattice_file('t1b3_2p4ma.lat',ring) endif enddo ! end loop over e+ bunch current ! Calculate dynamic aperture for maximum current and close dynap file ! if last bunch ! Don't calculate if integer current, since it was already done inside ! the current loop lastentry = ( itrain .eq. itrainmax .and. icar .eq. icarmax ) if ( calc_dynap .and. okbeta .and. .not. error .and. & mod(current,1.) .ne. 0. ) then orb0 = co(0) call dynap ( ring, orb0, jobid, lastentry, itrain, icar, lastcurrent, mode ) endif ! Write file of orbit differences for this maximum current setting ! Loop over s in steps of length STEP (meters) orbit = co do nstep = 1, ncopts s = (nstep-1)*step call twiss_and_track_at_s (ring, s, dumele, orbit, here) do j=1,6 difforb(nstep,j) = here%vec(j) - firstorb(nstep,j) enddo write(lundat2,6000)itrain,icar,lastcurrent,s,(difforb(nstep,j),j=1,6), (firstorb(nstep,j),j=1,6) 6000 format(2i5,f9.4,f9.4,2(6e14.5)) enddo ! At highest current, print out all BBI, including IP BBI, and their strengths ! write(6,4000)itrain,icar 4000 format(/' Beam-beam elements for e- train ',i2,' bunch',i2,':') do i=1,ring%n_ele_track if ( ring%ele(i)%key == beambeam$ ) then ! write(6,4100)ring%ele(i)%name,ring%ele(i)%s, & ! ring%ele(i)%value(charge$), & ! ring%ele(i)%value(bbi_const$) 4100 format(a12,f10.2,' charge=',f6.3,' bbi_const= ',e14.6) endif enddo 8500 continue deallocate(co) deallocate(diff_co) deallocate(orbit) enddo ! end loop over test car enddo ! end loop over test train close ( lundat ) close ( lundat2 ) close ( lundat3 ) close ( lundat4 ) deallocate (closed_orbit) call date_and_time(date, time) write(6,9100)date,time 9100 format(' ==> Exiting Program BUNCH on ',a8,' at ',a10) end program bunch_bbi !################################################################################### subroutine readjobid ! ! Read job ID number from file ! and store it in module JOBID ! use jobidmod character*100 fname fname='/nfs/acc/temp/critten/bunch_bbi/job_id_number' open(26,file=fname,err=100,status='old') read(26,1000)jobid 1000 format(i10) write(6,1100)fname,jobid 1100 format(' Routine READJOBID opened file ',a100, & ' Job number is ',i4) close(26) goto 900 100 continue print *,' Routine READJOBID could not open file ',fname jobid=1 900 continue ! end subroutine readjobid !################################################################################### subroutine writejobid(jobid) ! ! Write job number from file ! character*100 fname integer jobid fname='/nfs/acc/temp/critten/bunch_bbi/job_id_number' fname='job_id_number' open(26,file=fname,err=100,status='old') write(26,1000)jobid 1000 format(i10) close(26) write(6,1100)jobid,fname 1100 format(' Routine WRITEJOBID incremented job number to ',i4, & ' in file ',a20) goto 999 100 continue print *,' Routine WRITEJOBID could not open file ',fname ! 999 continue ! end subroutine writejobid !################################################################################### subroutine set_custom_pretzel (ring, co, off, prz1, prz13 ) ! ! Set horizontal separator kick values ! Leave at design values if set value is zero. ! This routine finds the horizontal separators ! in the ring given as input. ! It makes assumptions about how much the kicks ! are changed in radians by prz1 and prz13 in cu. ! The kicks are scaled with the crossing angle, ! which is why the closed orbit is required as ! an input paramter. ! If OFF is true, the prz values are ignored and ! the separator kicks are set to zero. The kicks ! are stored in order to be restored on the ! next call with OFF set to false. ! ! OFFCALL is a flag indicating that the previous ! call had OFF equal true. ! ! 1 Sep 2004 JAC ! ! Created a local index so that the local array ! of separators always has the order 8W,45W,45E,8E. ! 2 nov 2004 jac !=============================================== use bmad implicit none type (lat_struct) ring type (coord_struct), allocatable :: co(:) ! type (coord_struct) co(0:ring%n_ele_max) logical off real(rp) prz1,prz13 integer ihs, nhs, i integer ix(4) real(rp) hkick(4),scale(4) real(rp), save :: xangle_ip, kdiff real(rp), save :: storedkick(4) real(rp), save :: designkick(4) logical, save :: offcall = .false. logical, save :: firstcall = .true. print *,' SET_CUSTOM_PRETZEL called with off, prz1,prz13=',off, prz1,prz13 ! get location of separators nhs = 0 do i=1,ring%n_ele_max if ( ring%ele(i)%key /= elseparator$) cycle if( ring%ele(i)%name(1:1) /= 'H' ) cycle if (.not. attribute_free (ring%ele(i), 'HKICK')) cycle if ( nhs .lt. 4) then nhs = nhs + 1 if( ring%ele(i)%name(7:9) == '08W' ) then ihs=1 elseif( ring%ele(i)%name(7:9) == '45W' ) then ihs=2 elseif( ring%ele(i)%name(7:9) == '45E' ) then ihs=3 elseif( ring%ele(i)%name(7:9) == '08E' ) then ihs=4 else ihs=1 print *,'SET_CUSTOM_PRETZEL Warning: Unknown separator found:', ring%ele(i)%name endif ix(ihs) = i else print *,' SET_CUSTOM_PRETZEL Warning: More than 4 horiz seps found!' endif print *,' Found separator ',ring%ele(i)%name, & ' ring index= ',i, & ' local index= ',ihs if (offcall) then hkick(ihs) = storedkick(ihs) else hkick(ihs) = ring%ele(i)%value(hkick$) endif enddo if( nhs /= 4 ) then print *,' Rtn SET_CUSTOM_PRETZEL:": Did not find four horizontal separators' stop endif if (firstcall) then ! This is the design crossing angle. designkick=hkick xangle_ip = co(0)%vec(2) print *,' Rtn SET_CUSTOM_PRETZEL finds design horiz crossing angle at IP:',xangle_ip firstcall = .false. endif ! Now set kicks according to input values for prz1 and prz13 ! Scale the separator strengths (seps 1, 2, 5, 6 ) ! Separator 8W (smaller than 8E for negative PRZ13) scale(1:4) = 1. if ( off ) then scale(1:4) = 0. offcall = .true. storedkick = hkick else offcall = .false. if ( prz1 /= 0. .and. prz13 /= 0. ) then scale(1) = (prz1 + prz13) / ( 1.e6 * xangle_ip ) scale(2) = prz1 / ( 1.e6 * xangle_ip) scale(3) = prz1 / ( 1.e6 * xangle_ip) scale(4) = (prz1 - prz13) / ( 1.e6 * xangle_ip) elseif ( prz1 /= 0. .and. prz13 == 0. ) then ! Change only PRZ1: Preserve 8W/8E kick difference kdiff = designkick(2) - designkick(3) scale(1) = kdiff / ( 2.*hkick(1) ) + prz1 / ( 1.e6 * xangle_ip) scale(2) = prz1 / ( 1.e6 * xangle_ip) scale(3) = prz1 / ( 1.e6 * xangle_ip) scale(4) = - kdiff / ( 2.*hkick(4) ) + prz1 / ( 1.e6 * xangle_ip) ! Lack of symmetry complicates calculation. For now, assume PRZ13-asymmetry is zero, ! as is true for present lattice design conventions. ! Then all separator kicks scale the same. 3 May 2005 JAC ! If xangle_ip=-3.5 mrad, then prz1=3500 will leave the angle unchanged. scale(1) = - prz1 / ( 1.e6 * xangle_ip) scale(2) = - prz1 / ( 1.e6 * xangle_ip) scale(3) = - prz1 / ( 1.e6 * xangle_ip) scale(4) = - prz1 / ( 1.e6 * xangle_ip) elseif ( prz1 == 0. .and. prz13 /= 0. ) then ! Change only PRZ13: Change only difference between 8W and 8E scale(1) = 1 + prz13 / ( 1.e6 * xangle_ip) scale(4) = 1 - prz13 / ( 1.e6 * xangle_ip) endif endif print *,'Scale is ',scale(1:4) do i=1,4 !********************************************** ! print *,'22 Aug 06 Emergency hardwire prz13 for ELCSR' !********************************************** if(i.eq.1)then ! scale(i)=0.233e-3/designkick(i) endif if(i.eq.4)then ! scale(i)=-0.221e-3/designkick(i) endif ring%ele( ix(i) )%value(hkick$) = scale(i) * designkick(i) print *,ring%ele( ix(i) )%name,' Design kick = ', designkick(i), & ' New kick = ', scale(i) * designkick(i) ! call type_ele(ring%ele( ix(i) ),.false.,0,.false., 0, .true., ring ) call lat_make_mat6(ring, ix(i) , co ) ! needed in case pc is in sep enddo end subroutine set_custom_pretzel !################################################################################### subroutine dynap ( ring, orb0, jobid, lastentry, itrain, icar, current, mode ) ! ! Calculate dynamic aperture for this ring ! ! J.A.Crittenden 8/16/04 !=================================================================== use dynamic_aperture_mod use bmad use bmadz_mod use bmadz_interface implicit none ! Input arguments type (lat_struct) ring type (coord_struct) orb0,orbip logical lastentry integer jobid integer itrain, icar real(rp) current type (normal_modes_struct) mode type (coord_struct), allocatable :: orbit(:) type (aperture_data_struct) aperture type (aperture_param_struct) aperture_param ! integer, save :: i_e, i_xy, it, i, turn_lost, ixr, i_e_max, ix integer, save :: n_turn, n_xy_pts, point_range(2), n_energy_pts real(rp), save :: x_init, y_init real(rp), save :: energy(10) real(rp), save :: accuracy, aperture_multiplier real(rp) e_init, theta ! real(rp), save :: sigx, sigy, sige, sigxp real(rp), save :: alpha_ip, beta_ip ! character fname*40 character datfname*40 integer, save :: lun, lundat, lundat2 integer n ! integer ienter /0/ ! namelist / dynap_setup / n_turn, n_xy_pts, point_range, n_energy_pts, & x_init, y_init, & energy, accuracy, aperture_multiplier write(6,1000) jobid, lastentry, itrain, icar, current, orb0%vec(1:6) 1000 format(' Rtn DYNAP called with jobid, lastentry = ',i5,l, & ' itrain, icar, current = ', 2i4,f7.2,/10x, & ' closed orbit = ',6f14.5) if ( ienter .eq. 0 ) then ienter= ienter + 1 ! Read in setup info lun = lunget() fname='dynap_setup.in' open(lun,file=fname, status= 'old') read(lun, nml=dynap_setup) close(lun) ! aperture_param%n_turn = n_turn aperture_param%x_init = x_init aperture_param%y_init = y_init aperture_param%accuracy = accuracy ! write(6,1200)fname, n_turn, n_xy_pts, point_range, & x_init, y_init, & accuracy, aperture_multiplier, & n_energy_pts, (energy(i),i=1,n_energy_pts) 1200 format(' ==> Routine DYNAP initialized with file ',a40/, & ' Number of test turns: ',i5/, & ' Number of X,Y points: ',i4/, & ' Point range: ',2i6/, & ' First X value: ',f7.3/, & ' First Y value: ',f7.3/, & ' Chosen accuracy: ', f9.5/, & ' Aperture multiplier:', f4.1/, & ' Number of energy points: ',i4/, & ' Energy offsets: ',10f9.5) !======================================================= ! For first call, store N_TURN electron orbits ! with no X,Y offset wrt the closed_orbit, but with ! the first chosen energy offset if there is one. ! Open output orbit file lundat2=lunget() write(datfname,1050)jobid 1050 format('dat/dynap_orbit_',i4.4,'.dat') open(lundat2, file=datfname, status = 'unknown') ! allocate storage call reallocate_coord (orbit, ring%n_ele_track) orbit(0) = orb0 print *,' DYNAP input argument ORB0%vec(6)',orb0%vec(6) ! Track with first chosen energy offset ! Without any offset from the closed orbit, the orbit ! will not change from turn to turn ! if ( n_energy_pts .gt. 0 ) orbit(0)%vec(6) = energy(1) ! ! Fixed bug which set the initial energy equal to the offset ! rather than to the closed_orbit energy plus the offset 15 Nov 2006 jac if ( n_energy_pts .gt. 1 ) orbit(0)%vec(6) = orb0%vec(6)+energy(2) ! Prevent termination of tracking on aperture violations bmad_com%aperture_limit_on=.false. do n=0,n_turn ! print *,' DYNAP: track_many called with initial energy =',orbit(0)%vec(6) ! Track electrons forward in reversed ring call track_many (ring, orbit, 0, 0, +1) do i=1, ring%n_ele_track ! Write complete orbit only for turns 0-9, otherwise only IP if ( i.eq.1 .or. n.lt. 10 ) & write(lundat2,100) n, ring%ele(i)%s,ring%ele(i)%key, & 1000.*orbit(i)%vec(1),1000.*orbit(i)%vec(2), & 1000.*orbit(i)%vec(3),1000.*orbit(i)%vec(4), & 1000.*orbit(i)%vec(5),1000.*orbit(i)%vec(6) 100 format(2x, i6, f10.4, i6, 2(1x,f10.5),2(1x,f12.7),1x,f8.3,1x,f11.6) enddo enddo close (lundat2) !================================================================== ! Open output data file lundat=lunget() write(datfname,1100)jobid 1100 format('dat/dynap_',i4.4,'.dat') open(lundat, file=datfname, status = 'unknown') ! endif !First entry ! bmad_com%aperture_limit_on = .true. ! ! Calculate beam size at IP and energy spread do i=1, ring%n_ele_max if( ring%ele(i)%name == 'IP_L0' ) then ! sigx = sqrt( ring%ele(i)%a%beta * mode%a%emittance ) ! sigy = sqrt( ring%ele(i)%b%beta * mode%b%emittance ) ! Include energy spread contribution and use assumption of full coupling for ! vertical beam size 27 dec 2005 jac sige = mode%sige_e sigx = sqrt( ring%ele(i)%a%beta * mode%a%emittance + (ring%ele(i)%a%eta * sige)**2 ) sigy = sqrt( ring%ele(i)%b%beta * mode%a%emittance / 2 + (ring%ele(i)%b%eta * sige)**2 ) sigxp = sqrt( mode%a%emittance / ring%ele(i)%a%beta ) alpha_ip = ring%ele(i)%a%alpha beta_ip = ring%ele(i)%a%beta endif enddo i_e_max = max(1, n_energy_pts) print *,' i_e_max=',i_e_max do i_e = 1, i_e_max ! Fixed bug which set the initial energy equal to the offset ! rather than to the closed_orbit energy plus the offset 15 Nov 2006 jac ! e_init = energy(i_e) orbip = orb0 e_init = orb0%vec(6)+energy(i_e) orbip%vec(6) = e_init ! Moved from dynamic_aperture. This should really be replaced with closed_orbit_calc. orbip%vec(1) = orbip%vec(1) + e_init*ring%ele(0)%x%eta orbip%vec(3) = orbip%vec(3) + e_init*ring%ele(0)%y%eta orbip%vec(2) = orbip%vec(2) + e_init*ring%ele(0)%x%etap orbip%vec(4) = orbip%vec(4) + e_init*ring%ele(0)%y%etap print *,'DYNAP calculating for initial energy =',e_init ! Fixed horizontal angle ! orb0%vec(2) = sigxp ! print *,' Routine DYNAP using 1-sigma horizontal launch angle: ',sigxp do i_xy = point_range(1), point_range(2) theta = (i_xy - 1) * pi / max(1, n_xy_pts - 1) !15 nov 06 ! call dynamic_aperture1 (ring, orb0, theta, aperture_param, aperture) call dynamic_aperture1 (ring%branch(0), orbip, theta, aperture_param, aperture) write ( lundat, '(2i4,f7.2,f9.5,2f11.6, i7, 3f14.9, i5, f10.3)') & itrain, icar, current, & e_init, & aperture%x, aperture%y, & aperture%i_turn, & sigx, sigy, sige, & aperture%plane, & ring%param%total_length-ring%ele(aperture%ix_lat)%s write (6, '(2i4,f7.2,f9.5,2f11.6, i7, 3f14.9, i5, f10.3)') & itrain, icar, current, & e_init, & aperture%x, aperture%y, & aperture%i_turn, & sigx, sigy, sige, & aperture%plane, & ring%param%total_length-ring%ele(aperture%ix_lat)%s enddo ! i_xy enddo ! i_e ! if ( lastentry ) close ( lundat ) ! end subroutine dynap !################################################################################### subroutine ringmod ( ring, fname, fact ) ! ! Modify ring elements according to information ! in file fname. The file has a custom format ! developed by MGB for the purpose of BBI ! phase corrections. ! ! Corrections are weighted by factor FACT ! ! 23 November 2004 JAC use bmad implicit none ! Input arguments type (lat_struct) ring character*40 fname real(rp) fact ! ! integer lun,nhead,i,istat,nmod ! character*16 ename real(rp)kcorr real(rp) kcorrv(ring%n_ele_track) character*80 rec ! !***************************************** ! Read file ! write(6,1000)fact,fname 1000 format('Rtn RINGMOD called with factor ',f5.2, & ' for quad correction file ',a40) kcorrv(1:ring%n_ele_track)=0. lun = lunget() open(lun,file=fname, status= 'unknown', err=950, iostat=istat) read (lun,'(i4)')nhead do i=1,nhead read (lun,'(a80)')rec enddo nmod=0 10 continue read(lun,2000,err=500,end=900, iostat=istat)ename,kcorr 2000 format(2x,a16,2x,e16.6) ! Store quad k value do i=1,ring%n_ele_track if (ring%ele(i)%name(1:4) .eq. ename) then kcorrv(i) = kcorr*fact if(kcorrv(i).ne.0.)then nmod=nmod+1 print *,' Ringmod correction',nmod,ename,kcorr*fact endif endif enddo goto 10 500 continue print *,' Rtn RINGMOD read error. istat=',istat 900 continue close(lun) write(6,9000)fname,nmod 9000 format(' Rtn RINGMOD closing file ',a20,& ' after storing ',i5,' modifications') goto 999 950 continue print *,' Rtn RINGMOD open file error. istat=',istat close(lun) 999 continue ! Modify ring elements according to stored array do i=1,ring%n_ele_track ring%ele(i)%value(k1$) = & ring%ele(i)%value(k1$) + kcorrv(i) enddo end subroutine ringmod !################################################################################### subroutine pcbetcor(ring, co, fact) ! ! Set PCBETING knobs if available ! ! 2 November 2005 JAC use bmad implicit none ! Input arguments type (lat_struct) ring type (coord_struct), allocatable :: co(:) ! Input argument e+ bunch current in mA real(rp) fact integer i integer ik,ix character*16 knobname integer, parameter :: maxknobs = 108 integer indknob(maxknobs) integer ienter /0/ integer, save :: nknobs integer knobnum(maxknobs) /9,10,11,12,13,14,15,16,100*0/ real knobcmd(maxknobs) /0.144, -0.056, 0.501, -0.110, -1.570, -0.187, 0.,0., 100*0./ ! real knobcmd(maxknobs) /0.,0.,0.,0.,0.,0./ interface subroutine readknobs(nknobs,knobnum,knobcmd) use bmad implicit none integer nknobs integer, intent(out) :: knobnum(:) real, intent(out) :: knobcmd(:) end subroutine readknobs end interface print *,'Entering Rtn PCBETCOR with correction factor ',fact ienter = ienter + 1 if ( ienter .eq. 1) call readknobs(nknobs,knobnum,knobcmd) if (nknobs.eq.0) then print *,'Rtn PCBETCOR exiting because no knobs found' goto 900 endif print *,' Rtn pcbetcor called with NKNOBS=',nknobs do ik=1,nknobs if(knobnum(ik).le.9)then write(knobname,2000)knobnum(ik) 2000 format('RAW_PCBETING_'i1) elseif(knobnum(ik).le.99)then write(knobname,2100)knobnum(ik) 2100 format('RAW_PCBETING_'i2) else write(knobname,2200)knobnum(ik) 2200 format('RAW_PCBETING_'i3) endif call element_locator(knobname,ring,ix) if(ix.le.0) then print *,'Rtn PCBETCOR could not find element ',knobname, & ' Exiting without correction' goto 900 endif !if(ik.eq.1)then ! Turn off BBI in East (i.e. West in reversed ring!) if RAW_PCBETING_9 found ! do i=1,ring%n_ele_track ! if ( ring%ele(i)%key == beambeam$ ) then ! if ( ring%ele(i)%s < ring%param%total_length/2 ) then ! ring%ele(i)%value(charge$) = 0. ! write(6,1000)ring%ele(i)%s !1000 format('Rtn PCBETCOR turning off BBI at s=',f10.2) ! endif ! endif ! enddo !endif indknob(ik)=ix enddo do ik=1,nknobs ring%ele(indknob(ik))%control%var(1)%value = knobcmd(ik) * fact ! print *,' Rtn PCBETCOR sets element ',ring%ele(indknob(ik))%name, & ! ' to cmd ', ring%ele(indknob(ik))%control%var(1)%value call lat_make_mat6( ring, -1, co ) enddo 900 continue end subroutine pcbetcor !################################################################################### subroutine readknobs(nknobs,knobnum,knobcmd) ! ! Read PCBETA knob command file ! ! Output: NKNOBS -- Number of knobs found ! KNOBNUM -- Knob numbers in knob name ! KNOBCMD -- Command values ! use bmad implicit none integer nknobs integer, intent(out) :: knobnum(:) real, intent(out) :: knobcmd(:) character*50 knobcomfn integer lun,i,nkn,icmd character*80 string integer arr(100) integer nrec lun = lunget() nknobs=0 ! knobcomfn='bbi_beta_00193.com' ! 20051220 lattice 8x4 Charge smeared over 9 trains ! knobcomfn='bbi_beta_00209.com' ! 20051220 lattice 8x4 Charge not smeared over 9 trains ! knobcomfn='bbi_beta_00211.com' ! 20051220 lattice 8x3 Charge smeared over 9 trains ! knobcomfn='bbi_beta_00212.com' ! ! 20060213 lattice 8x3 Charge smeared over 9 trains ! knobcomfn='bbi_beta_00217.com' ! Same as above except additional compensation for IP BBI added ! Delta's for 1-6 are -600, -700, 2050, -900, -900, 0. ! knobcomfn='bbi_beta_00217ip.com' ! ! 20060213 lattice 8x3 Charge smeared over 9 trains ! prz1=2850 ! knobcomfn='bbi_beta_00218.com' ! Same as above except additional compensation for IP BBI added ! Delta's for 1-6 are -600, -700, 2050, -900, -900, 0. ! as determined in April ! In July, we measured: -1315 -1350 5300 -3500 -3200 ! knobcomfn='bbi_beta_00218ip.com' ! ! 20060213 lattice 8x4 Charge smeared over 9 trains. correct for bunch 2. ! prz1=2850 ! knobcomfn='bbi_beta_00219.com' ! Same as above except additional compensation for IP BBI added ! Delta's for 1-6 are -600, -700, 2050, -900, -900, 0. ! knobcomfn='bbi_beta_00219ip.com' ! ! 20060615 lattice 8x4 Charge smeared over 9 trains. correct for bunch 2. ! Sextupoles re-optimized at 0 mA ! prz1=2850 (design value for this lattice) ! knobcomfn='bbi_beta_00226.com' ! Same as above except additional compensation for IP BBI added ! Delta's for 1-6 are -600, -700, 2050, -900, -900, 0. ! knobcomfn='bbi_beta_00226ip.com' ! ! 20060616 lattice 8x4 Charge smeared over 9 trains. correct for bunch 2. ! Sextupoles re-optimized at 2 mA ! prz1=2850 (design value for this lattice) ! knobcomfn='bbi_beta_00227.com' ! Same as above except additional compensation for IP BBI added ! Delta's for 1-6 are -600, -700, 2050, -900, -900, 0. ! knobcomfn='bbi_beta_00227ip.com' ! ! 20060213 lattice e+ t2-6.b1-5 for MS 22 Aug 2006 ! HEP condx ! knobcomfn='bbi_beta_00228.com' ! ! 20060213 lattice e+ t2-6.b1-5 for MS 22 Aug 2006 ! ELCSR condx ! ! fit_postop_20070116_qt_beta_seps. 8x3 design pretzel & tunes. For t1.b2. ! No IP correction ! knobcomfn='bbi_beta_00301.com' ! fit_postop_20070116_qt_beta_seps. 8x3 design pretzel & tunes. For t1.b2. ! With IP correction ! knobcomfn='bbi_beta_00301ip.com' ! scsk01_zero_20070115-1 lattice. 8x3 design pretzel and tunes. For t1.b2. ! knobcomfn='bbi_beta_00241.com' ! ! fit_postop_20070116. 8x3 design pretzel & tunes. For t1.b2. ! No IP correction ! knobcomfn='bbi_beta_00252.com' ! fit_postop_20070116_qt_beta_seps. 8x3 design pretzel & tunes. For t1.b2. ! With IP correction knobcomfn='bbi_beta_00252ip.com' ! With MGB's by-hand calculated IP correction ! knobcomfn='bbi_beta_00252ip_mgb.com' ! With HALF of MGB's by-hand calculated IP correction ! knobcomfn='bbi_beta_00252ip_mgb2.com' ! With HALF of MGB's by-hand calculated IP correction for horizontal ! and one quarter of the calculated correction for vertical ! knobcomfn='bbi_beta_00252ip_mgb3.com' ! Above mgb corrections have wrong sign in coeffs 4-6 ! Corrected files are mgb4 (full), mgb5 (half), mgb6 (half horiz, quarter vert) ! knobcomfn='bbi_beta_00252ip_mgb4.com' ! knobcomfn='bbi_beta_00252ip_mgb5.com' ! knobcomfn='bbi_beta_00252ip_mgb6.com' ! ! Testing IP coeffs alone. 4/11/07 ! mgb4: wrong sign in coeffs 4-6, and coeff 5 factor 10 too small ! knobcomfn='bbi_beta_ip_mgb4.com' ! mgb: mgb full ip correction coefficients as intended ! knobcomfn='bbi_beta_ip_mgb.com' ! ! 3770_inj_20070606_2p0. 8x3 design pretzel & tunes. For t1.b2. ! No IP correction, LRBBI correction alone ! knobcomfn='bbi_beta_00254.com' ! IP correction alone (from 1928 (3 mA) and 1935 (1 mA)) ! knobcomfn='bbi_beta_1935ip.com' ! IP and LRBBI corrections superposed ! knobcomfn='bbi_beta_00254ip.com' ! IP and LRBBI corrections superposed, IP correction half strength ! knobcomfn='bbi_beta_00254ip2.com' ! IP and LRBBI corrections superposed, IP correction quarter strength knobcomfn='bbi_beta_00254ip4.com' open(lun,file=knobcomfn,err=100,status='old') write(6,1100)knobcomfn 1100 format(' Routine READKNOBS opened file ',a50) nknobs=0 nrec=0 10 continue read(lun,1000,end=800)string read(lun,1000,end=800) 1000 format(a80) nrec=nrec+1 ! print *,' READKNOBS read string',string read(string(25:80),*,end=20)arr 20 continue ! print *,' READKNOBS reads ARR:',arr(1:10) nkn=arr(2)-arr(1)+1 do i=1,nkn icmd=arr(2+i) if ( icmd .ne. 0 .and. nknobs .lt. 108 ) then nknobs=nknobs+1 knobnum(nknobs)=arr(1)+i-1 knobcmd(nknobs)=arr(2+i)/1000. endif enddo ! Read 18 records maximum if( nrec .lt. 18 )goto 10 800 continue close(lun) write(6,8000)nknobs 8000 format(' Rtn READKNOBS reached EOF. NKNOBS=',i4) goto 900 100 continue print *,' Routine READKNOBS could not open file ',knobcomfn 900 continue end subroutine readknobs !################################################################################### subroutine collide (ring, i_dim) ! ! 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 use bmad_interface use bmadz_mod use bmadz_interface 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(rp) Qx, Qy, Qx_in , Qy_in real(rp) rgamma, den, xi_v, xi_h real(rp) r_0/2.81794092e-15/ logical error 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_in%vec(1) = 2.5e-4 final_pos_in%vec(1) = 4e-4 final_pos_out%vec(:) = 0 write(6,1000)final_pos_in%vec(1:6) 1000 format(' Rtn COLLIDE imposing closure criterion coordinate vector:', & 6(2x,e10.2)) 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, error) print*,'Rtn COLLIDE: after close pretzel but before close vertical: ' print '(1x,3(a9,f12.4))',' Qx = ',ring%a%tune/twopi,' Qy = ',ring%b%tune/twopi,' Qz = ',ring%z%tune/twopi ! No vert seps 4 Sep 2009 ! 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, error) print * print *,'Rtn COLLIDE: After CLOSE VERTICAL ' print '(1x,3(a9,f12.4))',' Qx = ',ring%a%tune/twopi,' Qy = ',ring%b%tune/twopi,' Qz = ', ring%z%tune/twopi print '(a18,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 print '(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, error) print '(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 !###################################################################################!+ ! subroutine LRBBI_SETUP_CUSTOM (lat_in, lat_out, particle, con,pc, j_car, current, rec_taylor, bbiwt) ! ! Subroutine to add lrbbi elements to lat at relevant crossing points from custom bunch patterns ! ! Input: ! LAT_IN -- lat_struct: Lat containing the lattice to be modified ! particle -- Integer : (+-1) for positrons or electrons ! j_car -- Integer : Car of test particle ! con -- constraint_struct: contains total number of trains and cars ! pc -- Pc_struct : contains the parasitic crossings information ! current -- real : bunch current in mA ! rec_taylor -- Logical : passed to stash_split_lat. save taylor as digested file ! bbiwt(:,:) -- Real, optional : Weights for the current in each ! strong bunch. Range 0 to 1. ! The dimensions should be defined by the calling ! routine to be (n_trains_Tot,n_cars) ! ! Output: ! LAT_OUT -- lat_struct: Lat_out includes Beam-beam elements ! at each of crossing points for bunch i_train, j_car ! subroutine lrbbi_setup_custom(lat_in, lat_out, particle, con, pc, j_car, current, rec_taylor) use bmad use bmadz_mod use cesr_mod use bunchcross_mod use constraints_mod implicit none interface subroutine MARK_LRBBI_ONLY(master_lat, master_lat_oppos, lat, crossings) use bmad_interface implicit none type (lat_struct) :: lat type (lat_struct) :: master_lat, master_lat_oppos real(rp), dimension(:,:) :: crossings end subroutine MARK_LRBBI_ONLY end interface type (lat_struct) lat_in, lat_out type (lat_struct), save:: master_lat, master_lat_oppos, lat_1 type (coord_struct), allocatable, save :: co(:) type (ele_struct) ele integer particle, i_train, j_car, n_trains_tot, n_cars integer ierr integer crossings_tot, bunch_tot, total integer, dimension(:), allocatable :: ix_lrbbi, ix_master integer k integer n integer i, j integer ib,it ! Allow adjustment of train and bunch spacing in call to cesr_crossings ! Default spacing integer :: train_spacing(1:10)= (/140, 140, 147,0,0,0,0,0,0,0/) integer :: n_car_spacing(1:10)=(/7,0,0,0,0,0,0,0,0,0/) ! Example of 6, 8 ns spacing ! integer :: n_car_spacing(1:10) /3, 4, 8*0/ ! Example of 2, 12 ns spacing ! integer :: n_car_spacing(1:10) /1, 6, 8*0/ real(rp), dimension(:,:), allocatable :: crossings_1 real(rp), dimension(:), allocatable :: cross_positions ! ADDED TO FIX cesr_crossings integer, dimension(:), allocatable :: ptrain integer, dimension(:), allocatable :: pcar real(rp) current logical rec_taylor type (pc_struct) pc type (constraint_struct) con integer cros1, cros2 call reallocate_coord( co, lat_in%n_ele_max ) ! init lat_in%param%particle = particle call twiss_at_start(lat_in) co(0)%vec = 0. call closed_orbit_calc(lat_in, co, 4) call track_all (lat_in, co) call lat_make_mat6(lat_in,-1,co) call twiss_at_start(lat_in) call twiss_propagate_all (lat_in) ! set up for lrbbi n_trains_tot = con%n_trains * con%n_cars ! total strong beam (positron) bunches n_cars = con%n_trains * con%n_cars ! total weak beam (electron) bunches i_train=1 ! historical reason to keep it if (j_car * n_cars * n_trains_tot * current == 0) then print *,' LRBBI_SETUP: bunch specifications incomplete ' stop endif lat_in%param%n_part = current*0.001 *(lat_in%param%total_length/c_light)/e_charge master_lat = lat_in master_lat_oppos = lat_in lat_1 = lat_in master_lat_oppos%param%n_part = lat_in%param%n_part master_lat_oppos%param%particle = -particle bunch_tot = n_trains_tot ! total strong beam (positron) bunches crossings_tot = bunch_tot * 2 total = crossings_tot allocate(cross_positions(1:crossings_tot), stat=ierr) if(ierr .ne. 0) then print*, "CROSS_POSTIONS: ALLOCATION REQUEST DENIED." call err_exit endif ! ADDED TO FIX cesr_crossings allocate(ptrain(1:crossings_tot), stat=ierr) if(ierr .ne. 0) then print*, "PTRAIN: ALLOCATION REQUEST DENIED." call err_exit endif allocate(pcar(1:crossings_tot), stat=ierr) if(ierr .ne. 0) then print*, "PCAR: ALLOCATION REQUEST DENIED." call err_exit endif allocate(crossings_1(1:total, 1:5), stat=ierr) if(ierr .ne. 0) then print*, "CROSSINGS_1: ALLOCATION REQUEST DENIED." call err_exit endif allocate(ix_lrbbi(1:crossings_tot), stat=ierr) allocate(ix_master(1:crossings_tot), stat=ierr) ! use custom bunch pattern do i=1, bunch_tot !positron bunches cros1=pc%posbunch_elebunch(i,j_car,1) cros2=pc%posbunch_elebunch(i,j_car,2) cross_positions((i-1)*2+1) = pc%cross(cros1)%ele%s /lat_in%param%total_length cross_positions(i*2) = pc%cross(cros2)%ele%s / lat_in%param%total_length enddo do j = 1, crossings_tot crossings_1(j,1) = cross_positions(j) crossings_1(j,2) = 1 crossings_1(j,3) = j crossings_1(j,4) = 0 crossings_1(j,5) = 0 enddo call MARK_LRBBI_only(master_lat, master_lat_oppos, lat_1, crossings_1) call stash_split_lat (i_train, j_car, particle, n_trains_tot, n_cars,master_lat, rec_taylor) call stash_split_lat (i_train, j_car, particle, n_trains_tot, n_cars,master_lat_oppos, rec_taylor) call stash_split_lat (i_train, j_car, particle, n_trains_tot, n_cars,lat_1, rec_taylor) do i=1, crossings_tot find: do k=1,total if(int(crossings_1(k,2)) == 1 .and. int(crossings_1(k,3)) == i)ix_lrbbi(i) = crossings_1(k,4) if(int(crossings_1(k,2)) == 1 .and. int(crossings_1(k,3)) == i)ix_master(i) = crossings_1(k,5) enddo find end do call make_lrbbi(master_lat_oppos, lat_1, ix_lrbbi, ix_master) lat_out = lat_1 deallocate(cross_positions, stat=ierr) deallocate(crossings_1, stat=ierr) deallocate(ix_lrbbi, stat=ierr) deallocate(ix_master, stat=ierr) deallocate (ptrain, stat=ierr) deallocate (pcar, stat=ierr) end subroutine lrbbi_setup_custom