module injlib use bmad use bmadz_mod use bmadz_interface use injparam_mod use z_tune_mod contains !------------------------------------------------------------ ! Series of calls to update the track information ! July 2, 2004 DTK ! ! Modified to recalculate closed orbit after ! redefining lrbbi kicks in beambeam_separation ! Zero n_part if called for positrons ! 26 July 2004 JAC ! ! Added difference orbit diff_co as argument ! 16 Dec 2005 JAC subroutine tracking_update (ring, orbit, i_dim, diff_co, err_flag) implicit none type ( lat_struct ) ring type (coord_struct), allocatable :: orbit(:) type (coord_struct), allocatable :: diff_co(:) type (coord_struct) delta_ip integer i_dim, i_err, status logical, optional :: err_flag logical stable, positron, error, err real(rdef) charge !------------------------------------------------------- error = .false. ! print *,' Rtn TRACKING_UPDATE called for particle print ',ring%param%particle if(ring%param%particle == positron$) then positron = .true. else positron = .false. endif ! call reallocate_coord(orbit,ring%n_ele_max) !allocate(orbit(ring%n_ele_max), STAT=i_err) ! orbit(0)%vec(:) = 0. ! call twiss_at_start( ring ) !call clear_lat_1turn_mats (ring) ! Added 7/9/04 by DTK ! call closed_orbit_calc (ring, orbit, i_dim ) ! call lat_make_mat6( ring, -1, orbit ) ring%param%particle = electron$ ! Always do beambeam_separation with electrons call beambeam_separation ( ring, delta_ip, i_dim, diff_co ) ! Added 7/2/04 by DTK to recalculate and store offsets ! Restore particle type as before call to BEAMBEAM_SEPARATION if(positron) then ring%param%particle = positron$ ! Zero BBI if calculating for positrons charge = ring%param%n_part ring%param%n_part = 0. endif ! Recalculate closed orbit after defining BBI offsets and zero-ing charge if for positrons call init_coord (orbit(0), ele = ring%ele(0), element_end = upstream_end$) ! Call closed_orbit_calc here again even though we just ! did it in beambeam_separation, so that we can test ! the status flag. call closed_orbit_calc (ring, orbit, i_dim, err_flag = error) call lat_make_mat6( ring, -1, orbit ) call twiss_at_start(ring, status) if (status /= ok$) error = .true. ! To ensure errors are caught even if later routines work call twiss_propagate_all(ring, err_flag = err) if (err) error = .true. ! Restore strong beam current if(positron) ring%param%n_part = charge if (present(err_flag)) err_flag = error end subroutine tracking_update !####################################################### subroutine lrbbi_ele_setup (ring, lrbbi_ip, lrbbi_nonip) ! Turn on and off LRBBI ring elements ! 3 June 2004 JAC and DTK implicit none type (coord_struct), allocatable, save :: orb0(:) type (lat_struct) ring logical lrbbi_ip, lrbbi_nonip, ip integer i, nbbi_on, nbbi_off ! allocate storage call reallocate_coord (orb0, ring%n_ele_max) ! Make sure we have a reasonable orbit for remaking the matrices call twiss_at_start(ring) call init_coord (orb0(0), ele = ring%ele(0), element_end = upstream_end$) call closed_orbit_calc(ring, orb0, 4) nbbi_on = 0 nbbi_off = 0 do i=1, ring%n_ele_max if(ring%ele(i)%key /= beambeam$) cycle ip = ring%ele(i)%name == 'IP_COLLISION' !.or. ring%ele(i)%name == 'IP_L0' .or. ring%ele(i)%name == 'IP_L0_END' if((ip .and. lrbbi_ip) .or. (.not. ip .and. lrbbi_nonip)) then ring%ele(i)%is_on = .true. nbbi_on=nbbi_on+1 endif if((ip .and. .not. lrbbi_ip) .or. (.not. ip .and. .not. lrbbi_nonip)) then ring%ele(i)%is_on = .false. nbbi_off=nbbi_off+1 endif call lat_make_mat6(ring, i, orb0) enddo ! orb0%vec(6)=0. ! call twiss_and_track ( ring, orb0 ) write(6,1000)nbbi_on, nbbi_off 1000 format (' Rtn LRBBI_ELE_SETUP has turned on ' i4, ' ring elements.'/, & ' Rtn LRBBI_ELE_SETUP has turned off ' i4, ' ring elements.') ! write(6,1000)nbbi_off !1010 format (' Rtn LRBBI_ELE_SETUP has turned off ' i4, ' ring elements.') end subroutine lrbbi_ele_setup !####################################################### ! Turn om LRBBI ring elements ! 18 May 2004 JAC subroutine lrbbi_ele_on (ring) implicit none type (coord_struct), allocatable, save :: orb0(:) ! type (lat_struct) ring integer i, nbbi ! allocate storage call reallocate_coord (orb0, ring%n_ele_max) nbbi=0 do i=1,ring%n_ele_max if( ring%ele(i)%key == beambeam$ ) then ring%ele(i)%is_on = .true. nbbi=nbbi+1 endif call lat_make_mat6( ring, i, orb0 ) enddo write(6,1000)nbbi 1000 format (' Rtn LRBBI_ON has turned on ',i4,' ring elements.') end subroutine lrbbi_ele_on !####################################################### ! Turn off LRBBI ring elements ! 18 May 2004 JAC subroutine lrbbi_ele_off (ring) implicit none type (coord_struct), allocatable, save :: orb0(:) ! type (lat_struct) ring integer i, nbbi ! allocate storage call reallocate_coord (orb0, ring%n_ele_max) nbbi=0 do i=1,ring%n_ele_max if( ring%ele(i)%key == beambeam$ ) then ring%ele(i)%is_on = .false. nbbi=nbbi+1 endif call lat_make_mat6( ring, i, orb0 ) enddo write(6,1000)nbbi 1000 format (' Rtn LRBBI_OFF has turned off ',i4,' ring elements.') end subroutine lrbbi_ele_off !####################################################### subroutine toggle_lrbbi( ring, is_off, n_lr, ix_lr ) implicit none ! include '[cesr.bmad.code]bmad.inc' ! type (lat_struct) ring type (coord_struct), allocatable, save :: orb0(:) integer i, ix, n_lr, ix_lr(110) logical is_off ! allocate storage call reallocate_coord (orb0, ring%n_ele_max) do i=1,n_lr ix = ix_lr(i) if( is_off ) then ! turn off ring%ele(ix)%is_on = .false. else ring%ele(ix)%is_on = .true. endif call lat_make_mat6( ring, ix, orb0 ) enddo end subroutine toggle_lrbbi !####################################################### subroutine init_chromaticity( ring, xqune ) implicit none type (lat_struct) ring type (xqune_struct) xqune type (coord_struct), allocatable, save :: coord0(:) integer i real(rdef) xix0, xiy0, xix1, xiy1, xix2, xiy2 real dx, dy real(rdef) de real xihcmd, xivcmd ! allocate storage call reallocate_coord (coord0, ring%n_ele_max) ! Zero coord0 do i=0,ring%n_ele_max coord0(i)%vec(:) = 0. enddo ! find Xquneing 1,2 xqune%ix_xiv = 0 xqune%ix_xih = 0 do i=1,ring%n_ele_max ! print *, ' RingName=',ring%ele(i)%name if(ring%ele(i)%name == 'RAW_XQUNE1' .or. & ring%ele(i)%name == 'RAW_XQUNE_1' ) then ! Vertical Chromaticity xqune%ix_xiv = i endif if(ring%ele(i)%name == 'RAW_XQUNE2' .or. & ring%ele(i)%name == 'RAW_XQUNE_2' ) then ! Horizontal Chromaticity xqune%ix_xih = i endif enddo if( xqune%ix_xiv == 0 .or. xqune%ix_xih == 0 ) then print *,' Cant find Xquening 1 or 2' stop endif ! Get the starting chromaticity de = 0.0001 call twiss_at_start (ring) call twiss_propagate_all (ring) call envelope_chrom_calc(ring, de, xix0, xiy0 ) print *,' Xih = ',xix0,' Xiv = ',xiy0 ! Calibrate Horizontal Chromaticity knob xihcmd = 4. xivcmd = 0. ring%ele(xqune%ix_xih)%control%var(1)%value = xihcmd ring%ele(xqune%ix_xiv)%control%var(1)%value = xivcmd call lat_make_mat6(ring, xqune%ix_xih, coord0) call lat_make_mat6(ring, xqune%ix_xiv, coord0) call twiss_at_start (ring) call twiss_propagate_all (ring) call envelope_chrom_calc(ring, de, xix1, xiy1 ) dx = xix1-xix0 dy = xiy1-xiy0 xqune%xih_coef = dx/xihcmd print *,' H Chromaticity Measurement: dXix = ',dx, & ' dXiy = ',dy,' coef = ',xqune%xih_coef if( abs(dy) > 0.1 * abs(dx) ) then print *,' *** WARNING ***' print *,' *** Significant Vertical Chromaticity knob coupling !!' print *,' ***************' endif ! Calibrate Vertical Chromaticity knob xihcmd = 0. xivcmd = 4. ring%ele(xqune%ix_xih)%control%var(1)%value = xihcmd ring%ele(xqune%ix_xiv)%control%var(1)%value = xivcmd call lat_make_mat6(ring, xqune%ix_xih, coord0) call lat_make_mat6(ring, xqune%ix_xiv, coord0) call twiss_at_start (ring) call twiss_propagate_all (ring) call envelope_chrom_calc(ring, de, xix2, xiy2 ) dx = xix2-xix0 dy = xiy2-xiy0 xqune%xiv_coef = dy/xivcmd print *,' V Chromaticity Measurement: dXiy = ',dy, & ' dXix = ',dx,' coef = ',xqune%xiv_coef if( abs(dx) > 0.1 * abs(dy) ) then print *,' *** WARNING ***' print *,' *** Significant Horizontal Chromaticity knob coupling !!' print *,' ***************' endif ! Now reset to design ring%ele(xqune%ix_xih)%control%var(1)%value = 0. ring%ele(xqune%ix_xiv)%control%var(1)%value = 0. call lat_make_mat6(ring, xqune%ix_xih, coord0) call lat_make_mat6(ring, xqune%ix_xiv, coord0) call twiss_at_start (ring) call twiss_propagate_all (ring) call envelope_chrom_calc(ring, de, xix2, xiy2 ) print *,' Reset Chromatiticites: ',xix2, xiy2 end subroutine init_chromaticity !####################################################### subroutine init_tonality( ring, xqune, tunes, pretzel ) implicit none type (lat_struct) ring type (xqune_struct) xqune type (coord_struct), allocatable, save :: orb0(:) type (coord_struct), allocatable, save :: orbit(:) type (pretzel_struct) pretzel type (tunes_struct) tunes integer i, j real xix0, xiy0, xix1, xiy1, xix2, xiy2 real qx0, qx1, qx2, qy0, qy1, qy2 real dx, dy real de, thcmd, tvcmd, xqcmd ! allocate storage call reallocate_coord (orbit, ring%n_ele_max) call reallocate_coord (orb0, ring%n_ele_max) ! find Xquneing 3,4 xqune%ix_tonh = 0 xqune%ix_tonv = 0 do i=1,ring%n_ele_max if(ring%ele(i)%name == 'RAW_XQUNE3' .or. & ring%ele(i)%name == 'RAW_XQUNE_3' ) then ! Vertical Tonality xqune%ix_tonv = i endif if(ring%ele(i)%name == 'RAW_XQUNE4' .or. & ring%ele(i)%name == 'RAW_XQUNE_4' ) then ! Horizontal Tonality xqune%ix_tonh = i endif enddo if( xqune%ix_tonv == 0 .or. xqune%ix_tonh == 0 ) then print *,' Cant find Xquening 3 or 4' stop endif ! print *,' Enter XQ command to used for tonality calibration' ! print *,' 10. for 4s, 1. for 3s or CESRc' ! read *, xqcmd ! for 4s lattice ! xqcmd = 10. ! for 3s lattice xqcmd = 1.0 print *,' XQ command value for tonality calibration: ',xqcmd ring%param%particle = positron$ ! Turn BBI elements off for positrons call lrbbi_ele_setup ( ring, .false., .false. ) do i = 1,2 ! do both species call twiss_and_track( ring, orbit ) ! call envelope_chrom_calc(ring, de, xix0, xiy0 ) print *,' orb 1 : orbit= ',(orbit(0)%vec(j),j=1,4) qx0 = ring%a%tune/(2.*pi) qy0 = ring%b%tune/(2.*pi) print *,' qx0 = ',qx0,' qy0= ',qy0 if( ring%param%particle == electron$ ) then tunes%q0h_ele = qx0 tunes%q0v_ele = qy0 else tunes%q0h_pos = qx0 tunes%q0v_pos = qy0 endif thcmd = xqcmd tvcmd = 0. ring%ele(xqune%ix_tonh)%control%var(1)%value = thcmd ring%ele(xqune%ix_tonv)%control%var(1)%value = tvcmd call lat_make_mat6( ring, xqune%ix_tonh, orbit ) call lat_make_mat6( ring, xqune%ix_tonv, orbit ) print *,' Calling TWISS_AT_START' call twiss_at_start(ring) print *,' Calling TWISS_PROPAGATE_ALL' call twiss_propagate_all(ring) ! call envelope_chrom_calc(ring, de, xix1, xiy1 ) qx1 = ring%a%tune/(2.*pi) qy1 = ring%b%tune/(2.*pi) dx = qx1 - qx0 dy = qy1 - qy0 if( ring%param%particle == electron$ ) then xqune%ele_tonh_coef = dx/thcmd print *,' e- H Tonality Measurement: dQx = ',dx, & ' dQy = ',dy,' coef = ',xqune%ele_tonh_coef else xqune%pos_tonh_coef = dx/thcmd print *,' e+ H Tonality Measurement: dQx = ',dx, & ' dQy = ',dy,' coef = ',xqune%pos_tonh_coef endif if( abs(dy) > 0.1 * abs(dx) ) then print *,' *** WARNING ***' print *,' *** Significant (>10%) Vertical Tonality knob coupling !!' print *,' ***************' endif ! Now do the Vertical Tonality knob thcmd = 0. tvcmd = -xqcmd ring%ele(xqune%ix_tonh)%control%var(1)%value = thcmd ring%ele(xqune%ix_tonv)%control%var(1)%value = tvcmd call lat_make_mat6( ring, xqune%ix_tonh, orbit ) call lat_make_mat6( ring, xqune%ix_tonv, orbit ) call twiss_at_start(ring) call twiss_propagate_all(ring) ! call envelope_chrom_calc(ring, de, xix1, xiy1 ) qx2 = ring%a%tune/(2.*pi) qy2 = ring%b%tune/(2.*pi) dx = qx2 - qx0 dy = qy2 - qy0 if( ring%param%particle == electron$ ) then xqune%ele_tonv_coef = dy/tvcmd print *,' e- V Tonality Measurement: dQy = ',dy, & ' dQx = ',dx,' coef = ',xqune%ele_tonv_coef else xqune%pos_tonv_coef = dy/tvcmd print *,' e+ V Tonality Measurement: dQy = ',dy, & ' dQx = ',dx,' coef = ',xqune%pos_tonv_coef endif if( abs(dx) > 0.1 * abs(dy) ) then print *,' *** WARNING ***' print *,' *** Significant (>10%) Horizontal Tonality knob coupling !!' print *,' ***************' endif ! Now Reset thcmd = 0. tvcmd = 0. ring%ele(xqune%ix_tonh)%control%var(1)%value = thcmd ring%ele(xqune%ix_tonv)%control%var(1)%value = tvcmd call lat_make_mat6( ring, xqune%ix_tonh, orbit ) call lat_make_mat6( ring, xqune%ix_tonv, orbit ) call twiss_at_start(ring) call twiss_propagate_all(ring) if( i == 1 ) then ring%param%particle = electron$ ! Turn BBI elements on for electrons call lrbbi_ele_setup ( ring, .true., .true. ) else ring%param%particle = positron$ ! Turn BBI elements off for positrons call lrbbi_ele_setup ( ring, .false., .false. ) endif enddo end subroutine init_tonality !####################################################### subroutine init_qtuneing( ring, qtune ) implicit none type (lat_struct) ring type (qtune_struct) qtune type (coord_struct), allocatable, save :: coord0(:) integer i real qx0, qx1, qx2, qy0, qy1, qy2 real qtxcmd, qtycmd ! allocate storage call reallocate_coord (coord0, ring%n_ele_max) print *,'Entering Rtn INIT_QTUNEING ...' ! Zero coord0 do i=0,ring%n_ele_max coord0(i)%vec(:) = 0. enddo ! find qtune5, 6 qtune%ix_qtv = 0 qtune%ix_qth = 0 do i=1,ring%n_ele_max if(ring%ele(i)%name == 'RAW_QTUNE5' .or. & ring%ele(i)%name == 'RAW_QTUNE_5' ) then ! Vertical Tune qtune%ix_qtv = i endif if(ring%ele(i)%name == 'RAW_QTUNE6' .or. & ring%ele(i)%name == 'RAW_QTUNE_6' ) then ! Horizontal Tune qtune%ix_qth = i endif enddo if( qtune%ix_qtv == 0 .or. qtune%ix_qth == 0 ) then print *,' Cant find qtune 5 or 6' stop endif ! Get the tunes (the pretzel is now off) call twiss_at_start (ring) call twiss_propagate_all (ring) qx0 = ring%a%tune/(2.*pi) qy0 = ring%b%tune/(2.*pi) ! Now adjust Qh and measure the tunes qtxcmd = 0.02 ring%ele(qtune%ix_qth)%control%var(1)%value = qtxcmd call lat_make_mat6(ring, qtune%ix_qth, coord0) call twiss_at_start (ring) call twiss_propagate_all (ring) qx1 = ring%a%tune/(2.*pi) qy1 = ring%b%tune/(2.*pi) qtune%hcoef = (qx1-qx0)/qtxcmd print *,' Htune Measurement: dQx = ',qx1-qx0,' dQy = ',qy1-qy0, & ' coef = ',qtune%hcoef if( abs(qy1-qy0) > 0.1 * abs(qx1-qx0) ) then print *,' *** WARNING ***' print *,' *** Significant Vertical Tune knob coupling !!' print *,' ***************' endif ! Now adjust Qv and measure the tunes qtycmd = 0.02 ring%ele(qtune%ix_qth)%control%var(1)%value = 0.0 ring%ele(qtune%ix_qtv)%control%var(1)%value = qtycmd call lat_make_mat6(ring, qtune%ix_qth, coord0) call lat_make_mat6(ring, qtune%ix_qtv, coord0) call twiss_at_start (ring) call twiss_propagate_all (ring) qx2 = ring%a%tune/(2.*pi) qy2 = ring%b%tune/(2.*pi) qtune%vcoef = (qy2-qy0)/qtycmd print *,' Vtune Measurement: dQy = ',qy2-qy0,' dQx = ',qx2-qx0, & ' coef = ',qtune%vcoef if( abs(qx2-qx0) > 0.1 * abs(qy2-qx0) ) then print *,' *** WARNING ***' print *,' *** Significant Horizontal Tune knob coupling !!' print *,' ***************' endif ! Now reset to design tunes ring%ele(qtune%ix_qth)%control%var(1)%value = 0.0 ring%ele(qtune%ix_qtv)%control%var(1)%value = 0.0 call lat_make_mat6(ring, qtune%ix_qth, coord0) call lat_make_mat6(ring, qtune%ix_qtv, coord0) call twiss_at_start (ring) call twiss_propagate_all (ring) print *,' Finished Qtune knob setup:' print *,' Final tunes: qx = ', ring%a%tune/(2.*pi) print *,' Final tunes: qy = ', ring%b%tune/(2.*pi) print *,'Exiting Rtn INIT_QTUNEING' end subroutine init_qtuneing !####################################################### subroutine set_tunes( ring, pretzel, qtune, xqune, nbeam, & use_quads, use_sexts, species, qhtarg, qvtarg, tol ) implicit none type (lat_struct) ring type (pretzel_struct) pretzel type (qtune_struct) qtune type (xqune_struct) xqune type (coord_struct), allocatable, save :: orb1_(:), orb2_(:) integer nbeam, iter integer species(2), ie, ip real tol real qhtarg(2), qvtarg(2) real qx0(2), qy0(2), qx1(2), qy1(2), cmdh, cmdv, dqthcmd, dqtvcmd real dqhp, dqhe, dqvp, dqve, dxqhcmd, dxqvcmd logical use_quads, use_sexts, converged ! allocate storage call reallocate_coord (orb1_, ring%n_ele_max) call reallocate_coord (orb2_, ring%n_ele_max) converged = .false. iter = 0 ! Get orbit and tunes for species 1 ring%param%particle = species(1) orb1_(0)%vec(:) = 0. ! call closed_orbit_at_start(ring, orb1_(0),4,.true.) ! call track_all(ring, orb1_) call closed_orbit_calc (ring, orb1_,4 ) call lat_make_mat6(ring, -1, orb1_) call twiss_at_start(ring) call twiss_propagate_all(ring) qx0(1) = ring%a%tune/(2.*pi) qy0(1) = ring%b%tune/(2.*pi) print *,' Target tunes for 1st beam: ',qhtarg(1),qvtarg(1) print *,' Initial tunes for 1st beam: ',qx0(1),qy0(1) if( nbeam == 2 ) then ! get the tune of the other beam ring%param%particle = species(2) orb2_(0)%vec(:) = 0. ! call closed_orbit_at_start(ring, orb2_(0),4,.true.) ! call track_all(ring, orb2_) call closed_orbit_calc(ring, orb2_,4 ) call lat_make_mat6(ring, -1, orb2_) call twiss_at_start(ring) call twiss_propagate_all(ring) qx0(2) = ring%a%tune/(2.*pi) qy0(2) = ring%b%tune/(2.*pi) ring%param%particle = species(1) print *,' Target tunes for 2nd beam: ',qhtarg(2),qvtarg(2) print *,' Initial tunes for 2nd beam: ',qx0(2),qy0(2) endif print *,' Routine SET_TUNES: use_quads, use_sexts=',use_quads,use_sexts do while( .not. converged ) iter = iter + 1 if( use_quads .and. .not. use_sexts ) then ! Only Qtuneing if( nbeam /= 1 ) then print *,' cant adjust tunes of both beams only with quads' stop else ! single beam, just qtune dqthcmd = ( qhtarg(1) - qx0(1) )/qtune%hcoef dqtvcmd = ( qvtarg(1) - qy0(1) )/qtune%vcoef endif cmdh = ring%ele(qtune%ix_qth)%control%var(1)%value ring%ele(qtune%ix_qth)%control%var(1)%value = cmdh + dqthcmd cmdv = ring%ele(qtune%ix_qtv)%control%var(1)%value ring%ele(qtune%ix_qtv)%control%var(1)%value = cmdv + dqtvcmd call lat_make_mat6(ring, qtune%ix_qth, orb1_) call lat_make_mat6(ring, qtune%ix_qtv, orb1_) call twiss_at_start (ring) call twiss_propagate_all (ring) qx1(1) = ring%a%tune/(2.*pi) qy1(1) = ring%b%tune/(2.*pi) converged = abs(qhtarg(1)-qx1(1)) <= tol & .and. abs(qvtarg(1)-qy1(1)) <= tol print *,' cmdh = ',cmdh,' cmdv= ',cmdv if( converged ) then print *,' qtuneing converged in ', iter,'iters: qx, qy= ',qx1(1),qy1(1) else qx0(1) = qx1(1) qy0(1) = qy1(1) endif else if( .not. use_quads .and. use_sexts ) then ! Only Xquneing if( nbeam /= 1 ) then print *,' cant adjust tunes of both beams only with sextupoles' stop endif print *,' cant do this yet' stop else ! use both if( nbeam /= 2 ) then print *,' Dont know how to adjust one tune with quad and sexts' stop else if( species(1) == electron$ ) then ie = 1 ip = 2 else ie = 2 ip = 1 endif dqhe = qhtarg(ie) - qx0(ie) dqhp = qhtarg(ip) - qx0(ip) dqve = qvtarg(ie) - qy0(ie) dqvp = qvtarg(ip) - qy0(ip) dxqhcmd = ( dqhp - dqhe )/( xqune%pos_tonh_coef - xqune%ele_tonh_coef ) dxqvcmd = ( dqvp - dqve )/( xqune%pos_tonv_coef - xqune%ele_tonv_coef ) dqthcmd = ( dqhp - xqune%pos_tonh_coef*dxqhcmd )/qtune%hcoef dqtvcmd = ( dqvp - xqune%pos_tonv_coef*dxqvcmd )/qtune%vcoef print *,'Tunes this iteration: qxe=',qx0(ie),' qxp=',qx0(ip), & ' qye=',qy0(ie),' qyp=',qy0(ip) print *, 'Desired tonality change: dqh=',dqhp-dqhe,' dqv=',dqvp-dqve print *,'Increment cmd values: dxqh= ',dxqhcmd,' dxqv= ',dxqvcmd print *,' dqth= ',dqthcmd,' dqtv= ',dqtvcmd ! Load qtuneing knobs cmdh = ring%ele(qtune%ix_qth)%control%var(1)%value ring%ele(qtune%ix_qth)%control%var(1)%value = cmdh + dqthcmd cmdv = ring%ele(qtune%ix_qtv)%control%var(1)%value ring%ele(qtune%ix_qtv)%control%var(1)%value = cmdv + dqtvcmd call lat_make_mat6(ring, qtune%ix_qth, orb1_) call lat_make_mat6(ring, qtune%ix_qtv, orb1_) ! Load xquneing knobs cmdh = ring%ele(xqune%ix_tonh)%control%var(1)%value ring%ele(xqune%ix_tonh)%control%var(1)%value = cmdh + dxqhcmd cmdv = ring%ele(xqune%ix_tonv)%control%var(1)%value ring%ele(xqune%ix_tonv)%control%var(1)%value = cmdv + dxqvcmd call lat_make_mat6(ring, xqune%ix_tonh, orb1_) call lat_make_mat6(ring, xqune%ix_tonv, orb1_) call lat_make_mat6(ring, -1, orb1_) ! redo everything for sp1 orbit call twiss_at_start (ring) call twiss_propagate_all (ring) qx1(1) = ring%a%tune/(2.*pi) qy1(1) = ring%b%tune/(2.*pi) ! Get tune of the other beam (species 2) call lat_make_mat6(ring, -1, orb2_) call twiss_at_start(ring) call twiss_propagate_all(ring) qx1(2) = ring%a%tune/(2.*pi) qy1(2) = ring%b%tune/(2.*pi) converged = abs(qhtarg(1)-qx1(1)) <= tol & .and. abs(qvtarg(1)-qy1(1)) <= tol & .and. abs(qhtarg(2)-qx1(2)) <= tol & .and. abs(qvtarg(2)-qy1(2)) <= tol if( converged ) then print *,' qtuneing converged in ', iter,'iters:' print *,' Beam1: qx, qy= ',qx1(1),qy1(1) print *,' Beam2: qx, qy= ',qx1(2),qy1(2) ! Reset ring for species 1 ring%param%particle = species(1) orb1_(0)%vec(:) = 0. ! call closed_orbit_at_start(ring, orb1_(0),4,.true.) ! call track_all(ring, orb1_) call closed_orbit_calc(ring, orb1_,4 ) call lat_make_mat6(ring, -1, orb1_) call twiss_at_start(ring) call twiss_propagate_all(ring) else qx0(1) = qx1(1) qy0(1) = qy1(1) qx0(2) = qx1(2) qy0(2) = qy1(2) endif endif endif enddo end subroutine set_tunes !####################################################### subroutine set_qtuneing( ring, qtune, qhtarg, qvtarg, tol ) implicit none type (lat_struct) ring type (qtune_struct) qtune type (coord_struct), allocatable, save :: orb1_(:) integer iter integer ie, ip real tol real qhtarg, qvtarg real qx0, qy0, qx1, qy1, cmdh, cmdv, dqthcmd, dqtvcmd logical converged ! allocate storage call reallocate_coord (orb1_, ring%n_ele_max) converged = .false. iter = 0 ! Get orbit and tunes for species 1 orb1_(0)%vec(:) = 0. ! call closed_orbit_at_start(ring, orb1_(0),4,.true.) ! call track_all(ring, orb1_) call closed_orbit_calc(ring, orb1_,4 ) call lat_make_mat6(ring, -1, orb1_) call twiss_at_start(ring) call twiss_propagate_all(ring) qx0 = ring%a%tune/(2.*pi) qy0 = ring%b%tune/(2.*pi) print *,' Initial tunes: ',qx0,qy0 do while( .not. converged ) iter = iter + 1 dqthcmd = ( qhtarg - qx0 )/qtune%hcoef dqtvcmd = ( qvtarg - qy0 )/qtune%vcoef cmdh = ring%ele(qtune%ix_qth)%control%var(1)%value ring%ele(qtune%ix_qth)%control%var(1)%value = cmdh + dqthcmd cmdv = ring%ele(qtune%ix_qtv)%control%var(1)%value ring%ele(qtune%ix_qtv)%control%var(1)%value = cmdv + dqtvcmd call lat_make_mat6(ring, qtune%ix_qth, orb1_) call lat_make_mat6(ring, qtune%ix_qtv, orb1_) call twiss_at_start (ring) call twiss_propagate_all (ring) qx1 = ring%a%tune/(2.*pi) qy1 = ring%b%tune/(2.*pi) converged = abs(qhtarg-qx1) <= tol & .and. abs(qvtarg-qy1) <= tol if( converged ) then print *,' qtuneing converged in ', iter,'iters: qx, qy= ',qx1,qy1 else qx0 = qx1 qy0 = qy1 endif enddo end subroutine set_qtuneing !####################################################### subroutine set_pretzel( state, ring, pretzel ) implicit none type (lat_struct) ring type (pretzel_struct) pretzel type (coord_struct), allocatable, save :: orbit(:), coord0(:) type (coord_struct), allocatable :: diff_co(:) integer state, i real kick ! allocate storage call reallocate_coord (orbit, ring%n_ele_max) call reallocate_coord (coord0, ring%n_ele_max) call reallocate_coord (diff_co, ring%n_ele_max) ! write (6,1000)state,(i,ring%ele( pretzel%ix_hsep(i) )%value(hkick$),i=1,4) 1000 format(' Entering Rtn SET_PRETZEL with state',i2/, & 4(' Kick',i2,' is ',e13.6/)) if( state == 0 ) then ! Turn off separators do i=1,4 ring%ele( pretzel%ix_hsep(i) )%value(hkick$) = 0.0 ! call lat_make_mat6(ring, pretzel%ix_hsep(i), coord0) enddo else if( state == 1 ) then ! Put pretzel at saved value do i=1,4 ring%ele( pretzel%ix_hsep(i) )%value(hkick$) = & pretzel%hkick(i) ! call lat_make_mat6(ring, pretzel%ix_hsep(i), coord0) enddo else if( state == -1 ) then ! Flip separator signs do i=1,4 kick = ring%ele( pretzel%ix_hsep(i) )%value(hkick$) ring%ele( pretzel%ix_hsep(i) )%value(hkick$) = -kick ! call lat_make_mat6(ring, pretzel%ix_hsep(i), coord0) enddo endif ! Recalculate closed orbit, transfer matrices and twiss functions ! orbit%vec(6)=0. ! call twiss_and_track( ring, orbit ) ! New simplified routine added 7/2/04 by DTK ! This updates all the transfer matrices including bbi call tracking_update(ring, orbit, 4, diff_co) write (6,5000)state,(i,ring%ele( pretzel%ix_hsep(i) )%value(hkick$),i=1,4) 5000 format(' Exiting Rtn SET_PRETZEL called with state =',i4,/4(' Kick',i2,' is ',e13.6/)) end subroutine set_pretzel !####################################################### subroutine init_pretzel( ring, pretzel, prz_orbit ) implicit none type (lat_struct), target :: ring type (ele_struct), pointer :: ele type (coord_struct) prz_orbit(0:ring%n_ele_max) type (pretzel_struct) pretzel integer nhs, i, ihs ! Get element indices of separators ! If a separator has been split by a beam-beam element, ! get the index of the superlord which controls the slave ! elements. ! since the split separator elements have control type ! super_slave. 6/22/04 jac nhs = 0 do i=1,ring%n_ele_max ele => ring%ele(i) if (ele%key == elseparator$ .and. (ele%slave_status == free$ .or. & ele%slave_status == minor_slave$ .or. ele%lord_status == super_lord$)) then ! if( ele%name(1:1) == 'H' ) then ! print *,' Rtn: INIT_PRETZEL: Found separator ',ele%name,' index= ',i,' Control type',ele%control_type ! print *,' Kick: ',ele%value(hkick$) ihs = 0 if ( ele%name=='H_SEP_08W' .or. & ele%name=='H_SEP_09W') then ihs = 1 nhs = nhs + 1 elseif ( ele%name=='H_SEP_45W') then ihs = 2 nhs = nhs + 1 elseif ( ele%name=='H_SEP_45E') then ihs = 3 nhs = nhs + 1 elseif ( ele%name=='H_SEP_08E') then ihs = 4 nhs = nhs + 1 endif ! print *,' Storing index', i,' for sep nr',ihs if(ihs.ge.1)then pretzel%ix_hsep(ihs) = i pretzel%hkick_theory(ihs) = ele%value(hkick$) pretzel%hkick(ihs) = pretzel%hkick_theory(ihs) endif endif endif enddo ! Insist on finding exactly four horizontal separators if( nhs /= 4 ) then print *,' Did not find exactly FOUR horizontal separators!' print *,' Found ',nhs,' separators' if (nhs.gt.0)then write (6,5000)nhs,(pretzel%ix_hsep(i), i, & ring%ele( pretzel%ix_hsep(i) )%name, & ring%ele( pretzel%ix_hsep(i) )%value(hkick$),i=1,nhs) endif ! stop endif write (6,5000)nhs,(pretzel%ix_hsep(i), i, & ring%ele( pretzel%ix_hsep(i) )%name, & ring%ele( pretzel%ix_hsep(i) )%value(hkick$),i=1,nhs) 5000 format(' Rtn INIT_PRETZEL finds',i2, & ' horizontal separators:', & /4(' Element index',i5,': Kick',i2,& '(',a10,') is ',e13.6/)) write(6,6000) 1000* prz_orbit(0)%vec(1), 1000*prz_orbit(0)%vec(2) 6000 format(' Rtn INIT_PRETZEL - At the IP, e+ position = ',f6.2, & ' mm angle = ',f6.3,' mrad') ! Horizontal angle at IP pretzel%prz1_theory = prz_orbit(0)%vec(2) end subroutine init_pretzel !####################################################### ! If PRZ1 nonzero, use both PRZ1 and PRZ13 values to calculate sep kicks! If PRZ1 zero and PRZ13 nonzero, modify sep kicks to give IP separation! If both PRZ1 and PRZ13 zero, do nothing (keep default pretzel) ! Note that successive calls have an additive effect on the kicks, ! so calling once with prz13=1000 and again with prz13=-1000 will ! result in the kicks of 8E and 8W unchanged. ! 26 July 2004 JAC subroutine make_inj_pretzel( ring, pretzel, prz1, prz13 ) implicit none type (lat_struct) ring type (pretzel_struct) pretzel type (coord_struct) delta_ip type (coord_struct), allocatable, save :: orb0(:) type (coord_struct), allocatable :: diff_co(:) real prz1, prz13 integer i, ix real kick ! allocate storage call reallocate_coord (orb0, ring%n_ele_max) call reallocate_coord (diff_co, ring%n_ele_max) print *, 'Rtn MAKE_INJ_PRETZEL called with prz1, prz13=',prz1,prz13 if ( prz1 .ne. 0 ) then ! Now scale the separator strengths (seps 1, 2, 5, 6 ) ! Separator 8W (smaller than 8E for negative PRZ13) do i=1,4 ix = pretzel%ix_hsep(i) if( i == 1 ) then ! 8W kick = (prz1 + prz13)*pretzel%hkick_theory(i)/(1.e6*pretzel%prz1_theory) else if( i == 2 .or. i == 3 ) then !45W/E ! Scale the prz1 setting with the horizontal angle at the IP kick = prz1*pretzel%hkick_theory(i)/(1.e6*pretzel%prz1_theory) else if( i == 4 ) then ! 8E kick = (prz1 - prz13)*pretzel%hkick_theory(i)/(1.e6*pretzel%prz1_theory) endif ring%ele( ix )%value(hkick$) = kick pretzel%hkick(i) = kick print *,ring%ele(ix)%name,' Rtn MAKE_INJ_PRETZEL: theory sep kick = ',pretzel%hkick_theory(i), & ' Sep Kick = ',pretzel%hkick(i) call type_ele(ring%ele(ix),.false.,0,.false., radians$, .true.) enddo ! New simplified routine added 7/2/04 by DTK call tracking_update(ring, orb0, 4, diff_co) elseif ( prz13 .ne. 0 ) then ! If prz1 is zero and prz13 is not, then modify only 8E/W kicks do i=1,4,3 ix = pretzel%ix_hsep(i) kick = ring%ele( ix )%value(hkick$) if( i == 1 ) then ! 8W kick = kick + prz13*pretzel%hkick_theory(i)/(1.e6*pretzel%prz1_theory) else if( i == 4 ) then ! 8E kick = kick - prz13*pretzel%hkick_theory(i)/(1.e6*pretzel%prz1_theory) endif ring%ele( ix )%value(hkick$) = kick pretzel%hkick(i) = kick print *,ring%ele(ix)%name,' Rtn MAKE_INJ_PRETZEL: theory sep kick = ',pretzel%hkick_theory(i), & ' Sep Kick = ',pretzel%hkick(i) call type_ele(ring%ele(ix),.false.,0,.false., radians$, .true.) enddo ! New simplified routine added 7/2/04 by DTK call tracking_update(ring, orb0, 4, diff_co) endif end subroutine make_inj_pretzel !####################################################### subroutine make_mag_pretzel( ring, pretzel, wmag, emag ) implicit none type (lat_struct) ring type (pretzel_struct) pretzel type (coord_struct) delta_ip type (coord_struct), allocatable, save :: coord0(:) real wmag, emag integer i, ix, j real kick ! allocate storage call reallocate_coord (coord0, ring%n_ele_max) ! Zero coord0 do i=0,ring%n_ele_max do j=1,6 coord0(i)%vec(j) = 0.0 enddo enddo ! Do west first do i=1,8 ix = pretzel%ix_mag(i) if (i<=4) then ! West kick = wmag*pretzel%mag_hkick_theory(i)/(1.e6*pretzel%prz1_theory) else ! East kick = emag*pretzel%mag_hkick_theory(i)/(1.e6*pretzel%prz1_theory) endif ring%ele( ix )%value(hkick$) = kick pretzel%mag_hkick(i) = kick ! print *,ring%ele(ix)%name,' Magpie Theory = ', & ! pretzel%mag_hkick_theory(i),' Kick = ',pretzel%mag_hkick(i) call lat_make_mat6(ring, ix, coord0) enddo end subroutine make_mag_pretzel !####################################################### subroutine init_mag_pretzel( ring, pretzel, prz_orbit) implicit none type (lat_struct) ring type (coord_struct) prz_orbit(0:ring%n_ele_max) type (coord_struct), allocatable, save :: orbtest_(:), coord0(:) type (pretzel_struct) pretzel integer i, nmp, kk integer im1, im2, im3, im4 real a, b, c, d, k4, k6, k44, k46 ! ! print *,' Entering Rtn INIT_MAG_PRETZEL' ! allocate storage call reallocate_coord (orbtest_, ring%n_ele_max) call reallocate_coord (coord0, ring%n_ele_max) ! Zero coord structures do i=0,ring%n_ele_max coord0(i)%vec(:)=0 orbtest_(i)%vec(:)=0 enddo ! get indices of correctors used for magpies nmp = 0 im2 = 0 im3 = 0 im4 = 0 do i=1,ring%n_ele_max if ( ring%ele(i)%key == quadrupole$ .and. & ring%ele(i)%name == 'Q43W' ) then im2 = i endif if ( ring%ele(i)%key == quadrupole$ .and. & ring%ele(i)%name == 'Q43E' ) then im3 = i endif if ( ring%ele(i)%key == quadrupole$ .and. & ring%ele(i)%name == 'Q09E' ) then im4 = i endif if ( ring%ele(i)%key == overlay$ ) then if( ring%ele(i)%name == 'HB01' ) then pretzel%ix_mag(1) = i nmp = nmp+1 ! call type_ele(ring.ele_(i),.false.,0,.false., 0, .true.) else if( ring%ele(i)%name == 'HB02' ) then pretzel%ix_mag(2) = i nmp = nmp+1 ! call type_ele(ring.ele_(i),.false.,0,.false., 0, .true.) else if( ring%ele(i)%name == 'H44W' .or. & ring%ele(i)%name == 'H44') then pretzel%ix_mag(3) = i nmp = nmp+1 else if( ring%ele(i)%name == 'H46W' .or. & ring%ele(i)%name == 'H46') then pretzel%ix_mag(4) = i nmp = nmp+1 else if( ring%ele(i)%name == 'H46E' .or. & ring%ele(i)%name == 'H53') then pretzel%ix_mag(5) = i nmp = nmp+1 else if( ring%ele(i)%name == 'H44E' .or. & ring%ele(i)%name == 'H55') then pretzel%ix_mag(6) = i nmp = nmp+1 else if( ring%ele(i)%name == 'HB05' ) then pretzel%ix_mag(7) = i nmp = nmp+1 else if( ring%ele(i)%name == 'HB06' ) then pretzel%ix_mag(8) = i nmp = nmp+1 endif endif enddo if( nmp /= 8 ) then print *,' Cant find all correctors for magnetic pretzels' stop endif if( im2 == 0 ) then print *,' Cant find Q43W' stop endif if( im3 == 0 ) then print *,' Cant find Q43E' stop endif if( im4 == 0 ) then print *,' Cant find Q09E' stop endif ! Turn off separators for now do i=1,4 ring%ele( pretzel%ix_hsep(i) )%value(hkick$) = 0.0 enddo ! Now figure out the magpie kicks - West first ring%ele( pretzel%ix_mag(1) )%value(hkick$) = 0.001 call lat_make_mat6(ring, -1, coord0) call init_coord(orbtest_(0), ele = ring%ele(0), element_end = upstream_end$) call track_all (ring, orbtest_) im1 = pretzel%ix_hsep(1) ! match point is HSEP 8W a = orbtest_(im1)%vec(1)/0.001 c = orbtest_(im1)%vec(2)/0.001 ! print *,' orb = ',orbtest_(im1)%vec(1), orbtest_(im1)%vec(2) ring%ele( pretzel%ix_mag(1) )%value(hkick$) = 0.0 ring%ele( pretzel%ix_mag(2) )%value(hkick$) = 0.001 orbtest_(0)%vec(:) = 0. call lat_make_mat6(ring, -1, coord0) call track_all (ring, orbtest_) b = orbtest_(im1)%vec(1)/0.001 d = orbtest_(im1)%vec(2)/0.001 ! print *,' orb = ',orbtest_(im1)%vec(1), orbtest_(im1)%vec(2) k6 = (prz_orbit(im1)%vec(1)/a - prz_orbit(im1)%vec(2)/c)/(b/a - d/c) ! print *,' k6 = ',k6 k4 = prz_orbit(im1)%vec(1)/a - k6*b/a ! print *,' k4 = ',k4 ring%ele( pretzel%ix_mag(1) )%value(hkick$) = k4 ring%ele( pretzel%ix_mag(2) )%value(hkick$) = k6 orbtest_(0)%vec(:) = 0. call lat_make_mat6(ring, -1, coord0) call track_all (ring, orbtest_) ! print *,' ***** At match point: X = ',orbtest_(im1)%vec(1), & ! ' Desired = ',prz_orbit(im1)%vec(1) ! print *,' Xp = ',orbtest_(im1)%vec(2), & ! ' Desired = ',prz_orbit(im1)%vec(2) ring%ele( pretzel%ix_mag(1) )%value(hkick$) = 0.0 ring%ele( pretzel%ix_mag(2) )%value(hkick$) = 0.0 ring%ele( pretzel%ix_mag(4) )%value(hkick$) = 0.001 orbtest_(0)%vec(:) = 0. call lat_make_mat6(ring, -1, coord0) ring%param%particle = antiparticle(ring%param%particle) call track_many (ring, orbtest_, 0, 0, -1) ring%param%particle = antiparticle(ring%param%particle) a = orbtest_(im2)%vec(1)/0.001 c = orbtest_(im2)%vec(2)/0.001 ring%ele( pretzel%ix_mag(4) )%value(hkick$) = 0.0 ring%ele( pretzel%ix_mag(3) )%value(hkick$) = 0.001 orbtest_(0)%vec(:) = 0. call lat_make_mat6(ring, -1, coord0) ring%param%particle = antiparticle(ring%param%particle) call track_many (ring, orbtest_, 0, 0, -1) ring%param%particle = antiparticle(ring%param%particle) b = orbtest_(im2)%vec(1)/0.001 d = orbtest_(im2)%vec(2)/0.001 ! change sign of angle k44 = (prz_orbit(im2)%vec(1)/a - (-prz_orbit(im2)%vec(2)/c))/(b/a - d/c) k46 = prz_orbit(im2)%vec(1)/a - k44*b/a ring%ele( pretzel%ix_mag(1) )%value(hkick$) = k4 ring%ele( pretzel%ix_mag(2) )%value(hkick$) = k6 ring%ele( pretzel%ix_mag(3) )%value(hkick$) = k44 ring%ele( pretzel%ix_mag(4) )%value(hkick$) = k46 call lat_make_mat6(ring, -1, coord0) ! and get the orbit orbtest_(0)%vec(:) = 0. ! call closed_orbit_at_start (ring, orbtest_(0), 4, .true.) ! call track_all (ring, orbtest_) call closed_orbit_calc (ring, orbtest_, 4 ) pretzel%mag_hkick_theory(1) = k4 pretzel%mag_hkick_theory(2) = k6 pretzel%mag_hkick_theory(3) = k44 pretzel%mag_hkick_theory(4) = k46 ! print *,' West magpie kicks:' ! print *,' K4 = ',k4,' K6= ',k6,' K44 = ',k44,' k46= ',k46 ! Now do the East magpie ring%ele( pretzel%ix_mag(1) )%value(hkick$) = 0. ring%ele( pretzel%ix_mag(2) )%value(hkick$) = 0. ring%ele( pretzel%ix_mag(3) )%value(hkick$) = 0. ring%ele( pretzel%ix_mag(4) )%value(hkick$) = 0. ring%ele( pretzel%ix_mag(5) )%value(hkick$) = 0.001 call lat_make_mat6(ring, -1, coord0) orbtest_(0)%vec(:) = 0. call track_all (ring, orbtest_) a = orbtest_(im3)%vec(1)/0.001 c = orbtest_(im3)%vec(2)/0.001 ! print *,' orb = ',orbtest_(im3)%vec(1), orbtest_(im3)%vec(2) ring%ele( pretzel%ix_mag(5) )%value(hkick$) = 0.0 ring%ele( pretzel%ix_mag(6) )%value(hkick$) = 0.001 orbtest_(0)%vec(:) = 0. call lat_make_mat6(ring, -1, coord0) call track_all (ring, orbtest_) b = orbtest_(im3)%vec(1)/0.001 d = orbtest_(im3)%vec(2)/0.001 k44 = (prz_orbit(im3)%vec(1)/a - prz_orbit(im3)%vec(2)/c)/(b/a - d/c) ! print *,' k44 = ',k44 k46 = prz_orbit(im3)%vec(1)/a - k44*b/a ! print *,' k46 = ',k46 ring%ele( pretzel%ix_mag(5) )%value(hkick$) = k46 ring%ele( pretzel%ix_mag(6) )%value(hkick$) = k44 orbtest_(0)%vec(:) = 0. call lat_make_mat6(ring, -1, coord0) call track_all (ring, orbtest_) ! print *,' ***** At match point 1: X = ',orbtest_(im3)%vec(1), & ! ' Desired = ',prz_orbit(im3)%vec(1) ! print *,' Xp = ',orbtest_(im3)%vec(2), & ! ' Desired = ',prz_orbit(im3)%vec(2) ! print *,' ***** At match point 1: X = ',orbtest_(im4)%vec(1), & ! ' Desired = ',prz_orbit(im4)%vec(1) ! print *,' Xp = ',orbtest_(im4)%vec(2), & ! ' Desired = ',prz_orbit(im4)%vec(2) ring%ele( pretzel%ix_mag(5) )%value(hkick$) = 0.0 ring%ele( pretzel%ix_mag(6) )%value(hkick$) = 0.0 ring%ele( pretzel%ix_mag(8) )%value(hkick$) = 0.001 orbtest_(0)%vec(:) = 0. call lat_make_mat6(ring, -1, coord0) ring%param%particle = antiparticle(ring%param%particle) call track_many (ring, orbtest_, 0, 0, -1) ring%param%particle = antiparticle(ring%param%particle) a = orbtest_(im4)%vec(1)/0.001 c = orbtest_(im4)%vec(2)/0.001 ring%ele( pretzel%ix_mag(8) )%value(hkick$) = 0.0 ring%ele( pretzel%ix_mag(7) )%value(hkick$) = 0.001 orbtest_(0)%vec(:) = 0. call lat_make_mat6(ring, -1, coord0) ring%param%particle = antiparticle(ring%param%particle) call track_many (ring, orbtest_, 0, 0, -1) ring%param%particle = antiparticle(ring%param%particle) b = orbtest_(im4)%vec(1)/0.001 d = orbtest_(im4)%vec(2)/0.001 ! change sign of angle k6 = (prz_orbit(im4)%vec(1)/a - (-prz_orbit(im4)%vec(2)/c))/(b/a - d/c) k4 = prz_orbit(im4)%vec(1)/a - k6*b/a ring%ele( pretzel%ix_mag(7) )%value(hkick$) = k6 ring%ele( pretzel%ix_mag(8) )%value(hkick$) = k4 orbtest_(0)%vec(:) = 0. call lat_make_mat6(ring, -1, coord0) ring%param%particle = antiparticle(ring%param%particle) call track_many (ring, orbtest_, 0, 0, -1) ring%param%particle = antiparticle(ring%param%particle) ! print *,' ***** At match point 2: X = ',orbtest_(im4)%vec(1), & ! ' Desired = ',prz_orbit(im4)%vec(1) ! print *,' Xp = ',orbtest_(im4)%vec(2), & ! ' Desired = ',prz_orbit(im4)%vec(2) ! print *,' ***** At match point 1: X = ',orbtest_(im3)%vec(1), & ! ' Desired = ',prz_orbit(im3)%vec(1) ! print *,' Xp = ',orbtest_(im3)%vec(2), & ! ' Desired = ',prz_orbit(im3)%vec(2) ring%ele( pretzel%ix_mag(5) )%value(hkick$) = k46 ring%ele( pretzel%ix_mag(6) )%value(hkick$) = k44 ring%ele( pretzel%ix_mag(7) )%value(hkick$) = k6 ring%ele( pretzel%ix_mag(8) )%value(hkick$) = k4 call lat_make_mat6(ring, -1, coord0) ! and get the orbit orbtest_(0)%vec(:) = 0. ! call closed_orbit_at_start (ring, orbtest_(0), 4,.true.) ! call track_all (ring, orbtest_) call closed_orbit_calc (ring, orbtest_, 4 ) pretzel%mag_hkick_theory(5) = k46 pretzel%mag_hkick_theory(6) = k44 pretzel%mag_hkick_theory(7) = k6 pretzel%mag_hkick_theory(8) = k4 ! print *,' East magpie kicks:' ! print *,' K4 = ',k4,' K6= ',k6,' K55 = ',k44,' k53= ',k46 ! turn off magpies do i=1,8 ring%ele( pretzel%ix_mag(i) )%value(hkick$) = 0.0 enddo call lat_make_mat6(ring, -1, coord0) ! ! print *,' Exiting Rtn INIT_MAG_PRETZEL' end subroutine init_mag_pretzel !####################################################### ! choose phase space coordinates for an injected particle. ! if not random, then inject at nsig in x, y, dE ! if random, then select a particle at random from the phase space dist ! out to abs(nsig) ! subroutine inject_particle( ring, random, style, nsig, weight, inj, synch_beam, orbit) implicit none type (lat_struct) ring type (inj_struct) inj type (synch_beam_struct) synch_beam type (coord_struct) orbit(0:ring%n_ele_max) type (normal_modes_struct) mode integer i, style real nsig, weight, drele real vrannorm(6), voffset(6), sigma(3) logical random if( random ) then ! select randomly from phase space distribution sigma(1) = synch_beam%a%sigma ! H beamsize sigma(2) = synch_beam%b%sigma ! V beamsize sigma(3) = synch_beam%length ! bunch length call random_phase_space( style, vrannorm, sigma, nsig, weight ) ! Convert smearing vector from normalized coordinates to real 6-d coordinates ! X Error voffset(1) = vrannorm(1) ! X-prime Error voffset(2) = (vrannorm(2) - synch_beam%a%alpha*vrannorm(1)) / synch_beam%a%beta ! Y Error voffset(3) = vrannorm(3) ! Y-prime Error voffset(4) = (vrannorm(4) - synch_beam%b%alpha*vrannorm(3)) / synch_beam%b%beta ! Phase Error (z position) voffset(5) = vrannorm(5) ! Energy error voffset(6) = synch_beam%dele/synch_beam%length * vrannorm(6) ! Now offset phase space centroid according to energy error drele = synch_beam%centroid(6) + voffset(6) ! Energy offset voffset(1) = voffset(1) + drele*synch_beam%a%eta voffset(2) = voffset(2) + drele*synch_beam%a%etap voffset(3) = voffset(3) + drele*synch_beam%b%eta voffset(4) = voffset(4) + drele*synch_beam%b%etap else ! not random ! Inject particle offset from centroid by NSIG sigma ! Added dispersion contribution on 17 Feb 2015 - jac ! Only negative energy offsets will inject drele = synch_beam%centroid(6) - nsig*synch_beam%dele voffset(1) = - nsig*synch_beam%a%sigma voffset(2) = - synch_beam%a%alpha*voffset(1) / synch_beam%a%beta voffset(3) = nsig*synch_beam%b%sigma voffset(4) = - synch_beam%b%alpha*voffset(3) / synch_beam%b%beta voffset(5) = 0. voffset(6) = - nsig*synch_beam%dele ! Dispersion contribution voffset(1) = voffset(1) + drele*synch_beam%a%eta voffset(2) = voffset(2) + drele*synch_beam%a%etap voffset(3) = voffset(3) + drele*synch_beam%b%eta voffset(4) = voffset(4) + drele*synch_beam%b%etap endif do i = 1,6 orbit(inj%ix_injpt)%vec(i) = synch_beam%centroid(i) + voffset(i) enddo end subroutine inject_particle !####################################################### subroutine random_phase_space( style, vec, sigma, nsig, weight ) ! ! style = 0 gives uniform sampling of the 2-d phase space within a square of +-nsig ! style = 1 gives uniform sampling of the 2-d phase space within a radius of nsig ! style = 2 gives uniform sampling in r of the 2-d phase space ! style = 3 gives unweighted real 2-d phase space. Limit to +-nsig implicit none external hexist logical hexist real vec(6), nsig, weight, sigma(3), r, u, up, y, y1, y2 real phi, f, rnsig integer seed, i, style logical init_needed / .true. / save seed, init_needed if( init_needed ) then init_needed = .false. ! print *,' initializing random number generator' ! New random number generator y = myran( .true. ) ! true to init endif weight = 1.0 do i=1,3 if( style == 0 ) then! randomly select a point within a square of +-nsig y1 = myran( .false. ) u = -nsig + 2.*y1*nsig y2 = myran( .false. ) up = -nsig + 2.*y2*nsig r = sqrt(u**2 + up**2) vec(2*i-1) = u*sigma(i) vec(2*i) = up*sigma(i) weight = weight * 0.5*nsig**2 * exp( -r**2/2. ) elseif( style == 1 ) then! randomly select a point within a radius r = nsig r = 10.*nsig do while (r > nsig) y1 = myran( .false. ) u = -nsig + 2.*y1*nsig y2 = myran( .false. ) up = -nsig + 2.*y2*nsig r = sqrt(u**2 + up**2) enddo vec(2*i-1) = u*sigma(i) vec(2*i) = up*sigma(i) weight = weight * 0.5*nsig**2 * exp( -r**2/2. ) else if( style == 2 ) then ! sample uniformly in r y1 = myran( .false. ) r = y1*nsig y2 = myran( .false. ) phi = 6.28319 * y2 u = r * cos( phi ) up = r * sin( phi ) vec(2*i-1) = u*sigma(i) vec(2*i) = up*sigma(i) weight = weight * r * nsig * exp( -r**2/2. ) else if( style == 3 ) then ! pick randomly from phase-space (no weighting) ! ! Remove obsolete (pre-2009) penny-throwing algorithm ! 27 Feb 2015 jac ! y2 = 1.1 ! f = 1. ! do while( y2 > f ) ! acceptance-rejection on a gaussian ! y1 = myran( .false. ) ! r = y1*nsig ! f = exp( -r**2/2. ) ! y2 = myran( .false. ) ! enddo ! y2 = myran( .false. ) ! phi = 6.28319 * y2 ! u = r * cos( phi ) ! up = r * sin( phi ) ! ! Throw between +/- 6 sigma ! 30 nov 2009 jac Modify to limit range to +-nsig ! if (.not.hexist(3333)) call hbfun2(3333,' ',1000,-6.,6.,1000.,-6.,6.,g2d) ! Default to nsig=6 for backward compatibility if(nsig.le.0)nsig=6 rnsig=real(nsig) if (.not.hexist(3333))then print *,' Rtn RANDOM_PHASE_SPACE: Throw gauss-weighted over +-',nsig,' sigma' call hbfun2(3333,' ',1000,-rnsig,rnsig,1000.,-rnsig,rnsig,g2d) endif call hrndm2(3333,u,up) vec(2*i-1) = u*sigma(i) vec(2*i) = up*sigma(i) weight = 1.0 endif enddo end subroutine random_phase_space !####################################################### real function g2d(x,y) real x,y g2d=exp(-(x**2+y**2)/2.)/sqrt(twopi) end function !####################################################### subroutine pulsed_bump( ring, n, inj ) implicit none type (lat_struct) ring type (inj_struct) inj type (synch_beam_struct) synch_beam type (coord_struct), allocatable, save :: coord0(:) type (coord_struct), allocatable :: diff_co(:) character*40 fname integer n, i real arg ! allocate storage call reallocate_coord (coord0, ring%n_ele_max) call reallocate_coord (diff_co, ring%n_ele_max) do i= 1,3 arg = inj%omega*( real(n) *ring%param%total_length/c_light + inj%t0(i) ) if(arg>0.0 .and. arg< pi) then ring%ele(inj%ix_bumper(i))%value(hkick$) = inj%pb_coef(i)*sin(arg) else ring%ele(inj%ix_bumper(i))%value(hkick$) = 0.0 endif ! print *,' turn ',n,' bump ',i,' arg= ',arg,' kick = ', & ! ring%ele(inj%ix_bumper(i))%value(hkick$) ! print *,' t0 = ',inj%t0(i),' omega= ',inj%omega ! call lat_make_mat6(ring, inj%ix_bumper(i), coord0) call lat_make_mat6(ring, inj%ix_bumper(i)) enddo ! Now do the pinger if( n == inj%pingturn ) then ring%ele(inj%ix_pinger)%value(hkick$) = inj%pingkick else ring%ele(inj%ix_pinger)%value(hkick$) = 0.0 endif !call lat_make_mat6(ring, inj%ix_pinger, coord0) call lat_make_mat6(ring, inj%ix_pinger) ! Added 7/19/04 by DTK ! call tracking_update(ring, coord0, 4, diff_co) end subroutine pulsed_bump !####################################################### !NEW VERSION July 24,2001 subroutine pulsed_element_setup( ring, orbit, inj, synch_beam, & pb26e, pb28e, pb36e, dt, ntping, kping, pingcu) implicit none type (lat_struct) ring type (coord_struct) orbit(0:ring%n_ele_max) type (coord_struct), allocatable, save :: orbtest(:) type (coord_struct), allocatable, save :: coord0(:) type (coord_struct), allocatable, save :: diff_co(:) type (normal_modes_struct) mode type (twiss_struct) q34e_x, q34e_y type (twiss_struct) q32e_x, q32e_y type (twiss_struct) q07w_x, q07w_y type (inj_struct) inj type (synch_beam_struct) synch_beam type (ele_struct) ele type (pretzel_struct) pretzel integer i, j, ix1 integer bumploc(8), pingloc(2) integer nbump, nping integer q07wloc, q32eloc, q34eloc ! Pulsed bump setting in computer units real pb26e, pb28e, pb36e real r26,r28,r36,pblimit real bumpcu2rad real kping, dt integer ntping real pingcu real sigh, epsx, dele, sig_z real ph26, bx26, xk26, s26 real ph28, bx28, xk28, s28 real ph36, bx36, xk36, s36 real ph34, bx34, amp34, amp34_2 real ph32, bx32, amp32, amp32_2 real prtzl_pb_amp real max_amp, orbit_xamp real xinj, xpinj, pbamp, siginj, epsinj real l1, l2, xbamp ! allocate storage call reallocate_coord (orbtest, ring%n_ele_max) call reallocate_coord (coord0, ring%n_ele_max) call reallocate_coord (diff_co, ring%n_ele_max) ! ! Zero coord0 do i=0,ring%n_ele_max coord0(i)%vec(:) = 0. enddo ! Hardwired parameters inj%septthick = 0.003 inj%injerr = 0.002 inj%duration = 10.8e-6 ! 10.8 microseconds inj%omega = twopi/(2*inj%duration) ! Use 5 sigma clearance 2.6.03 jac inj%nsigcesr = 5.0 ! inj%nsigcesr = 7.5 inj%nsigsync = 2.6 ! Reduce clearances as per MGB PB lifetime calculations ! 1 Jan 2015 jac inj%nsigcesr = 4.5 inj%nsigsync = 2.0 ! 7 Jan 2015 jac inj%nsigcesr = 4.0 ! ! 4/9/03 Calibration for bumper is 3.593 mrad per 1000 cu at 1.885 GeV ! (3.593 mrad*GeV/1000 cu) ! 10/7/2003 Beam energy changed from gev to ev ! 10 sep 09 ! bumpcu2rad=3.593e+3*1.885/ring%E_TOT ! bumpcu2rad=3.593e+6*1.885/ring%ele(0)%value(E_TOT$) ! Easier to see energy correction. Result is rad/1000cu bumpcu2rad=3.593e-3*1.885e+9/ring%ele(0)%value(E_TOT$) ! Correct to 1 cu 1 Nov 2011 jac bumpcu2rad=bumpcu2rad/1000. !=========================================== ! Set up pinger. If a pinger setting in computer units ! is specified (pingcu), use it. Otherwise use the ! kick given in radians (kping). ! ! 4/4/03 pingcu2rad=(0.768/10**6)*1.885/ring%param%energy ! 10/7/2003 Beam energy changed from gev to ev ! 7/9/04 The calibration factor is about 0.8 mrad per 1000 cu at 1.89 GeV ! The software limit in the control system at present is 550 cu ! 10 sep 09 ! pingcu2rad=0.768e+3*1.885/ring%E_TOT ! print *,' PULSED_ELEMENT_SETUP: beam energy=',ring%E_TOT inj%pingcu2rad=0.768e+6*1.885/ring%ele(0)%value(E_TOT$) print *,' PULSED_ELEMENT_SETUP: beam energy=',ring%ele(0)%value(E_TOT$) if ( pingcu == 0 ) then inj%pingkick = kping/1000. else inj%pingkick = inj%pingcu2rad * pingcu endif print *,' PULSED_ELEMENT_SETUP: Input parameter KPING, ping kick(rad)=', & kping, inj%pingkick print *,' PULSED_ELEMENT_SETUP: pingcu2rad, ping kick(cu)=', & inj%pingcu2rad, inj%pingkick/inj%pingcu2rad inj%pingturn = ntping inj%delta_t = dt !! call emitt_calc( ring, all$, mode ) !! Replace emmit_calc with radiation_integrals. 20 june 2003 jac call radiation_integrals( ring, orbit, mode ) epsx = mode%a%emittance dele = mode%sige_e sig_z = mode%sig_z ! print *,' H emittance = ', epsx ! print *,' dE = ', dele,' bunch length = ',sig_z ! get locations of Q34E, bumpers and pinger nbump = 0 nping = 0 inj%ix_injpt = 0 inj%ix_injsex = 0 do i=1,ring%n_ele_max if ( ring%ele(i)%key == quadrupole$ .and. & ring%ele(i)%name == 'Q34E' ) then inj%ix_injpt = i q34eloc = i q34e_x%beta = ring%ele(i)%a%beta q34e_x%alpha = ring%ele(i)%a%alpha q34e_x%gamma = ring%ele(i)%a%gamma q34e_x%phi = ring%ele(i)%a%phi q34e_x%eta = ring%ele(i)%a%eta q34e_x%etap = ring%ele(i)%a%etap q34e_x%sigma = & sqrt( q34e_x%beta*epsx + dele**2 * q34e_x%eta**2) q34e_y%beta = ring%ele(i)%b%beta q34e_y%alpha = ring%ele(i)%b%alpha q34e_y%gamma = ring%ele(i)%b%gamma q34e_y%phi = ring%ele(i)%b%phi q34e_y%eta = ring%ele(i)%b%eta q34e_y%etap = ring%ele(i)%b%etap q34e_y%sigma = ring%ele(i)%b%sigma ! print *,' PULSED_ELEMENT_SETUP found quadrupole ',ring%ele(i)%name else if ( ring%ele(i)%key == quadrupole$ .and. & ring%ele(i)%name == 'Q32E' ) then q32eloc = i q32e_x%beta = ring%ele(i)%a%beta q32e_x%alpha = ring%ele(i)%a%alpha q32e_x%gamma = ring%ele(i)%a%gamma q32e_x%phi = ring%ele(i)%a%phi q32e_x%eta = ring%ele(i)%a%eta q32e_x%etap = ring%ele(i)%a%etap q32e_x%sigma = & sqrt( q32e_x%beta*epsx + dele**2 * q32e_x%eta**2) q32e_y%beta = ring%ele(i)%b%beta q32e_y%alpha = ring%ele(i)%b%alpha q32e_y%gamma = ring%ele(i)%b%gamma q32e_y%phi = ring%ele(i)%b%phi q32e_y%eta = ring%ele(i)%b%eta q32e_y%etap = ring%ele(i)%b%etap q32e_y%sigma = ring%ele(i)%b%sigma ! print *,' PULSED_ELEMENT_SETUP found quadrupole ',ring%ele(i)%name else if ( ring%ele(i)%key == quadrupole$ .and. & ring%ele(i)%name == 'Q07W' ) then q07wloc = i q07w_x%beta = ring%ele(i)%a%beta q07w_x%alpha = ring%ele(i)%a%alpha q07w_x%gamma = ring%ele(i)%a%gamma q07w_x%phi = ring%ele(i)%a%phi q07w_x%eta = ring%ele(i)%a%eta q07w_x%etap = ring%ele(i)%a%etap q07w_x%sigma = & sqrt( q07w_x%beta*epsx + dele**2 * q07w_x%eta**2) q07w_y%beta = ring%ele(i)%b%beta q07w_y%alpha = ring%ele(i)%b%alpha q07w_y%gamma = ring%ele(i)%b%gamma q07w_y%phi = ring%ele(i)%b%phi q07w_y%eta = ring%ele(i)%b%eta q07w_y%etap = ring%ele(i)%b%etap q07w_y%sigma = ring%ele(i)%b%sigma ! print *,' PULSED_ELEMENT_SETUP found quadrupole ',ring%ele(i)%name else if( ring%ele(i)%key == sextupole$ .and. & ring%ele(i)%name == 'SEX_34E' ) then print *,' Found sextupole at 34E i=',i inj%ix_injsex = i else if( ring%ele(i)%key == kicker$ .and. & ring%ele(i)%name(1:4) == 'BUMP' ) then if( ring%ele(i)%name == 'BUMP_26E' ) then inj%ix_bumper(1) = i nbump = nbump + 1 else if( ring%ele(i)%name == 'BUMP_28E' ) then inj%ix_bumper(2) = i nbump = nbump + 1 else if( ring%ele(i)%name == 'BUMP_36E' ) then inj%ix_bumper(3) = i nbump = nbump + 1 endif ! print *,' Found bumper ',ring%ele(i)%name else if( ring%ele(i)%key == kicker$ .and. & ! ring%ele(i)%name(1:10) == 'PINGER_36E' ) then ! Changed so that we are sure to get the superlord ring%ele(i)%name == 'PINGER_36E' ) then nping = nping + 1 inj%ix_pinger = i ! print *,' Found pinger ',ring%ele(i)%name endif enddo ! Did we find everything? if( inj%ix_injpt == 0 ) then print *,' Didnt find the injection point...exiting' stop else if( inj%ix_injsex == 0 ) then print *,' Didnt find sextupole s34E...exiting' stop else if( nbump /= 3 ) then print *,' Found ',nbump,' bumpers...exiting' stop else if( nping /= 1 ) then print*,' Found ',nping,' pingers...exiting' stop endif ! Injection point parameters write(6, 10) q34e_x%beta, q34e_y%beta, q34e_x%eta 10 format(2x,'Routine PULSED_ELEMENT_SETUP: betax = ',f6.3, & ' betay= ',f6.3,' etax= ',f6.3) write(6,20) 1000.*q34e_x%sigma 20 format(2x,'Routine PULSED_ELEMENT_SETUP: Horizontal beamsize = ', & f7.4,' mm ') ! Calculate pulsed bump coefficients in pretzel-off optics !-------------------------------------------------------------- ! 30 may 2007 - jac - This obsolete attempt does not allow for ! a wiggler orbit, nor does it turn off BBI ! call lat_make_mat6( ring, -1, coord0) ! call twiss_at_start( ring ) ! call twiss_propagate_all( ring ) ! !-------------------------------------------------------------- ring%param%particle = electron$ ! Turn BBI elements off call lrbbi_ele_setup ( ring, .false., .false. ) ! Turn pretzel off ! SET_PRETZEL recalculates the transfer matrices and twiss functions call set_pretzel(0, ring, pretzel) ! Check if these are super-lords (means a pc landed in this element) do i=1,3 ele = ring%ele(inj%ix_bumper(i)) ! 10 sep 09 ! if( ele%control_type == super_lord$ ) then ! just use first slave if( ele%lord_status == super_lord$ ) then ! just use first slave ix1 = ring%control(ele%ix1_slave)%slave%ix_ele print *, 'Found a super lord for bumper, using ', & ring%ele(ix1)%name else ix1 = inj%ix_bumper(i) endif if( i == 1 ) then !26 bx26 = ring%ele(ix1)%a%beta ph26 = ring%ele(ix1)%a%phi s26 = ring%ele(ix1)%s else if( i == 2 ) then !28 bx28 = ring%ele(ix1)%a%beta ph28 = ring%ele(ix1)%a%phi s28 = ring%ele(ix1)%s else if( i == 3 ) then !36 bx36 = ring%ele(ix1)%a%beta ph36 = ring%ele(ix1)%a%phi s36 = ring%ele(ix1)%s endif enddo ! end loop over bumpers to get s, betas and phases print *,' s36= ',s36,' s28= ',s28,' s26= ',s26 ! Put optics back to pretzel-on !------------------------------------------------------------------------- ! 30 may 2007 - jac. Modified to account for wiggler orbit and BBI ! call lat_make_mat6( ring, -1, orbit) ! call twiss_at_start( ring ) ! call twiss_propagate_all( ring ) ! ! Turn pretzel back on ! SET_PRETZEL recalculates the transfer matrices and twiss functions call set_pretzel(1, ring, pretzel) ! Turn BBI elements back on call lrbbi_ele_setup ( ring, .true., .true. ) ! Put orbit in coord0 call tracking_update (ring, coord0, 4, diff_co) !------------------------------------------------------------------------- ! Preliminary calculation of relative bumper kicks xk26 = 1.0 xk28 = -xk26*sqrt(bx26)*sin(ph36-ph26)/(sqrt(bx28)*sin(ph36-ph28)) xk36 = -(xk26*sqrt(bx26)*cos(ph36-ph26) & +xk28*sqrt(bx28)*cos(ph36-ph28))/sqrt(bx36) print *, 'Relative Lattice Bump coefficients from Twiss functions: 26, 28, 36:', xk26,xk28,xk36 ! write(13,30) 'Bump coefficients: 26, 28, 36:', xk26,xk28,xk36 30 format(a,3f9.4) !================================================================ ! If pulsed bump settings are available in the initialization ! data file (pb26e <> 0), use them. Otherwise calculate and ! use the bump magnitude which yields the desired clearance. !================================================================ if (pb26e == 0 ) then bx32 = q32e_x%beta ph32 = q32e_x%phi bx34 = q34e_x%beta ph34 = q34e_x%phi ! Assume arbitrary pb kick values are in mrad (xk26 = 1 mrad) amp32 = xk26*sqrt(bx26*bx32)*sin(-(ph32-ph26))+ & xk28*sqrt(bx28*bx32)*sin(-(ph32-ph28)) amp32_2 = xk36*sqrt(bx36*bx32)*sin(ph32-ph36) amp32=amp32/1000. amp32_2=amp32_2/1000. amp34 = xk26*sqrt(bx26*bx34)*sin(-(ph34-ph26))+ & xk28*sqrt(bx28*bx34)*sin(-(ph34-ph28)) amp34_2 = xk36*sqrt(bx36*bx34)*sin(ph34-ph36) amp34=amp34/1000. amp34_2=amp34_2/1000. ! Now figure out the bump timings ! Assume bump timings are set to be simultaneous for positrons l1 = abs( s36 - s28 ) l2 = abs( s28 - s26 ) inj%t0(1) = -2.*(l1+l2)/c_light + inj%delta_t + pi/(2.*inj%omega) inj%t0(2) = -2.*l1/c_light + inj%delta_t + pi/(2.*inj%omega) inj%t0(3) = inj%delta_t + pi/(2.*inj%omega) ! Now get pretzel amplitude at Q34E print *,' Pretzel amplitude (m) and angle (radians) at Q34E = ', & orbit(inj%ix_injpt)%vec(1), orbit(inj%ix_injpt)%vec(2) ! Now add the pulsed bump to it ! Check amplitude at injection with pretzel and pulsed_bump using ! these kicks calculated with XK26=1.0. ! If these still-unscaled kick values result in a pretzel+bump amplitude ! at 34E which is greater than the value of MAX_AMP, then recalcalculate ! them to limit the orbit to MAX_AMP. This will allow to get an orbit ! with pretzel and bump on. That orbit will be used to normalize the ! kicks to satisfy the NSIGCESR criterion for the distance to the ! inner vacuum chamber wall. orbit_xamp = orbit(inj%ix_injpt)%vec(1) ! Maximum amplitude allowed for unscaled kick calculation. Otherwise ! kicks are rescaled to reduce the orbit at 34E to allow orbit calculation ! for final scaling of pulsed bump kicks. max_amp = 0.040 prtzl_pb_amp = amp34_2 + orbit_xamp print *,'Before bump optimization: amp34 ',amp34,' amp34_2 = ',amp34_2, '(m)' print *,'orbit_xamp without bump',orbit_xamp print *,'prztl_pb_amp (before bump optimization)',prtzl_pb_amp orbit_xamp=abs(orbit_xamp) if( abs(prtzl_pb_amp) > abs(max_amp) ) then ! Algebraic mistake in calculation of xk26 corrected ! 13 October 2003 JAC ! Derivation: Calculate linear relationship between xk36 and xk26 ! using bump equalities (see above for relative calculation) ! Substitute xk36 value which gives max_amp, as shown here below, ! to get xk26. ! xk26 = -(max_amp - orbit_xamp)/(sqrt(bx26*bx34) * & ! sin(ph36-ph34)*(cos(ph36-ph26) + (sign mistake here) sin(ph36-ph26) * & ! cos(ph36-ph28) / sin(ph36-ph28))) xk26 = -(max_amp - orbit_xamp)/(sqrt(bx26*bx34) * & sin(ph36-ph34)*(cos(ph36-ph26) - sin(ph36-ph26) * & cos(ph36-ph28) / sin(ph36-ph28))) xk28 = -xk26*sqrt(bx26)*sin(ph36-ph26) / & (sqrt(bx28)*sin(ph36-ph28)) xk36 = (max_amp - orbit_xamp)/(sqrt(bx34*bx36)*sin(ph36-ph34)) print *,' Renormalized kick values to limit pb contribution to orbit at 34W to ',max_amp, 'm' print *,' New XK26 = ',xk26, 'rad' print *,' New XK28 = ',xk28, 'rad' print *,' New XK36 = ',xk36, 'rad' else ! convert to radians xk26 = xk26/1000. xk28 = xk28/1000. xk36 = xk36/1000. endif amp34 = xk26*sqrt(bx26*bx34)*sin(-(ph34-ph26))+ & xk28*sqrt(bx28*bx34)*sin(-(ph34-ph28)) amp34_2 = xk36*sqrt(bx36*bx34)*sin(ph34-ph36) amp32 = xk26*sqrt(bx26*bx32)*sin(-(ph32-ph26))+ & xk28*sqrt(bx28*bx32)*sin(-(ph32-ph28)) amp32_2 = xk36*sqrt(bx36*bx32)*sin(ph32-ph36) print *,'Prior to final bumper kick normalization::' print *,' Linear approx PB Amplitude at Q34E (amp_34, from 26 & 28): ',amp34, 'm' print *,' Linear approx PB Amplitude at Q34E (amp34_2, from 36): ',amp34_2, 'm' prtzl_pb_amp = amp34_2 + orbit(inj%ix_injpt)%vec(1) print *,' Pretzel + Pulsed Bump Amplitude(amp34_2+orbit): ',prtzl_pb_amp, 'm' print *,' Linear approx Pulsed Bump Amplitude at Q32E: ',amp32_2, 'm' prtzl_pb_amp = amp32_2 + orbit(inj%ix_injpt)%vec(1) print *,' Pretzel + Pulsed Bump Amplitude(amp32_2+orbit): ',prtzl_pb_amp, 'm' !================================================================================ ! Final bumper kick normalization. ring%ele(inj%ix_bumper(1))%value(hkick$) = xk26 ring%ele(inj%ix_bumper(2))%value(hkick$) = xk28 ring%ele(inj%ix_bumper(3))%value(hkick$) = xk36 ! 30 may 2007 - jac - Redo the pretzel- on optics with the new bumper kicks. ! Include BBI effects. do i=1,3 call lat_make_mat6(ring, inj%ix_bumper(i), coord0) enddo ! and get the orbit. Don't enforce aperture limits for this test orbit. bmad_com%aperture_limit_on=.false. orbtest(0)%vec(:) = 0. ! call closed_orbit_at_start (ring, orbtest(0), 4,.true.) ! call track_all (ring, orbtest) call closed_orbit_calc (ring, orbtest, 4 ) bmad_com%aperture_limit_on=.true. ! and measure the pulsed bump amplitude xbamp = orbtest(inj%ix_injpt)%vec(1) - orbit(inj%ix_injpt)%vec(1) print *,'Tracked Pulsed Bump Amplitude before scaling (xbamp=orbtest-orbit): ',xbamp, 'm' ! Now calculate the maximum pulsed bump amplitude which satsisfies the clearance criterion pbamp = -0.045 - (orbit(inj%ix_injpt)%vec(1) - inj%nsigcesr*q34e_x%sigma) print *,'nsigcesr = ',inj%nsigcesr,' q34_x%sigma = ',q34e_x%sigma print *,'Pulsed bump amplitude satisfying clearance criterion (pbamp):',pbamp, 'm' ! and scale appropriately xk26 = pbamp/xbamp*xk26 xk28 = pbamp/xbamp*xk28 xk36 = pbamp/xbamp*xk36 !----------------------------------------------------------------------------- ! Impose control system limits on bumper kick amplitudes. ! 1 November 2011 jac pb26e = xk26/bumpcu2rad pb28e = xk28/bumpcu2rad pb36e = xk36/bumpcu2rad ! Hardware limit from MGB 1 Nov 2011 pblimit=600. print *,'Pulsed bump kicks (rad) and cu settings before hardware limit: ' write(6,1100)xk26,pb26e, & xk28,pb28e, & xk36,pb36e 1100 format(' 26E:',e10.3,2x,e10.3/, & ' 28E:',e10.3,2x,e10.3/, & ' 36E:',e10.3,2x,e10.3 & ) r26 = abs(pb26e/pblimit) if (r26.gt.1.)then print *,' PB26E exceeds hardware limit. Scaling all settings by factor ',1/r26 pb26e = pb26e / r26 pb28e = pb28e / r26 pb36e = pb36e / r26 endif r28 = abs(pb28e/pblimit) if (r28.gt.1.)then print *,' PB28E exceeds hardware limit. Scaling all settings by factor ',1/r28 pb26e = pb26e / r28 pb28e = pb28e / r28 pb36e = pb36e / r28 endif r36 = abs(pb36e/pblimit) if (r36.gt.1.)then print *,' PB36E exceeds hardware limit. Scaling all settings by factor ',1/r36 pb26e = pb26e / r36 pb28e = pb28e / r36 pb36e = pb36e / r36 endif xk26 = bumpcu2rad * pb26e xk28 = bumpcu2rad * pb28e xk36 = bumpcu2rad * pb36e print *,'After pulsed bump optimization and hardware limit:' amp34_2 = xk36*sqrt(bx36*bx34)*sin(ph34-ph36) print *,' Linear approx Pulsed Bump Amplitude at Q34E: ',amp34_2, 'm' prtzl_pb_amp = amp34_2 + orbit(inj%ix_injpt)%vec(1) print *,' Pretzel + Pulsed Bump Amplitude AT Q34E: ',prtzl_pb_amp, 'm' amp32_2 = xk36*sqrt(bx36*bx32)*sin(ph32-ph36) print *,' Linear approx Pulsed Bump Amplitude at Q32E: ',amp32_2, 'm' prtzl_pb_amp = amp32_2 + orbit(inj%ix_injpt)%vec(1) print *,' Pretzel + Pulsed Bump Amplitude AT Q32E: ',prtzl_pb_amp, 'm' !----------------------------------------------------------------------------- ! Load final pulsed bump coefficients inj%pb_coef(1) = xk26 inj%pb_coef(2) = xk28 inj%pb_coef(3) = xk36 !================================================================================ else ! Use pulsed-bump settings (cu) from setup file print *,' Using pulsed-bump settings (cu) from setup file ' Print *,' Bump calibration factor is ',bumpcu2rad,' cu/rad' ! bumpcu2rad gives radians from cu xk26 = bumpcu2rad*pb26e xk28 = bumpcu2rad*pb28e xk36 = bumpcu2rad*pb36e ! Now figure out the bump timings l1 = abs( s36 - s28 ) l2 = abs( s28 - s26 ) inj%t0(1) = -2.*(l1+l2)/c_light + inj%delta_t + pi/(2.*inj%omega) inj%t0(2) = -2.*l1/c_light + inj%delta_t + pi/(2.*inj%omega) inj%t0(3) = inj%delta_t + pi/(2.*inj%omega) ! bx34 = q34e_x%beta ph34 = q34e_x%phi bx32 = q32e_x%beta ph32 = q32e_x%phi amp34_2 = xk36*sqrt(bx36*bx34)*sin(ph34-ph36) print *,'Linear approx Pulsed Bump Amplitude at Q34E: ',amp34_2, 'm' prtzl_pb_amp = amp34_2 + orbit(inj%ix_injpt)%vec(1) print *,'Pretzel + Pulsed Bump Amplitude AT Q34E: ',prtzl_pb_amp, 'm' amp32_2 = xk36*sqrt(bx36*bx32)*sin(ph32-ph36) print *,'Linear approx Pulsed Bump Amplitude at Q32E: ',amp32_2, 'm' prtzl_pb_amp = amp32_2 + orbit(inj%ix_injpt)%vec(1) print *,'Pretzel + Pulsed Bump Amplitude AT Q32E: ',prtzl_pb_amp, 'm' ! Enter the pulsed bump coefficients inj%pb_coef(1) = xk26 inj%pb_coef(2) = xk28 inj%pb_coef(3) = xk36 endif ! Endif statement for choosing type of pulsed bump calculation print *,'Pulsed bump Coefficients (kicks in radians) and control system setting (cu): ' ! do i=1,3 ! print *,'pb_coef number ',i,' = ',inj%pb_coef(i) write(6,1111)inj%pb_coef(1),pb26e, & inj%pb_coef(2),pb28e, & inj%pb_coef(3),pb36e 1111 format(' 26E:',e10.3,2x,e10.3/, & ' 28E:',e10.3,2x,e10.3/, & ' 36E:',e10.3,2x,e10.3 & ) ! enddo ! Check the pulsed bump ring%ele(inj%ix_bumper(1))%value(hkick$) = inj%pb_coef(1) ring%ele(inj%ix_bumper(2))%value(hkick$) = inj%pb_coef(2) ring%ele(inj%ix_bumper(3))%value(hkick$) = inj%pb_coef(3) do i=1,3 call lat_make_mat6(ring, inj%ix_bumper(i), coord0) enddo ! and get the orbit orbtest(0)%vec(:) = 0. ! call closed_orbit_at_start (ring, orbtest(0), 4,.true.) ! call track_all (ring, orbtest) call closed_orbit_calc (ring, orbtest, 4 ) print *, '*** Final tracked pulsed bump calculation:' print *,'Pulsed bump amplitude= ', & orbtest(inj%ix_injpt)%vec(1) - orbit(inj%ix_injpt)%vec(1) print *,'Stored Beam x, x-prime at pulsed bump maximum: ', & orbtest(inj%ix_injpt)%vec(1), orbtest(inj%ix_injpt)%vec(2) print *,'Stored beam y, yprime at pulsed-bump maximum:', & orbtest(inj%ix_injpt)%vec(3),orbtest(inj%ix_injpt)%vec(4) print *,' Pretzel + pulsed bump + ',inj%nsigcesr,' sigma= ', & orbtest(inj%ix_injpt)%vec(1) - inj%nsigcesr*q34e_x%sigma ! print *,'Tracked Stored-Beam x position at Q32E: ', & orbtest(q32eloc)%vec(1) ! Record the Stored beam x, x-prime at pulsed bump maximum inj%xpb_stored(1) = orbtest(inj%ix_injpt)%vec(1) inj%xpb_stored(2) = orbtest(inj%ix_injpt)%vec(2) ! Record the Stored beam y, y-prime at pulsed bump maximum inj%ypb_stored(1) = orbtest(inj%ix_injpt)%vec(3) inj%ypb_stored(2) = orbtest(inj%ix_injpt)%vec(4) ! Record the Stored beam x, x-prime at bumper 36E inj%x36_stored(1) = orbtest(inj%ix_bumper(3))%vec(1) inj%x36_stored(2) = orbtest(inj%ix_bumper(3))%vec(2) ! Print out the stored beam X position at Q07W print *,' PULSED_ELEMENT_SETUP: s, x at Q07W:', & ring%ele(q07wloc)%s, orbtest(q07wloc)%vec(1) ! Now turn off the pulsed bump ring%ele(inj%ix_bumper(1))%value(hkick$) = 0. ring%ele(inj%ix_bumper(2))%value(hkick$) = 0. ring%ele(inj%ix_bumper(3))%value(hkick$) = 0. do i=1,3 call lat_make_mat6(ring, inj%ix_bumper(i), coord0) enddo ! and the pinger ring%ele(inj%ix_pinger)%value(hkick$) = 0. call lat_make_mat6(ring, inj%ix_pinger, coord0) end subroutine pulsed_element_setup !####################################################### subroutine synch_beam_setup( ring, inj, mode, synch_beam, use_default_centroid, & synch_centroid, use_default_twiss, & synch_twiss, synch_epsx, synch_epsy, & ntestturns, angmin, angstep, nangsteps) implicit none type (lat_struct) ring type (inj_struct) inj type (normal_modes_struct) mode type (synch_beam_struct) synch_beam type (twiss_struct) tw real synch_epsx, synch_epsy, dele, sig_z, xinj, xpinj real synch_twiss(6), synch_centroid(6) integer i logical use_default_centroid, use_default_twiss ! integer ntestturns, nangsteps real angmin, angstep real inj_angle logical error ! Now define the incoming synchrotron beam synch_beam%epsx = synch_epsx synch_beam%epsy = synch_epsy ! Use CESR bunch length and energy spread ! 11 March 2015 jac dele = mode%sige_e sig_z = mode%sig_z synch_beam%dele = dele synch_beam%length = sig_z if( use_default_twiss ) then ! take injection point twiss parameters ! Assume Synch/XFR/CESR match is perfect for the moment tw = ring%ele(inj%ix_injpt)%a synch_beam%a%beta = tw%beta synch_beam%a%alpha = -tw%alpha ! flip sign since e- in other direc synch_beam%a%gamma = tw%gamma synch_beam%a%phi = tw%phi synch_beam%a%eta = tw%eta synch_beam%a%etap = -tw%etap ! flip sign synch_beam%a%sigma = & sqrt( synch_beam%a%beta*synch_beam%epsx + synch_beam%dele**2 * & synch_beam%a%eta**2) tw = ring%ele(inj%ix_injpt)%b synch_beam%b%beta = tw%beta synch_beam%b%alpha = -tw%alpha ! flip sign synch_beam%b%gamma = tw%gamma synch_beam%b%phi = tw%phi synch_beam%b%eta = 0.0 synch_beam%b%etap = 0.0 ! flip sign synch_beam%b%sigma = sqrt( synch_beam%b%beta*synch_beam%epsy ) else ! use provided twiss parameters synch_beam%a%beta = synch_twiss(1) synch_beam%a%alpha = synch_twiss(2) synch_beam%a%gamma = (1. + synch_beam%a%alpha**2)/synch_beam%a%beta synch_beam%a%phi = 0. synch_beam%a%eta = synch_twiss(5) synch_beam%a%etap = synch_twiss(6) synch_beam%a%sigma = & sqrt( synch_beam%a%beta*synch_beam%epsx + synch_beam%dele**2 * & synch_beam%a%eta**2) synch_beam%b%beta = synch_twiss(3) synch_beam%b%alpha = synch_twiss(4) synch_beam%b%gamma = (1. + synch_beam%b%alpha**2)/synch_beam%b%beta synch_beam%b%phi = 0.0 synch_beam%b%eta = 0.0 synch_beam%b%etap = 0.0 ! flip sign synch_beam%b%sigma = sqrt( synch_beam%b%beta*synch_beam%epsy ) endif print *,' synch beam X: emitt, beta, alpha, eta, etap', & synch_beam%epsx, synch_beam%a%beta, synch_beam%a%alpha, & synch_beam%a%eta, synch_beam%a%etap print *,' synch beam Y: emitt, beta, alpha, eta, etap', & synch_beam%epsy, synch_beam%b%beta, synch_beam%b%alpha, & synch_beam%b%eta, synch_beam%b%etap print *,' synch beam: sig_E/E = ',synch_beam%dele,' sig_z= ', & synch_beam%length print *,' Synch beam sigmas (x,y,z): ',synch_beam%a%sigma, & synch_beam%b%sigma, synch_beam%length print *,' Synch beam sigmas (xp,yp,zp): ', & synch_beam%a%sigma/synch_beam%a%beta, & synch_beam%b%sigma/synch_beam%b%beta, synch_beam%dele ! Now calculate the amplitudes if( use_default_centroid ) then do i=1,6 synch_beam%centroid(i) = 0.0 enddo ! ! Always use custom energy offset, even when default centroid chosen synch_beam%centroid(6) = synch_centroid(6) ! ! Set vertical injection and position equal to stored, bumped beam values ! Inj angle is negative. synch_beam%centroid(3) = inj%ypb_stored(1) synch_beam%centroid(4) = - inj%ypb_stored(2) ! Choose injection position to clear septum wall by NSIGSYNC sigma synch_beam%centroid(1) = -0.045 - inj%nsigsync*synch_beam%a%sigma - & inj%septthick - inj%injerr ! The difference from the stored pulsed-bump orbit is the amplitude ! of the injection oscillation xinj = synch_beam%centroid(1) - inj%xpb_stored(1) ! One can use the CESR twiss functions at the injection point to ! calculate the injection angle, even if matching is not assumed ! and custom synch beam twiss values are used. ! However, this approximation does not take into account an offset ! in the injected beam energy. Also, it does not account for the ! pretzel amplitude at bumper 36E. The optimal injection angle ! puts the injected orbit at the position of the stored pretzel&pb ! orbit at bumper 36E corrected for the energy offset. 11 Aug 2003 jac ! xpinj = -xinj*( -ring%ele(inj%ix_injpt)%a%alpha/ring%ele(inj%ix_injpt)%a%beta ) synch_beam%centroid(2) = -inj%xpb_stored(2) + xpinj print *,' Xprime calculated using CESR twiss = ',synch_beam%centroid(2) ! ! Use the injection angle determined via scanning survival efficiency, ! if available, even if it does not satisfy the turns criterion. ! Otherwise, use the calculation by INJECTION_ANGLE. ! If the latter is also not available, use the twiss function approximation call injangle_scan ( ring, inj, synch_beam, .true., & ntestturns, angmin, angstep, nangsteps, & inj_angle, error ) if ( inj_angle.ne.0. ) then synch_beam%centroid(2) = inj_angle xpinj = synch_beam%centroid(2) + inj%xpb_stored(2) else ! ! Rtn INJECTION_ANGLE calculates the injection angle which puts the injected ! orbit at the position of the stored pretzel&pb orbit at bumper 36E call injection_angle ( ring, inj, synch_beam, & inj_angle, error ) if ( .not. error ) then ! inj_angle is negative. xpb_stored(2) is positive. synch_beam%centroid(2) = inj_angle xpinj = synch_beam%centroid(2) + inj%xpb_stored(2) else print *, 'Error calculating injection angle. Use CESR twiss approximation.' endif endif ! print *,'Amplitudes of injection oscillation: Xinj = ',xinj,' Xpinj = ',xpinj ! ! synch_beam%centroid(2) = -3.6e-3 else ! Here if using custom sync beam centroid values from initialization file do i=1,6 synch_beam%centroid(i) = synch_centroid(i) enddo endif ! Final values used for injection angle and position print *, ' After SYNCH_BEAM_SETUP:' print *,' X, Xprime of injected bunch used to inject = ', & synch_beam%centroid(1), synch_beam%centroid(2) print *,' Y, Yprime of injected bunch used to inject = ', & synch_beam%centroid(3), synch_beam%centroid(4) end subroutine synch_beam_setup !####################################################### !+ ! Subroutine injection_angle ( ring, inj, synch_beam, inj_angle, error ) ! ! Calculate injection angle by scanning through angle ! and checking survival efficiency 15 October 2003 JAC ! ! Modules Needed: ! ! Input: ! ring -- lat_struct: Ring ! inj -- injection parameters ! synch_beam -- synch beam parameters ! ntesturns, angmin, angstep, nangsteps -- Angle scan parameters ! ! Output: ! inj_angle -- Real: injection angle for electrons (radians, negative) ! error -- Logical: true if calculation unsuccessful subroutine injangle_scan ( ring, inj, synch_beam, write_ntuple, & ntestturns, angmin, angstep, nangsteps, & inj_angle, error ) implicit none type (lat_struct) ring type (inj_struct) inj type (synch_beam_struct) synch_beam logical write_ntuple integer ntestturns, nangsteps real angmin, angstep real inj_angle logical error ! type (coord_struct), allocatable, save :: inj_orbit(:) ! real injang character datfname*40 logical lostit integer i, n, ix_start, lun5 integer indprint integer norbwrt1,norbwrt2 integer ixlost, t_state integer nstep, npts, itrackdir real weight ! real angmax real survangmin, survangmax ! integer ntbest, ixbest ! print *,'Entering Rtn INJANGLE_SCAN ...' ! Allocate storage call reallocate_coord (inj_orbit, ring%n_ele_max) ! allocate (inj_orbit(0:ring%n_ele_max)) bmad_com%aperture_limit_on = .true. call inject_particle( ring, .false., 0, 0., weight, inj, & synch_beam, inj_orbit ) angmax= angmin + nangsteps*angstep print *,' Rtn INJANGLE_SCAN called with min, max, step:',& angmin,angmax,angstep survangmin=0. survangmax=0. ! Keep track of furthest progress if ! number of turns criterion not satisfied ! in order to return best angle ixbest = ring%n_ele_track ntbest = -1 nstep = 0 if ( angstep .eq. 0 ) then angmax = angmin angstep = 1. endif ! Initialized optimized angle inj_angle = 0. injang = angmin - angstep do injang = injang + angstep if (injang > angmax) exit nstep = nstep + 1 ! lun5=lunget() ! write(datfname,4500) 4500 format('dat/inj_angscan_orbits.dat') ! open(lun5, file=datfname, status = 'unknown') ! Set central injected angle for zero energy offset ! Rtn INJECT_PARTICLE will return in INJ_ORBIT the value of the ! injection angle including the dispersive contribution ! The optimized value of the injection angle will not include ! the dispersive part. synch_beam%centroid(2) = injang call inject_particle( ring, .false., 0, 0., weight, inj, & synch_beam, inj_orbit ) ! Take the beam to the IP ix_start = inj%ix_injpt norbwrt2 = inj%ix_injpt n = 0 lostit = .false. ! Set apertures for Q34E and S34E larger on first turn for injection channel ! 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 ! print *,i,ring%ele( i )%name(1:7) 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 ! print *,' 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 TRACK_MANY, 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) do while( n < ntestturns .and. .not. lostit ) 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 = t_state) ring%param%particle = antiparticle(ring%param%particle) ! do i=0, norbwrt2 ! write(lun5,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 ! Fill track ntuple for first two turns ! Track electrons backwards itrackdir = -1 npts = mod ( ix_start + ring%n_ele_track -1 , ring%n_ele_track ) + 1 if (t_state /= moving_forward$) then ! Confusing nomenclature here. exit_end and entrance_end are defined ! to be for positrons We are tracking electrons here, so ! tracking backward. If a particle is lost at the "exit_end" (which ! is really the entrance end for electrons), then there is no ! orbit data for the following element and we don't include it ! in the ntuple. if (inj_orbit(t_state)%location .eq. exit_end$ ) then npts = npts - t_state + 1 else npts = npts - t_state + 2 endif endif if ( write_ntuple .and. n < 2 .and. npts > 0 ) call track_nt(0, ring, inj, nstep, n, inj_orbit, ix_start, npts, itrackdir) ! For turns after turn zero, track around entire ring ix_start = 0 norbwrt2 = ring%n_ele_track ! Increment turn number n = n+1 ! Set apertures for Q34E and S34E back to normal for turns after turn 0 ! 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 ! print *,' Resetting aperture limit to 4.5 cm for element ',i,ring%ele(i)%name endif enddo if (t_state /= moving_forward$) then lostit=.true. ! Jump out of turns loop if particle lost ixlost = t_state endif enddo ! End of loop over turns ! close(lun5) if( lostit ) then write(6,5100) n-1, 1000*injang, ixlost, & ring%ele(ixlost)%name, ring%ele(ixlost)%s, & 100*inj_orbit(ixlost)%vec(1), & 100*inj_orbit(ixlost)%vec(3), & 100*ring%ele(ixlost)%value(x1_limit$), & 100*ring%ele(ixlost)%value(y1_limit$) 5100 format(' Lost on turn ',i2,' for injection angle ',f6.2,& ' mrad at element ', i4, & a9,' s:',f6.1,' X:',f5.2,' Y:',f5.2, & ' Xapert:',f5.2,' Yapert:',f5.2) ! print *,' Lost on turn ',n-1,' for angle= ',injang ! print *,' elem ',ixlost,ring%ele( ixlost )%name ! print *,' x,y lim = ',ring%ele( ixlost )%value(x1_limit$), & ! ring%ele( ixlost )%value(y1_limit$) ! Return angle of furthest progress if turns criterion ends up not being satisfied if( n-1.gt.ntbest .or. (n-1.eq.ntbest .and. ixlost.lt.ixbest) ) then ! Set returned angle to angle of furthest progress ! Will be overridden if turns criterion met inj_angle = injang ntbest = n-1 ixbest = ixlost write(6,5150)1000*injang,ring%ele(ixlost)%name,ntbest 5150 format(' New angle of furthest progress:',f6.2, & ' mrad at element ',a9,' on turn ',i3) endif else print *,' Beam survived ',ntestturns,' turns for angle= ',injang if ( survangmin .eq. 0. ) survangmin = injang survangmax = injang endif enddo ! End of loop over angles ! ! Close track ntuple file if ( write_ntuple ) call track_nt(1, ring, inj, 0, 0, inj_orbit, 0, 0, 0) ! ! If survival angle found, return center of survival range as ! best injection angle if ( survangmin .eq. 0. ) then error = .true. print *,'Rtn INJANGLE_SCAN finds no survival injection angle for ', & ntestturns, ' turns' else error = .false. inj_angle = ( survangmax + survangmin ) / 2. print *,'Rtn INJANGLE_SCAN finds survival range (mrad):',survangmin,survangmax print *,'Rtn INJANGLE_SCAN returns optimal injection angle:',inj_angle endif ! Make sure pulsed bump and pinger are off ! If the orbit did not survive on a turn where the ! bump or pinger were on, they could still be left ! on by PULSED_BUMP at this point. ring%ele(inj%ix_bumper(1))%value(hkick$) = 0. ring%ele(inj%ix_bumper(2))%value(hkick$) = 0. ring%ele(inj%ix_bumper(3))%value(hkick$) = 0. do i=1,3 call lat_make_mat6(ring, inj%ix_bumper(i)) enddo ring%ele(inj%ix_pinger)%value(hkick$) = 0. call lat_make_mat6(ring, inj%ix_pinger) end subroutine injangle_scan !####################################################### !+ ! Subroutine injection_angle ( ring, inj, synch_beam, inj_angle, error ) ! ! Calculate injection angle 12 August 2003 JAC ! ! Modules Needed: ! ! Input: ! ring -- lat_struct: Ring ! inj -- injection parameters ! synch_beam -- synch beam parameters ! ! Output: ! inj_angle -- Real: injection angle for electrons (radians, negative) ! error -- Logical: true if calculation unsuccessful subroutine injection_angle ( ring, inj, synch_beam, inj_angle, error ) implicit none type (lat_struct) ring type (inj_struct) inj type (synch_beam_struct) synch_beam real inj_angle logical error real(rp) tmat(6,6) real xstored, xpstored, x36e integer i,nt ! need to make the transfer matrices for the inj orbit here !!! ! see make_mat6 ! Calculate transfer matrix in forward direction from bumper 36E to ! injection point ! Multiply transfer matrices from 36E (don't include) to Q34E (include) tmat = ring%ele(inj%ix_bumper(3)+1)%mat6 ! if ( any( ring%ele(inj%ix_bumper(3)+1)%vec0 /= 0. ) ) print *,' 0th order nonzero!',ring%ele(i)%vec0 nt=1 do i = inj%ix_bumper(3)+2, inj%ix_injpt nt = nt +1 tmat = matmul ( ring%ele(i)%mat6, tmat ) ! Check zero-th order ! if ( any( ring%ele(i)%vec0 /= 0. ) ) print *,' 0th order nonzero!',ring%ele(i)%vec0 end do ! Invert matrix call mat_inverse ( tmat, tmat ) ! Desired injected orbit position at bumper 36E is stored orbit ! position corrected for injected energy offset ! Get stored beam position at 36E. Use tracking to account from pretzel (rather ! than using transfer matrix) xstored = inj%x36_stored(1) xpstored = inj%x36_stored(2) x36e = xstored + synch_beam%centroid(6) * ring%ele(inj%ix_bumper(3))%a%eta inj_angle = ( tmat(1,6) * synch_beam%centroid(6) + tmat(1,1) * synch_beam%centroid(1) - x36e ) / tmat(1,2) write(6,9000)100*xstored, 1000*xpstored, 100*synch_beam%centroid(6),ring%ele(inj%ix_bumper(3))%a%eta,100*x36e, nt, 1000*inj_angle 9000 format(' Rtn INJECTION_ANGLE calculation:'/ & 5x,' Stored beam position (cm) and angle (mrad) at bumper 36E (cm):',2f10.3,/ & 5x,' Injected beam energy offset (%): ',es9.2,/ & 5x,' Eta value at bumper 36E (cm/%):',f8.3,/ & 5x,' Desired injected beam position at bumper 36E (cm):',f8.3,/ & 5x,' Number of elements between Bump36E and Q34E:',i8,/ & 5x,' Returned injection angle (mrad): ',f8.2) error = .false. end subroutine injection_angle !####################################################### subroutine init_injection(ring, pretzel, inj, synch_beam, fname, use_lrbbi) implicit none type (lat_struct) ring type (coord_struct), allocatable, save :: orbit(:), prz_theory(:) type (coord_struct), allocatable :: diff_co(:) type (coord_struct) delta_ip type (inj_struct) inj type (normal_modes_struct) mode type (synch_beam_struct) synch_beam type (pretzel_struct) pretzel logical use_lrbbi real prz1, prz13, dt, wmag, emag logical scan_ping /.false./ real pingcu, kping, synch_epsx, synch_epsy integer ntping real pb26e, pb28e, pb36e real synch_centroid(6), synch_twiss(6) real inj_angle integer ntestturns, nangsteps real angmin, angstep character*40 pingfname real pingchoice, bestangle, besteff real kping_in, bestangle_in real xp36e, xinj, x34e, x36e, beta34e, beta36e integer i_dim integer i integer lun character fname*40 logical write_inj_prz_orbit, write_mag_prz_orbit, write_inj_orbit logical write_pulsed_bump_orbits, use_default_centroid, use_default_twiss real inxrms_ws,indelta_ws,inxrms_wpi,inxprms_wpi ! Injection envelope params data i_dim / 6 / namelist / inj_input / prz1, prz13, wmag, emag, scan_ping, pingcu, kping, dt, ntping, & pb26e, pb28e, pb36e, & write_inj_prz_orbit, write_mag_prz_orbit, write_inj_orbit, & write_pulsed_bump_orbits, use_default_centroid, synch_centroid, & use_default_twiss, synch_twiss, synch_epsx, synch_epsy, & inxrms_ws, indelta_ws, inxrms_wpi, inxprms_wpi, & ntestturns, angmin, angstep, nangsteps, pingfname ! allocate storage call reallocate_coord (orbit, ring%n_ele_max) call reallocate_coord (prz_theory, ring%n_ele_max) call reallocate_coord (diff_co, ring%n_ele_max) write(6,1000)fname 1000 format('Routine INIT_INJECTION opening injection input file ',a40) ! ! Initialize injection parameters prz1=0. prz13=0. wmag=0. emag=0. scan_ping=.false. pingcu=0. kping=0. dt=0. ntping=1 pb26e=0. pb28e=0. pb36e=0. lun = lunget() open(lun, file=fname, status = 'old') read(lun, nml=inj_input) close(lun) ! load injparam structure with envelope calculation parameters inj%xrms_ws = inxrms_ws inj%delta_ws = indelta_ws inj%xrms_wpi = inxrms_wpi inj%xprms_wpi = inxprms_wpi ! print *,' Injection envelope input params',inxrms_ws,indelta_ws,inxrms_wpi,inxprms_wpi ! Get theory closed orbit and twiss parameters ring%param%particle = positron$ ! Set to positrons for setup ! Turn BBI elements off for positrons call lrbbi_ele_setup ( ring, .false., .false. ) ! Update tracking - added 23 July 2004 in order ! to get the new closed orbit without bbi after turning the bbi off call tracking_update(ring, orbit, 4, diff_co) call twiss_at_start (ring) if ( .not. ring%param%stable ) print *,' Rtn INIT_INJECTION: Ring transfer matrix is unstable !!' !!!!!!!! call twiss_propagate_all (ring) prz_theory(0)%vec(:) = 0. ! call closed_orbit_at_start (ring, prz_theory(0), 4,.true.) ! call track_all (ring, prz_theory) call closed_orbit_calc (ring, prz_theory, 4 ) call lat_make_mat6( ring, -1, prz_theory ) call twiss_at_start (ring) call twiss_propagate_all (ring) call init_pretzel( ring, pretzel, prz_theory ) ! call init_mag_pretzel( ring, pretzel, prz_theory ) prz1 = -prz1 ! Want a negative crossing angle at IP for e+ ! 31.3.2003 JAC Change sign of PRZ13 so as to agree with CESRV prz13 = -prz13 ! make the injection pretzel with IP separation, ! iff manual setting of PRZ1 is nonzero or PRZ13 is nonzero 9 July 2003 JAC if ( prz1 .ne. 0 .or. prz13 .ne. 0) then ! Note that using a custom pretzel setting requires recalculating the BBI, ! which is necessarily an iterative procedure. print *,' Routine INIT_INJECTION using custom pretzel. prz1, prz13 =',prz1, prz13 call make_inj_pretzel( ring, pretzel, prz1, prz13 ) else print *,' Routine INIT_INJECTION using design pretzel' endif ! make the magnetic pretzel ! call make_mag_pretzel( ring, pretzel, wmag, emag ) ! Don't write pretzel-off orbits if LRBBI switched on 24 May 2004 JAC !*********************************************************************** if ( .not. use_lrbbi ) then ! Now turn off the pretzel call set_pretzel( 0, ring, pretzel ) ! turn pretzel off if( write_inj_prz_orbit ) then ring%param%particle = electron$ write(fname,1100) 1100 format('dat/inj_prz_orbit.dat') call write_orbit( ring, fname ) ring%param%particle = positron$ endif if( write_mag_prz_orbit ) then call set_pretzel( 0, ring, pretzel ) ! turn off for orbit fname = 'dat/mag_prz_orbit.dat' call write_orbit( ring, fname ) call set_pretzel( 1, ring, pretzel ) ! turn back on endif call set_pretzel( 1, ring, pretzel ) ! turn back on endif !********************************************************************** ! Get the final orbit and twiss parameters print *,' Rtn INIT_INJECTION getting e- stored orbit for injection' ! Calculate beam-beam separations at all elements (delta_ip) ! redefine BBI strengths, and print out separation at IP ! call beambeam_separation ( ring, delta_ip, 4 ) ! Switch to electrons and calculate orbit ring%param%particle = electron$ ! Turn BBI elements on for electrons call lrbbi_ele_setup ( ring, .true., .true. ) ! New simplified routine added 7/2/04 by DTK call tracking_update(ring, orbit, 4, diff_co) fname = 'dat/inj_orbit.dat' call write_orbit( ring, fname ) write(6,6000) 1000*orbit(0)%vec(1), 1000*orbit(0)%vec(2) 6000 format(' Rtn INIT_INJECTION - At the IP, horizontal e- position = ',f6.2, & ' mm angle = ',f6.3,' mrad') ! Finally, set up the injection elements print * print *,'***Pulsed_element_setup***' ! 31.3.2001 JAC Put signs so as to print out the same PRZ values as in setup file print 50,'prz1=',-prz1,' prz13=',-prz13 50 format(1x,a,f8.1,a,f8.1,a,f6.3) print *,' pb26e, pb28e, pb36e=',pb26e,pb28e,pb36e print *,'pingcu=',pingcu,' kping=',kping,' (mrad)' print *,'ntping=',ntping,' dt=',dt ! print *,'emag=',emag,'wmag=',wmag call pulsed_element_setup( ring, orbit, inj, synch_beam, & pb26e, pb28e, pb36e, dt, ntping, kping, pingcu) call twiss_at_start( ring ) orbit(0)%vec(:) = 0. !call closed_orbit_at_start (ring, orbit(0), 4,.true.) !call track_all (ring, orbit) call closed_orbit_calc (ring, orbit, 4 ) print *,' orbit at injpt after setting bump',ring%ele(inj%ix_injpt)%s,orbit(inj%ix_injpt)%vec(1) ! Load mode structure to pass CESR bunch length and energy spread to SYNCH_BEAM_SETUP ! 11 March 2015 jac call radiation_integrals( ring, orbit, mode ) print *,' ***Synch_Beam_Setup***' call synch_beam_setup( ring, inj, mode, synch_beam, use_default_centroid, & synch_centroid, use_default_twiss, & synch_twiss, synch_epsx, synch_epsy, & ntestturns, angmin, angstep, nangsteps) call twiss_at_start( ring ) orbit(0)%vec(:) = 0. !call closed_orbit_at_start (ring, orbit(0), 4,.true.) !call track_all (ring, orbit) call closed_orbit_calc (ring, orbit, 4 ) print *,' orbit at injpt after synch setup',ring%ele(inj%ix_injpt)%s,orbit(inj%ix_injpt)%vec(1) ! Now that we can calculate injection efficiencies, ! we can optimize the pinger !=========================================== ! Pinger Scan. ! The bump settings are independent of the pinger, ! but the optimal pinger kick depends on the ! bump settings. ! ! If scan_ping is true, ! determine the ping value by scanning ! injection efficiency. If not, use ! values from PULSED_ELEMENT_SETUP print *,' **** Pinger Scan ****' print *,' scan_ping, ping scan file=',scan_ping,pingfname if ( scan_ping ) then kping_in = inj%pingkick bestangle_in = synch_beam%centroid(2) call inj_pingscan( ring, inj, synch_beam, & scan_ping, pingfname, besteff, pingchoice, bestangle) if (besteff .gt. 0. ) then print *,' Rtn INIT_INJECTION setting pinger to ',pingchoice,' mrad, i.e. ', & pingchoice/(1000*inj%pingcu2rad),' cu' print *,' and injection angle to ',bestangle inj%pingkick = pingchoice/1000. synch_beam%centroid(2) = bestangle else inj%pingkick = kping_in synch_beam%centroid(2) = bestangle_in print *,' Rtn INJ_PINGSCAN fails to find efficient pinger setting.' print *,' Pinger kick and injection angle left unchanged.' print *,' at ', inj%pingkick, 'rad and ', synch_beam%centroid(2),' mrad' endif endif ! Final values used for injection angle and position print *,' Final values used for synch beam injection postion (m) and angle (rad):' print *,' X, Xprime of injected bunch = ', & synch_beam%centroid(1), synch_beam%centroid(2) print *,' Y, Yprime of injected bunch = ', & synch_beam%centroid(3), synch_beam%centroid(4) ! ! Following "A Description of CESR Injection", page, 7, Stu Henderson, 27 October 1995 x34e = inj%xpb_stored(1) x36e = orbit(inj%ix_bumper(3))%vec(1) xinj = synch_beam%centroid(1) - x34e beta34e = ring%ele(inj%ix_injpt)%a%beta beta36e = ring%ele(inj%ix_bumper(3))%a%beta xp36e = sqrt ( xinj**2 / ( beta34e*beta36e ) ) print *,' ----------------------- Pinger Orbit Correction -----------------------' write(6,'(a,2f8.2)')' Stored orbit X position at 34E with bump = ',100*x34e write(6,'(a,4f8.2)')' xinj (cm), x36e (cm), beta34e (m), beta36e (m) = ',100*xinj,100*x36e,beta34e,beta36e write(6,'(a,f8.2)')' Angle between inj and stored beams at 36E, if both beams on axis (mrad): ', -1000*xp36e write(6,'(a,2f8.3)')' Pinger kick (mrad,cu) = ', 1000*inj%pingkick, 1000*inj%pingkick/inj%pingcu2rad write(6,'(a,f8.3)')' Injected orbit oscillation reduction factor = ', 1 + ( inj%pingkick / xp36e ) print *,' ------------------------------------------------------------------------' orbit(0)%vec(:) = 0. ! call closed_orbit_at_start (ring, orbit(0), 4,.true.) call closed_orbit_calc (ring, orbit, 4 ) print '(a15,4e12.4)','After inj_pingscan: Closed orbit ', orbit(0)%vec(1:4) ! Added 7/19/04 by DTK call tracking_update(ring, orbit, 4, diff_co) if( write_pulsed_bump_orbits ) then call write_pb_orbits( ring, pretzel, inj ) endif call tracking_update(ring, orbit, 4, diff_co) end subroutine init_injection !####################################################### subroutine write_orbit( ring, fname ) implicit none integer i, j, lun character(*) fname type (lat_struct) ring type (coord_struct), allocatable, save :: orbit(:) type (coord_struct), allocatable :: diff_co(:) ! allocate storage call reallocate_coord (orbit, ring%n_ele_max) call reallocate_coord (diff_co, ring%n_ele_max) ! New simplified routine added 7/2/04 by DTK call tracking_update(ring, orbit, 4, diff_co) lun = lunget() print *,' WRITE_ORBIT writing filename = ',fname open(lun, file=fname, status = 'unknown' ) do i=0,ring%n_ele_track if (orbit(i)%state /= alive$) exit write(lun, 100) ring%ele(i)%s, (orbit(i)%vec(j),j=1,6) enddo 100 format(2x, f9.4, 6(1x,e12.5) ) close(lun) end subroutine write_orbit !####################################################### subroutine write_pb_orbits( ring, pretzel, inj ) implicit none integer i, j, k, n, lun1, lun2, t_state type (lat_struct) ring type (pretzel_struct) pretzel type (inj_struct) inj type (coord_struct), allocatable, save :: orbit(:) character*40 fname ! allocate storage call reallocate_coord (orbit, ring%n_ele_max) do k = 1,2 if( k == 1 ) then ! First do electron PB orbits with pretzel lun1 = lunget() fname='dat/ele_pb_prz_orbits.dat' open(lun1, file=fname, status = 'unknown' ) print *,' WRITE_PB_ORBITS writing file ',fname lun2 = lunget() fname = 'dat/stored_ele_at_injpt.dat' open(lun2, file=fname, status = 'unknown' ) print *,' WRITE_PB_ORBITS writing file ',fname else if( k == 2 ) then ! do electron without pretzel lun1 = lunget() call lrbbi_ele_setup ( ring, .false., .false. ) call set_pretzel( 0, ring, pretzel ) fname = 'dat/ele_pb_noprz_orbits.dat' open(lun1, file=fname, status = 'unknown' ) print *,' WRITE_PB_ORBITS writing file ', fname endif ! ! orbit(0)%vec(:) = 0. ! call closed_orbit_at_start (ring, orbit(0), 4,.true.) call twiss_at_start(ring) orbit(0)%vec(:) = 0. call closed_orbit_calc(ring, orbit, 4) print *,' WRITE_PB_ORBITS: closed orbit at IP',orbit(0)%vec(1:4) ! We are tracking backwards for the stored beam, so change sign of angles ! 16 April 2003 JAC NO, don't. Convention changed again. ! orbit(0)%vec(2) = -orbit(0)%vec(2) ! x-prime ! orbit(0)%vec(4) = -orbit(0)%vec(4) ! y-prime ! do n=-5,50 n=-5 t_state = moving_forward$ do while ( n < 50 .and. t_state == moving_forward$) call pulsed_bump(ring, n, inj ) ring%param%particle = antiparticle(ring%param%particle) call track_many (ring, orbit, 0, 0, -1, track_state = t_state) ring%param%particle = antiparticle(ring%param%particle) if ( n >-1 .and. n < 3) & write(6,2350)n,ring%ele(inj%ix_injpt)%s,100*orbit(inj%ix_injpt)%vec(1) 2350 format(' WRITE_PB_ORBITS: Turn',i2,' Orbit at Q34E (s=',f6.2,') is ',f6.3,' cm') if (t_state /= moving_forward$) then print *,' Rtn WRITE_PB_ORBITS: orbit lost while writing file',fname print *,' on turn',n,' at ',t_state do i= 0 , t_state - 1 orbit (i)%vec(:) = 0. enddo endif do i=0, ring%n_ele_track write(lun1,200) n, ring%ele(i)%s, & 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) 200 format(2x, i6, f10.4, 2(1x,f8.3),2(1x,f10.5),1x,f8.3,1x,f11.6) enddo if( k == 1 ) then write(lun2,200) n, ring%ele(inj%ix_injpt)%s, & 1000.*orbit(inj%ix_injpt)%vec(1),1000.*orbit(inj%ix_injpt)%vec(2), & 1000.*orbit(inj%ix_injpt)%vec(3),1000.*orbit(inj%ix_injpt)%vec(4), & 1000.*orbit(inj%ix_injpt)%vec(5),1000.*orbit(i)%vec(6) endif ! Increment turn number n = n + 1 enddo close(lun1) if( k == 1 ) close(lun2) enddo call set_pretzel( 1, ring, pretzel ) call lrbbi_ele_setup ( ring, .true., .true. ) end subroutine write_pb_orbits !####################################################### !+ ! Subroutine envelope_chrom_calc (ring, delta_e, chrom_x, chrom_y) ! ! Subroutine to calculate the chromaticities by computing the tune change ! when then energy is changed. ! ! Modules Needed: ! ! Input: ! ring -- lat_struct: Ring ! delta_e -- Real(rdef): Delta energy used for the calculation. ! If 0 then default of 1.0e-4 is used. ! ! Output: ! delta_e -- Real(rdef): Set to 1.0e-4 if on input DELTA_E =< 0. ! chrom_x -- Real(rdef): Horizontal chromaticity. ! chrom_y -- Real(rdef): Vertical chromaticity. !- !$Id: envelope_chrom_calc.f90,v 1.3 2002/02/23 20:32:12 dcs Exp $ !$Log: envelope_chrom_calc.f90,v $ !Revision 1.3 2002/02/23 20:32:12 dcs !Double/Single Real toggle added ! !Revision 1.2 2001/09/27 18:31:49 rwh24 !UNIX compatibility updates ! subroutine envelope_chrom_calc (ring, delta_e, chrom_x, chrom_y) implicit none type (lat_struct) ring type (lat_struct), save :: ring2 type (coord_struct), allocatable, save :: coord_(:) integer i, key real(rdef) high_tune_x, high_tune_y, low_tune_x, low_tune_y real(rdef) delta_e, chrom_x, chrom_y ! allocate storage call reallocate_coord (coord_, ring%n_ele_max) if (delta_e <= 0) delta_e = 1.0e-4 ring2 = ring ! lower energy tune ! coord_(0)%z%vel = -delta_e coord_(0)%vec(6) = -delta_e ! call closed_orbit_at_start (ring2, coord_(0), 4, .true.) ! call track_all (ring2, coord_) call closed_orbit_calc (ring2, coord_, 4 ) call lat_make_mat6 (ring2, -1, coord_) call twiss_at_start (ring2) low_tune_x = ring2%a%tune / twopi if (low_tune_x < 0) low_tune_x = 1 + low_tune_x low_tune_y = ring2%b%tune / twopi if (low_tune_y < 0) low_tune_y = 1 + low_tune_y ! higher energy tune ! coord_(0)%z%vel = delta_e coord_(0)%vec(6) = delta_e ! call closed_orbit_at_start (ring2, coord_(0), 4, .true.) ! call track_all (ring2, coord_) call closed_orbit_calc (ring2, coord_, 4 ) call lat_make_mat6 (ring2, -1, coord_) call twiss_at_start (ring2) high_tune_x = ring2%a%tune / twopi if (high_tune_x < 0) high_tune_x = 1 + high_tune_x high_tune_y = ring2%b%tune / twopi if (high_tune_y < 0) high_tune_y = 1 + high_tune_y ! compute the chrom chrom_x = (high_tune_x - low_tune_x) / (2 * delta_e) chrom_y = (high_tune_y - low_tune_y) / (2 * delta_e) end subroutine envelope_chrom_calc !####################################################### subroutine inj_pingscan( ring, inj, synch_beam, & scan_ping, pingfname, & besteff, pingchoice, bestangle) ! ! Rtn INJ_PINGSCAN scans ping values and returns optimal ! ping kick value PINGCHOICE (mrad) and injection angle ! BESTANGLE (mrad). The calculated injection efficiency ! BESTEFF is also returned. If no injection found, it ! is negative. ! implicit none type (lat_struct) ring type (coord_struct) delta_ip type (coord_struct), allocatable :: inj_orbit(:) type (coord_struct), allocatable, save :: orbit(:) type (coord_struct), allocatable, save :: closed_orbit(:) type (coord_struct), allocatable :: diff_co(:) type (inj_struct) inj type (synch_beam_struct) synch_beam type (normal_modes_struct) mode logical scan_ping real(rdef), allocatable :: dk1(:) logical ok real sigmax real xclsigma integer i,ii, 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 betax, etax real eff, epsx, epsy, dele logical lostit real pingval,kping character*40 pingfname integer nt,ntmax integer ixclmin, iangle real xclsig, xclsigmin real xclear, xclmin, xclsigt, xclsigtmin, dx integer iele integer itmin integer nangsteps real angmin,angstep,inj_angle,angmax logical error real besteff, pingchoice, bestangle, bestclear integer i_dim, t_state, track_state logical damping_on, fluctuations_on character datfname*40 character fname*40 real wmax integer imax ! Namelist variables integer nturns, n_part, style real pingmin,pingmax,pingstep, nsig logical use_lrbbi, use6D namelist / scan_input / pingmin,pingmax,pingstep, & use_lrbbi, & nturns, nsig, n_part, style, use6D ! allocate storage call reallocate_coord (inj_orbit, ring%n_ele_max) call reallocate_coord (orbit, ring%n_ele_max) call reallocate_coord (closed_orbit, ring%n_ele_max) call reallocate_coord (diff_co, ring%n_ele_max) ! Read namelist file lun = lunget() open(lun,file=pingfname, status= 'old') read(lun, nml=scan_input) close(lun) write(6,1000)pingmin, pingmax, pingstep, use_lrbbi, nturns, nsig, n_part, style, use6D 1000 format(' ==> Program INJ_PINGSCAN called:',/, & ' PINGMIN is: ',f6.3/, & ' PINGMAX is: ',f6.3/, & ' PINGSTEP is: ',f6.3/, & ' USE_LRBBI is: ', l/, & ' NTURNS is: ',i6/, & ' NSIG is: ',f6.3/, & ' N_PART is: ',i6/, & ' STYLE is: ',i6/, & ' 6D Tracking: ',l) lun = lunget() fname = 'dat/inj_pingscan.dat' open(lun, file=datfname, status = 'unknown') lun2 = lunget() fname = 'dat/ping_lost_array.dat' open(lun2, file=datfname, status = 'unknown') ! Make sure pulsed bump and pinger are off ring%ele(inj%ix_bumper(1))%value(hkick$) = 0. ring%ele(inj%ix_bumper(2))%value(hkick$) = 0. ring%ele(inj%ix_bumper(3))%value(hkick$) = 0. do i=1,3 call lat_make_mat6(ring, inj%ix_bumper(i)) enddo call twiss_at_start( ring ) closed_orbit(0)%vec(:) = 0. call closed_orbit_calc(ring, closed_orbit, 4) print *,' Orbit at injpt after entering INJ_PINGSCAN:',closed_orbit(inj%ix_injpt)%vec(1) ! Set up 4-D or 6-D tracking !--------------------------- if ( use6D ) then ! 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$) call set_z_tune(ring) i_dim = 6 print *,' Rtn INJ_PINGSCAN using 6-D tracking with Qz = ',ring%z%tune else i_dim = 4 print *,' Rtn INJ_PINGSCAN using 4-D tracking' endif ! Set up BBI !----------- if( use_lrbbi ) then print *,' Rtn INJ_PINGSCAN turning on BBI. ring%param%n_part= ',ring%param%n_part call lrbbi_ele_setup (ring, .true., .true. ) call tracking_update(ring, closed_orbit, i_dim, diff_co) else call tracking_update(ring, closed_orbit, i_dim, diff_co) endif print *,' Orbit at injpt after setting up BBI and tracking type:',closed_orbit(inj%ix_injpt)%vec(1) ! ! Print out tunes print * print *,' Rtn INJ_PINGSCAN: 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) print '(a33,4e12.4)',' Before ping scan: Closed orbit ', closed_orbit(0)%vec(1:4) ! Get emittance and energy spread ! Make sure this is the electron orbit call radiation_integrals( ring, closed_orbit, mode ) epsx = mode%a%emittance dele = mode%sige_e print *,' Rtn INJ_PINGSCAN finds epsx,dele=',epsx,dele ! Added 7/19/04 by DTK call tracking_update(ring, closed_orbit, i_dim, diff_co) print *,' After radiation_integrals: closed orbit at IP',closed_orbit(0)%vec(1:4) ! Enable aperture limits for tracking bmad_com%aperture_limit_on = .true. ! Set range over which injection phase space will be generated ! if type=2 nsig = abs(nsig) besteff = 0. pingchoice = 0. bestangle = 0. bestclear = 0. ! Set ring particle species to positron for tracking electrons backward in TRACK_MANY ring%param%particle = positron$ pingval = pingmin - pingstep do pingval = pingval + pingstep if (pingval > pingmax) exit ! Convert to rad kping=pingval/1000. inj%pingkick=kping ! Test clearance of stored electron orbit for several turns following ! the pinger here. Skip this pinger value if clearance criterion ! violated at any element. global_com%exit_on_error = .false. ntmax=5 nt=-1 xclsigmin = 100. xclmin = 0.045 ixclmin = -1 itmin = -1 t_state = moving_forward$ ! Track until centroid hits wall and then no further do while( nt < ntmax .and. t_state == moving_forward$) ! Increment turn number nt = nt+1 call pulsed_bump(ring, nt, inj ) ! track_many uses orbit as input for the starting point ! of the tracking, which is in this case element 0 call init_coord (orbit(0), closed_orbit(0)%vec, ring%ele(0), & downstream_end$, electron$, -1) call track_many (ring, orbit, 0, 0, -1, track_state = t_state) ! write(6,2300)nt,ring%ele(inj%ix_injpt)%s,100*orbit(inj%ix_injpt)%vec(1) 2300 format(' Turn',i2,' Orbit at inj point (s=',f6.2,') is ',f6.3,' cm') ! Now test clearance for this turn orbit ! Stop calculating worst aperture if centroid hits wall. if (t_state /= moving_forward$) then iele = t_state if (orbit(iele)%location == exit_end$ .or. iele .eq. 0) then xclear = ring%ele(iele)%value(x1_limit$) - abs(orbit(iele)%vec(1)) else xclear = ring%ele(iele)%value(x1_limit$) - abs(orbit(iele-1)%vec(1)) endif betax = ring%ele(iele)%a%beta etax = ring%ele(iele)%a%eta sigmax = sqrt( epsx * betax + ( dele * etax )**2 ) xclsigtmin = 100. if ( sigmax .gt. 0. ) xclsigtmin = xclear / sigmax else ! No aperture violation ! Loop over elements to determine narrowest aperture on this turn xclsigtmin = 100. do i=ring%n_ele_track,0,-1 if (pingval.eq.-1.1 .and. nt.eq.0 ) then print *,ring%ele(i)%name,ring%ele(i)%s, & 100*ring%ele(i)%value(x1_limit$), 100*orbit(i)%vec(1) endif if ( ring%ele(i)%value(x1_limit$) .gt. 0. ) then ! ! Clearance logic: if this element has aperture limits defined ! at the entrance end, calculate the clearance with the orbit ! value of the preceding element ! Note that $exit_end means high-s end even when tracking backward ! The orbit value is always given for the high-s end if ( ring%ele(i)%aperture_at .eq. exit_end$ ) then dx = ring%ele(i)%value(x1_limit$) - abs(orbit(i)%vec(1)) elseif ( ring%ele(i)%aperture_at .eq. entrance_end$ & .and. i .ne. 0 ) then dx = ring%ele(i)%value(x1_limit$) - abs(orbit(i-1)%vec(1)) elseif ( ring%ele(i)%aperture_at .eq. both_ends$ & .and. i .ne. 0 ) then dx = min ( & ring%ele(i)%value(x1_limit$) - abs(orbit(i-1)%vec(1)), & ring%ele(i)%value(x1_limit$) - abs(orbit(i)%vec(1)) & ) else dx=100. print *, ' Rtn INJ_PINGSCAN Warning! No aperture end type found for element ',ring%ele(i)%name endif ! betax = ring%ele(i)%a%beta etax = ring%ele(i)%a%eta ! print *,' betax, etax = ',ring%ele(i)%name,betax,etax sigmax = sqrt( epsx * betax + ( dele * etax )**2 ) xclsigt = 100. if ( sigmax .gt. 0. ) xclsigt = dx / sigmax ! if ( i .eq. 748 ) print *,' D318: dx,xclsigt=',dx,xclsigt !If minimum clearance found for this turn, store it if ( xclsigt .lt. xclsigtmin ) then xclsigtmin = xclsigt xclear = dx iele = i endif endif !end of existence test for aperture limit for this element enddo ! end loop over elements endif ! end of lost test ! Test if this is the narrowest aperture found yet. ! If negative, it will be, since we exit the turns loop on aperture violations ! If positive, minimum was already found in element loop, so put ".le." here ! print *,'turn,xclsigmin,xclear,xclsigtmin,element',nt,xclsigmin,xclear,xclsigtmin,ring%ele(iele)%name if ( xclsigtmin .lt. xclsigmin ) then ixclmin = iele xclmin = xclear xclsigmin = xclsigtmin sigmax = xclear/xclsigmin itmin = nt ! if (nt.eq.0) & write(6,1800)nt,ring%ele(iele)%name,ring%ele(iele)%s,& xclsigtmin,100*xclear,1000*sigmax 1800 format(6x,' Aperture limit on turn',i2,' at ',a10, & '(s=',f8.3,'): ',f6.1,' sigma, ',f6.2, & ' cm. Beam size:',f5.1,' mm') endif enddo ! End of loop over turns write ( 6, 2100 ) pingval, int(inj%pingturn), ntmax, xclsigmin, & 100*xclmin, & ring%ele(ixclmin)%name, itmin 2100 format('Rtn INJ_PINGSCAN: For a ping kick of',f6.2,' mrad on turn',i2, & ' the minimum wall clearance in first',i2, & ' turns for the stored beam is ', & f5.1,' sigma (',f6.2,' cm) at ',a10,' on turn',i2) ! Scan to find optimum injection angle ! angmin = -0.0049 angmin = -0.0040 angstep = 0.0001 nangsteps = 20 ! nangsteps = 1 angmax = angmin + nangsteps * angstep inj_angle = angmin - angstep do inj_angle = inj_angle + angstep if (inj_angle > angmax) exit synch_beam%centroid(2) = inj_angle ! Calculate injection efficiency ! do i=0,ring%n_ele_track lost_array(i) = 0. enddo totweight = 0.0 injweight = 0.0 do i=1,n_part call inject_particle( ring, .true., 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 ! End loop over particles eff = injweight/totweight ! Choose stored-beam wall clearance criterion here if ( eff .gt. besteff .and. xclsigmin .gt. 2.5 ) then besteff = eff pingchoice = pingval bestangle = inj_angle bestclear = xclsigmin endif wmax = -1 do i = 1, ring%n_ele_track if ( lost_array ( i ) .gt. wmax ) then imax = i wmax = lost_array ( i ) endif enddo ! Write out efficiency scan data write(lun,300) pingval, eff write(6,3100) 1000*inj_angle, totweight, injweight, 100*eff 3100 format(6x,' Rtn INJ_PINGSCAN: Inj_angle=',f6.2,' mrad', & ' Totwt=',e9.2,' Injwt=',e9.2, & ' Efficiency:',f6.1,'%') if ( wmax .gt. 0 )then write(6,3200) 100*wmax/totweight, ring%ele(imax)%name 3200 format(' Largest loss was',f6.1,'% at ',a11) else print *,' ' endif 300 format(1x,f7.4,1x,f7.4,1x,f7.4) ! Write out the lost array here do ii=0,ring%n_ele_track lost_array(ii) = lost_array(ii)/totweight write(lun2, 301) ii, ring%ele(ii)%s, lost_array(ii) 301 format(1x,i5, 1x, f10.4, 1x, f10.6 ) ! enddo enddo ! End loop over injection angle enddo ! End loop over ping values ! Turn off the pulsed bump ! ring%ele(inj%ix_bumper(1))%value(hkick$) = 0. ! ring%ele(inj%ix_bumper(2))%value(hkick$) = 0. ! ring%ele(inj%ix_bumper(3))%value(hkick$) = 0. ! do i=1,3 ! call lat_make_mat6(ring, inj%ix_bumper(i), orbit) ! enddo write (6,6000) 100*besteff, pingchoice, pingchoice/inj%pingcu2rad, 1000*bestangle, bestclear 6000 format ( ' Exiting Rtn INJ_PINGSCAN.'/ & ' Best efficiency of ',f6.1, & '% found for ping kick of ',f6.2,' mrad (',es11.4,' cu)', & ' Best injection angle is ',f6.2,' mrad', & ' Stored beam clearance is >',f6.2,' sigma') call tracking_update(ring, closed_orbit, i_dim, diff_co) print *,' After ping scan: closed orbit at IP',closed_orbit(0)%vec(1:4) ! Make sure pulsed bump and pinger are off ring%ele(inj%ix_bumper(1))%value(hkick$) = 0. ring%ele(inj%ix_bumper(2))%value(hkick$) = 0. ring%ele(inj%ix_bumper(3))%value(hkick$) = 0. do i=1,3 call lat_make_mat6(ring, inj%ix_bumper(i)) enddo ring%ele(inj%ix_pinger)%value(hkick$) = 0. call lat_make_mat6(ring, inj%ix_pinger) call twiss_at_start(ring) closed_orbit(0)%vec(:) = 0. call closed_orbit_calc(ring, closed_orbit, i_dim) print *,' After pb off: closed orbit at IP',closed_orbit(0)%vec(1:4) end subroutine inj_pingscan !####################################################### ! Sequence of calls needed to store electron tune values ! in ring structure following a change of ! the strong beam current value ! ! DTK July 2004 ! ! Modified to include beambeam_separation JAC 26 July 2004 subroutine update_tunes(ring, current, stable) implicit none type (lat_struct) ring type (coord_struct), allocatable :: orbit(:) type (coord_struct) delta_ip real(rdef) current integer status logical stable,positron, err_flag call reallocate_coord(orbit,ring%n_ele_max) ! Store particle type for later restoration if(ring%param%particle == positron$) then positron = .true. else positron = .false. endif ring%param%n_part = current*0.001 *(ring%param%total_length/c_light)/e_charge ! The new set of calls call clear_lat_1turn_mats (ring) ! orbit(0)%vec = 0. ! global_com%exit_on_error = .false. ! call closed_orbit_calc (ring, orbit, 4) ring%param%particle = electron$ ! Always do beambeam_separation with electrons call beambeam_separation ( ring, delta_ip, 4 ) ! Recalculate closed orbit for electrons after ! bbi kicks redefined by beambeam_separation orbit(0)%vec(:) = 0. call closed_orbit_calc (ring, orbit, 4, err_flag = err_flag) if (err_flag) stable = .false. call lat_make_mat6(ring,-1,orbit) call twiss_at_start(ring, status) if (status /= ok$) stable = .false. call twiss_propagate_all( ring ) if (.not.ring%param%stable) stable = .false. ! Restore particle type as before call to TRACKING_UPDATE if(positron) ring%param%particle = positron$ end subroutine update_tunes !####################################################### subroutine test_bbi(ring, current) ! Routine to test effect of strong beam current ! and PRZ13 on tunes ! ! DTK July 2004 ! ! Removed calls to beambeam_separation after modifying ! update_tunes to use beambeam_separation JAC 26 July 2004 implicit none type (lat_struct) ring type (lat_struct), save :: ring_in, ring_out type (pretzel_struct) pretzel type (coord_struct), allocatable, save :: prz_theory(:), co_(:) type (coord_struct) delta_ip ! integer particle, i_train, j_car, n_trains_tot, n_cars, slices, i logical stable real(rdef) current !, coupling_sb print *,' === Entering Routine TEST_BBI ===' print *,' Show tunes wrt strong beam current and PRZ13 (4D tracking)' call reallocate_coord (prz_theory, ring%n_ele_max) call reallocate_coord (co_, ring%n_ele_max) global_com%exit_on_error = .false. ! First type out default tunes write(6,1000)ring%a%tune/twopi,ring%b%tune/twopi 1000 format('Tunes upon entering TEST_BBI: Qx = ',f8.5, ' Qy = ',f8.5) ! All results for electrons ring%param%particle = electron$ ! Initialize pretzel stuff call twiss_at_start ( ring ) prz_theory(0)%vec(:) = 0. call closed_orbit_calc (ring, prz_theory, 4) call init_pretzel( ring, pretzel, prz_theory ) ! Tunes with IP LRBBI call lrbbi_ele_setup(ring, .true., .false.) stable = .true. call update_tunes(ring, current, stable) if(stable) then write(6,1100)ring%a%tune/twopi,ring%b%tune/twopi 1100 format('Tunes with IP BBI and without LRBBI: Qx = ',f8.5, ' Qy = ',f8.5) else print *, 'Not stable. IP BBI tunes not computed.' endif ! ring_sep = ring ! call beambeam_separation ( ring_sep, delta_ip, 4 ) stable = .true. call update_tunes(ring, current*3./4., stable) if(stable) then write(6,1150)ring%a%tune/twopi,ring%b%tune/twopi 1150 format('Tunes with IP BBI and 3/4 current: Qx = ',f8.5, ' Qy = ',f8.5) else print *, 'Not stable. IP BBI with 3/4 current tunes not computed.' endif ! ring_sep = ring ! call beambeam_separation ( ring_sep, delta_ip, 4 ) stable = .true. call update_tunes(ring, current/2, stable) if(stable) then write(6,1200)ring%a%tune/twopi,ring%b%tune/twopi 1200 format('Tunes with IP BBI and half current: Qx = ',f8.5, ' Qy = ',f8.5) else print *, 'Not stable. IP BBI with half current tunes not computed.' endif ! ring_sep = ring ! call beambeam_separation ( ring_sep, delta_ip, 4 ) stable = .true. call update_tunes(ring, current/4., stable) if(stable) then write(6,1250)ring%a%tune/twopi,ring%b%tune/twopi 1250 format('Tunes with IP BBI and 1/4 current: Qx = ',f8.5, ' Qy = ',f8.5) else print *, 'Not stable. IP BBI with 1/4 current tunes not computed.' endif ! ring_sep = ring ! call beambeam_separation ( ring_sep, delta_ip, 4 ) ! Tunes with full current and pretzel 13 at 1000cu call make_inj_pretzel( ring, pretzel, 0., 1000. ) call set_pretzel( 1, ring, pretzel ) stable = .true. call update_tunes(ring, current, stable) if(stable) then write(6,1300)ring%a%tune/twopi,ring%b%tune/twopi 1300 format('Tunes with IP BBI and pretzel 13 at 1000cu: Qx = ',f8.5, ' Qy = ',f8.5) else print *, 'Not stable. IP BBI with pretzel 13 tunes not computed.' endif ! ring_sep = ring ! call beambeam_separation ( ring_sep, delta_ip, 4 ) stable = .true. call update_tunes(ring, current*3./4., stable) if(stable) then write(6,1350)ring%a%tune/twopi,ring%b%tune/twopi 1350 format('Tunes with IP BBI, pretzel 13 at 1000cu and 3/4 current: Qx = ',f8.5, ' Qy = ',f8.5) else print *, 'Not stable. IP BBI with pretzel 13 and 3/4 current tunes not computed.' endif ! ring_sep = ring ! call beambeam_separation ( ring_sep, delta_ip, 4 ) stable = .true. call update_tunes(ring, current/2, stable) if(stable) then write(6,1400)ring%a%tune/twopi,ring%b%tune/twopi 1400 format('Tunes with IP BBI, pretzel 13 at 1000cu and half current: Qx = ',f8.5, ' Qy = ',f8.5) else print *, 'Not stable. IP BBI with pretzel 13 and half current tunes not computed.' endif ! ring_sep = ring ! call beambeam_separation ( ring_sep, delta_ip, 4 ) stable = .true. call update_tunes(ring, current*1./4., stable) if(stable) then write(6,1450)ring%a%tune/twopi,ring%b%tune/twopi 1450 format('Tunes with IP BBI, pretzel 13 at 1000cu and 1/4 current: Qx = ',f8.5, ' Qy = ',f8.5) else print *, 'Not stable. IP BBI with pretzel 13 and 1/4 current tunes not computed.' endif ! ring_sep = ring ! call beambeam_separation ( ring_sep, delta_ip, 4 ) ! Tunes with non-IP BBI ! We want to set the prz13 back to 0 print *,' Rtn TEST_BBI setting PRZ13 back to 0 ...' call make_inj_pretzel( ring, pretzel, 0., -1000. ) ! call set_pretzel( 0, ring, pretzel ) call lrbbi_ele_setup(ring, .false., .true.) stable = .true. call update_tunes(ring, current, stable) if(stable) then write(6,1500)ring%a%tune/twopi,ring%b%tune/twopi 1500 format('Tunes with non-IP BBI and prz13=0: Qx = ',f8.5, ' Qy = ',f8.5) else print *, 'Not stable. Non-IP BBI tunes not computed.' endif ! Electron tunes with full BBI call lrbbi_ele_setup(ring, .true., .true.) stable = .true. print *,' particle before last update_tunes in test_bbi:',ring%param%particle call update_tunes(ring, current, stable) ! Get closed orbit and tunes for electrons ! ring%param%particle = electron$ ! call twiss_at_start(ring) ! co_(0)%vec(:) = 0. ! call closed_orbit_at_start(ring, co_(0), 4, .true. ) ! call track_all (ring, co_) ! call lat_make_mat6(ring,-1,co_) ! call twiss_at_start(ring) ! call twiss_propagate_all(ring) if(stable) then write(6,1600)ring%a%tune/twopi,ring%b%tune/twopi 1600 format('Electron tunes with full BBI and prz13=0: Qx = ',f8.5, ' Qy = ',f8.5) else print *, 'Not stable. Full BBI tunes not computed.' endif print *,' === Exiting Routine TEST_BBI ===' end subroutine test_bbi !####################################################### subroutine track_nt(iflag, ring, inj, itrack, iturn, orbit, ifirst, npts, idir) ! ! Book, fill and close ntuple with track data ! ! Input arguments: ! IFLAG = 1 to close ntuple, ignored otherwise ! RING Ring data ! ITRACK Track ID number ! iTURN Turn number ! ORB Structure containing the track data ! IFIRST Index of first element with track data ! NPTS Total number of points on the track ! IDIR Direction of the track (+ for e+, - for e-) ! use bmad use injparam_mod implicit none type (lat_struct) ring type (coord_struct) orbit(0:ring%n_ele_max) type (inj_struct) inj type (synch_beam_struct) synch_beam ! Input arguments integer iflag integer iturn integer itrack, ifirst, npts, idir ! integer n, inc ! integer, save :: ntlun integer istat,icycle character ntname*40 ! integer i,k,id0,nvar character*12 htitle character*8 tags(11) real vv(11) ! data nvar/11/ data htitle/' Tracks '/ data tags/' id ','injxpos ','injxang','injde ', & ' turn ','element ', & ' s ',' x ',' xang ',' y ',' yang '/ ! integer ienter data ienter/0/ ! !----------------------------------------------------------------------- ! Allocate storage ! allocate (orbit(0:ring%n_ele_max)) ! ! write(6,1000)iflag, itrack, iturn, ifirst, npts, idir 1000 format(' Rtn TRACK_NT called with iflag, itrack, iturn, ifirst, npts, idir=', & 6i5) ! print *, 'ntlun=',ntlun ! ienter=ienter+1 if(ienter.eq.1)then ! Don't close unbooked ntuple ! Don't book ntuple on closing call if ( iflag .eq. 1 )goto 950 ! Open file for ntuple ntlun=lunget() ntname='ntuples/tracks.ntu' print *,' Rtn TRACK_NT opening ntuple file ',ntname call hrOPEN(ntlun,'NT2',ntname,'N',1024,istat) ! Book ntuple id0=10 call hbookn(id0,htitle,nvar,'//NT2',10000,tags) else ! Close ntuple if(iflag.eq.1)goto 900 endif ! ! Exit with error message if NPTS=0 if ( npts .le. 0 ) then print *,' Rtn TRACK_NT exits on error. Input argument NPTS=',npts goto 950 endif ! vv(1) = itrack vv(2) = orbit(inj%ix_injpt)%vec(1) vv(3) = orbit(inj%ix_injpt)%vec(2) vv(4) = orbit(inj%ix_injpt)%vec(6) vv(5) = iturn ! i=ifirst if ( i == 0 ) i = ring%n_ele_track inc=idir/abs(idir) n = 0 do while ( n < npts ) n=n+1 ! print *,' in loop. inc,i,n=',inc,i,n vv(6) = i vv(7) = ring%ele(i)%s vv(8) = orbit(i)%vec(1) vv(9) = orbit(i)%vec(2) vv(10) = orbit(i)%vec(3) vv(11) = orbit(i)%vec(4) ! ! print *,vv ! call hfn(id0,vv) ! i=mod( i + inc + ring%n_ele_track + 1 , ring%n_ele_track + 1 ) ! enddo ! goto 950 ! 900 continue !================================== ! Write out histograms and ntuples print *,' Rtn TRACK_NT writing out ntuple to file ', ntname call hcdir('//NT2',' ') CALL HROUT(0,ICYCLE,' ') CALL HREND('NT2') close(ntlun) ! 950 continue ! end subroutine track_nt subroutine track_injected_particle( ring, inj, species, inj_orbit, nturns, & writeout, nwrite, lun, np, weight, lostit, track_state) ! ! Track injected particle of specied given by input argument SPECIES ! ! Input arguments: RING ring structure ! INJ inj structure ! INJ_ORBIT injection orbit structure ! NTURNS number of turns to track ! WRITEOUT logical for output to unit LUN ! NWRITE write out every NWRITEth turn ! LUN unit number to write to ! NP particle index number ! WEIGHT weight associated with this particle ! Output arguments ! LOSTIT logical, true if particle lost ! TRACK_STATE status variable from TRACK_MANY ! Set to MOVING_FORWARD$ if not lost, element index where lost if lost use bmad use injparam_mod implicit none ! Input structures type (lat_struct) ring type (inj_struct) inj type (coord_struct) inj_orbit(0:ring%n_ele_max) ! Input arguments integer species, nturns, nwrite, lun, np, track_state real weight logical lostit, writeout ! Internal variables integer n, i, ix_start, j integer ixw integer ilostele real s_lost integer kk integer ring_species_in ix_start = inj%ix_injpt ! start at the injection point, go to IP n = 0 lostit = .false. ! 16 April 2003 JAC Impose a maximum aperture limit here to avoid ! sinh overflow errors from TRACK_MANY which occur with the 6-wiggler lattice ! The errors were avoided by lowering this maximum aperture from 1 m to 10 cm. bmad_com%max_aperture_limit = 0.1 ! Make sure injected particle is actually in the aperture at inj pt. ! if( bmad_com%aperture_limit_on ) then ! if( inj_orbit(inj%ix_injpt)%vec(1) > -0.045 - inj%septthick ) then ! lostit = .true. ! endif ! For backwards tracking in TRACK_MANY, set ring particle species opposite to ! that of tracked particle. Save to return RING with incoming species ring_species_in = ring%param%particle ring%param%particle = antiparticle(species) do while( n < nturns .and. .not. lostit ) call pulsed_bump(ring, n, inj ) if( n == 0 )then ! Expand aperture for selected elements on injection turn 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.10 ring%ele(j)%value(x2_limit$)=0.10 ! print *,' Expanding turn 0 aperture limit to 10 cm for element ',j,ring%ele(j)%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 this 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) call init_coord (inj_orbit(ix_start), inj_orbit(ix_start)%vec, ring%ele(ix_start), & downstream_end$, species, -1) call track_many (ring, inj_orbit, ix_start, 0, -1, track_state = track_state) ! do i=0, ring%n_ele_track ! if ( np == 1 ) write(6,101) n, ring%ele(i)%name, 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 101 format(' Turn ', i3, 2x, a10, f10.4, 2(1x,f8.3),2(1x,f10.5),1x,f8.3,1x,f11.6) ! After injection turn, set aperture near injection point back to ! beampipe wall position 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 ! print *,' Resetting aperture limit to 4.5 cm for element ',j,ring%ele(j)%name endif enddo else ! Track around full ring from element zero when not on injection turn (n=0) call track_many (ring, inj_orbit, 0, 0, -1, track_state = track_state) endif if( track_state /= moving_forward$) then lostit = .true. ! write(6,5000)np, n,ring%param%end_lost_at, ring%param%ix_lost, & ! ring%ele( ring%param%ix_lost )%name, & ! 100*inj_orbit(ring%param%ix_lost)%vec(1), & ! (inj_orbit(ix_start)%vec(kk),kk=1,6) 5000 format('Lost injected electron nr',i5, & ' on turn ',i2,' at end',i3,' of element ',i5,2x,a12,' x=',f6.2,' cm'/, & 5x, ' Injected phase space values:',6(1x,e12.5)) endif ! Collect statistics on undulator transfer functions ! by writing output file of phase space values for all orbits call und_inj_ps (ring, inj, inj_orbit, n, track_state, weight) n = n+1 ! Write out tracked particle coordinates at injection point on selected turns ! to ps_track.dat if( writeout .and. mod(n,nwrite) == 0 ) then ! Choose ring position at which to write out phasespace data ixw=inj%ix_injpt ! ixw=attribute_index ( ele, 'Q08W') ! ixw=72 ! if ( ixw .gt. 0 ) then ! Changed signs of angles and time to adapt to new convention ! 16 April 2003 JAC ! Added weight and element at which the particle was lost, if lost if(lostit)then ilostele = track_state s_lost = ring%ele(ilostele)%s else ilostele = -1 s_lost = -1. endif write(lun,301) np, n-1, inj_orbit(ixw)%vec(1), -inj_orbit(ixw)%vec(2), & inj_orbit(ixw)%vec(3), -inj_orbit(ixw)%vec(4), & -inj_orbit(ixw)%vec(5), inj_orbit(ixw)%vec(6), & weight,ilostele,s_lost 301 format(1x,i5,1x,i6,1x,7(1x,e12.5),i10,1x,f8.3) else print *,' Error in TRACK_INJECTED_PARTICLE: Ring index ixw=',ixw print *,' so phasespace file not written' endif else ! print *,' Phasespace record not written. particle, turn, writeout, mod(n,nwrite=', & ! np, n, writeout, mod(n,nwrite) endif enddo ! Restore incoming ring particle species ring%param%particle = ring_species_in !Delete this when done ! if( lostit ) then ! print *,'Lost at pingval - ',inj%pingkick,' on turn ',(n-1) ! print *,'inj%pingturn - ',inj%pingturn ! endif ! Make sure pulsed bump and pinger are off ! If the orbit did not survive on a turn where the ! bump or pinger were on, they could still be left ! on by PULSED_BUMP at this point. ring%ele(inj%ix_bumper(1))%value(hkick$) = 0. ring%ele(inj%ix_bumper(2))%value(hkick$) = 0. ring%ele(inj%ix_bumper(3))%value(hkick$) = 0. do i=1,3 call lat_make_mat6(ring, inj%ix_bumper(i)) enddo ring%ele(inj%ix_pinger)%value(hkick$) = 0. call lat_make_mat6(ring, inj%ix_pinger) end subroutine track_injected_particle !####################################################### subroutine und_inj_ps (ring, inj, inj_orbit, iturn, track_state, weight) ! Write file of transfer function statistics for wigglers and undulators ! from TRACK_INJECTED_PARTICLE ! ! 6 Feb 2015 jac ! ! Input arguments: RING ring structure ! INJ inj structure containing injection parameters ! INJ_ORBIT injection orbit structure for injection turn ITURN ! TRACK_STATE status variable from TRACK_MANY ! Set to MOVING_FORWARD$ if not lost, element index where lost if lost ! WEIGHT is a probability value from INJECT_PARTICLE for calculating efficiencies use bmad use injparam_mod implicit none ! Input structures type (lat_struct) ring type (inj_struct) inj type (coord_struct) inj_orbit(0:ring%n_ele_max) ! Input arguments integer iturn, track_state real weight ! Internal variables integer ienter /0/ save ienter integer ilimit /10000/ save ilimit integer lun save lun character*18 fname /'dat/und_inj_ps.dat'/ integer ix_inj save ix_inj integer i, ikey real(rp) s_lost ! Gatekeeper with file size circuit breaker if (ienter.lt.ilimit)then ienter = ienter + 1 elseif (ienter.eq.ilimit)then write(*,1000)ilimit 1000 format(' Rtn UND_INJ_PS reaches limit of ',i6, & ' calls. No more output will be written.') ienter = ienter + 1 else goto 999 endif if (ienter.eq.1) then lun = lunget() print *,' Rtn UND_INJ_PS opening output file ',fname open(lun, file=fname, status = 'unknown' ) ix_inj = inj%ix_injpt if (ix_inj.eq.0)then print *,' Rtn UND_INJ_PS finds IX_INJ=0. No output will be written.' close(lun) go to 999 endif else if (ix_inj.eq.0) goto 999 endif ! Write out injection coordinates and incoming and outgoing ! orbit vectors for all wiggler-type elements ! NB: element index i-1 is the outgoing end of element i ! since we are tracking backward do i=0,ring%n_ele_track ikey = ring%ele(i)%key if(track_state.ge.0.and.track_state.le.ring%n_ele_track)then s_lost = ring%ele(track_state)%s else s_lost=-1. endif if (ikey.eq.wiggler$ .and. ring%ele(i)%name(1:3).eq.'CCU') then write(lun,2000)iturn,s_lost,inj_orbit(ix_inj)%vec, & i,inj_orbit(i)%vec, inj_orbit(i-1)%vec, weight, & trim(ring%ele(i)%name) endif enddo 2000 format(i3,1x,f6.2,1x,6(e12.5,1x),i6,1x,6(e12.5,1x),7(e12.5,1x),a) 999 continue end subroutine und_inj_ps !####################################################### real function myran( init ) use nr use nrutil implicit none REAL :: y, secnds integer, save :: idum integer :: ncall, i logical :: init if( init ) then ! initialize and randomize starting point idum = -1 y = ran( idum ) y = secnds(0.0) ! seconds since midnight y=1. ncall = y ! print *,' calling ran ',ncall,' times ' do i=1,ncall ! start off the sequence in a random spot y = ran( idum ) enddo endif myran = ran( idum ) end function end module injlib