module calc_res_mod use bmad use cesr_basic_mod use radiation_mod use sim_utils use srdt_mod implicit none contains !============================================================================= !============================================================================= subroutine calc_h1_h2 (ring, h1, h2, order) type (lat_struct):: ring complex*8, intent(out) :: h1(10) complex*8, intent(out) :: h2(11) integer order type(summation_rdt_struct) srdt_sums integer n_slices_gen_opt integer n_slices_sxt_opt n_slices_gen_opt = 15 n_slices_sxt_opt = 30 call srdt_calc (ring, srdt_sums, order, n_slices_gen_opt, n_slices_sxt_opt) h1(1) = srdt_sums%h11001 h1(2) = srdt_sums%h00111 h1(3) = srdt_sums%h20001 h1(4) = srdt_sums%h00201 h1(5) = srdt_sums%h10002 h1(6) = srdt_sums%h21000 h1(7) = srdt_sums%h30000 h1(8) = srdt_sums%h10110 h1(9) = srdt_sums%h10020 h1(10) = srdt_sums%h10200 h2(:) = 0 if (order .ge. 2) then h2(1) = srdt_sums%h22000 ! dnu_x/dj_x h2(2) = srdt_sums%h11110 ! dnu_x,y/dj_y,x h2(3) = srdt_sums%h00220 ! dnu_y/dj_y h2(4) = srdt_sums%h31000 ! 2Qx h2(5) = srdt_sums%h11200 ! 2Qy h2(6) = srdt_sums%h40000 ! 4Qx h2(7) = srdt_sums%h20020 ! 2Qx-2Qy h2(8) = srdt_sums%h20110 ! 2Qx h2(9) = srdt_sums%h20200 ! 2Qx+2Qy h2(10) = srdt_sums%h00310 ! 2Qy h2(11) = srdt_sums%h00400 ! 4Qy end if end subroutine calc_h1_h2 !============================================================================= !============================================================================= subroutine calc_h(ring,h) type (lat_struct):: ring type(ele_struct) :: ele ! Lie generator terms: break into real and imaginary components real(rp) :: h11001r , h00111r , h20001r , h00201r , h10002r , & h21000r , h30000r , h10110r , h10020r , h10200r , & h11001i , h00111i , h20001i , h00201i , h10002i , & h21000i , h30000i , h10110i , h10020i , h10200i complex*8, intent(out) :: h(10) integer ix, jx real(rp) :: b2l,b3l, b1_grad, b2_grad, b1, b2 integer :: count = 0. h11001r = 0; h00111r = 0; h20001r = 0; h00201r = 0; h10002r = 0 h21000r = 0; h30000r = 0; h10110r = 0; h10020r = 0; h10200r = 0 h11001i = 0; h00111i = 0; h20001i = 0; h00201i = 0; h10002i = 0 h21000i = 0; h30000i = 0; h10110i = 0; h10020i = 0; h10200i = 0 b2l=0.;b3l=0. ! filter down to elements only having nonzero k1/k2/b1/b2 if (ring%ele(0)%ix_pointer .ne. -1) then call filter_eles(ring) endif ! loop over all elements: do ix = 1, ring%n_ele_track if ((ring%ele(ix)%ix_pointer .eq. 1) .or. (ring%ele(ix)%ix_pointer .eq. 2)) then if (match_reg(ring%ele(ix)%name, "SK")) cycle ele = ring%ele(ix) call twiss_at_element (ring%ele(ix), average=ele) ! 'bnl' notation here follows Bengtsson, where n starts at 1 rather than zero-- i.e, b2l(Bengtsson) = b1l(Bmad) b1_grad = value_of_attribute(ele,'B1_GRADIENT') if (b1_grad .eq. real_garbage$) b1_grad = 0. b1 = value_of_attribute(ele,'B1') if (b1 .eq. real_garbage$) b1 = 0. b2_grad = value_of_attribute(ele,'B2_GRADIENT') if (b2_grad .eq. real_garbage$) b2_grad = 0. b2 = value_of_attribute(ele,'B2') if (b2 .eq. real_garbage$) b2 = 0. b2l = (b1_grad+b1) * value_of_attribute(ele,'L') b3l = (b2_grad+b2) * value_of_attribute(ele,'L')/2. h11001r = h11001r + (1/4.)*(b2l-2.*b3l*ele%a%eta) * ele%a%beta h00111r = h00111r + (-1/4.)*(b2l-2.*b3l*ele%a%eta) * ele%b%beta h20001r = h20001r + (1/8.)*(b2l-2.*b3l*ele%a%eta) * ele%a%beta * cos(2.*ele%a%phi) h20001i = h20001i + (1/8.)*(b2l-2.*b3l*ele%a%eta) * ele%a%beta * sin(2.*ele%a%phi) h00201r = h00201r + (-1/8.)*(b2l-2.*b3l*ele%a%eta) * ele%b%beta * cos(2.*ele%b%phi) h00201i = h00201i + (-1/8.)*(b2l-2.*b3l*ele%a%eta) * ele%b%beta * sin(2.*ele%b%phi) h10002r = h10002r + (1/2.)*(b2l-b3l*ele%a%eta) * ele%a%eta * sqrt(ele%a%beta) * cos(ele%a%phi) h10002i = h10002i + (1/2.)*(b2l-b3l*ele%a%eta) * ele%a%eta * sqrt(ele%a%beta) * sin(ele%a%phi) h21000r = h21000r + (-1/8.)*b3l*(ele%a%beta)**(3/2.) * cos(ele%a%phi) h21000i = h21000i + (-1/8.)*b3l*(ele%a%beta)**(3/2.) * sin(ele%a%phi) h30000r = h30000r + (-1/24.)*b3l*(ele%a%beta)**(3/2.) * cos(3.*ele%a%phi) h30000i = h30000i + (-1/24.)*b3l*(ele%a%beta)**(3/2.) * sin(3.*ele%a%phi) h10110r = h10110r + (1/4.)*b3l*(ele%a%beta)**(1/2.) * ele%b%beta * cos(ele%a%phi) h10110i = h10110i + (1/4.)*b3l*(ele%a%beta)**(1/2.) * ele%b%beta * sin(ele%a%phi) h10020r = h10020r + (1/8.)*b3l*(ele%a%beta)**(1/2.) * ele%b%beta * cos(ele%a%phi - 2.*ele%b%phi) h10020i = h10020i + (1/8.)*b3l*(ele%a%beta)**(1/2.) * ele%b%beta * sin(ele%a%phi - 2.*ele%b%phi) h10200r = h10200r + (1/8.)*b3l*(ele%a%beta)**(1/2.) * ele%b%beta * cos(ele%a%phi + 2.*ele%b%phi) h10200i = h10200i + (1/8.)*b3l*(ele%a%beta)**(1/2.) * ele%b%beta * sin(ele%a%phi + 2.*ele%b%phi) endif enddo ! load values into "h" for return: h(1) = cmplx(h11001r,h11001i) h(2) = cmplx(h00111r,h00111i) h(3) = cmplx(h20001r,h20001i) h(4) = cmplx(h00201r,h00201i) h(5) = cmplx(h10002r,h10002i) h(6) = cmplx(h21000r,h21000i) h(7) = cmplx(h30000r,h30000i) h(8) = cmplx(h10110r,h10110i) h(9) = cmplx(h10020r,h10020i) h(10) = cmplx(h10200r,h10200i) end subroutine calc_h !============================================================================= !============================================================================= subroutine calc_xi(ring, xix, xiy) ! BROKEN as of 2016.03.24 or earlier type(lat_struct) :: ring type(ele_struct) :: ele real :: b2l=0.,b3l=0., b1_grad, b1, b2_grad, b2 integer ix real(rp), intent(out) :: xix, xiy xix=0. xiy=0. ! filter down to elements only having nonzero k1/k2/b1/b2 if (ring%ele(0)%ix_pointer .ne. -1) then call filter_eles(ring) endif do ix=1, ring%n_ele_track if ((ring%ele(ix)%ix_pointer .eq. 1) .or. (ring%ele(ix)%ix_pointer .eq. 2)) then if (match_reg(ring%ele(ix)%name, "SK")) cycle ele = ring%ele(ix) call twiss_at_element (ring%ele(ix), average = ele) b1_grad = value_of_attribute(ele,'B1_GRADIENT') if (b1_grad .eq. real_garbage$) b1_grad = 0. b1 = value_of_attribute(ele,'B1') if (b1 .eq. real_garbage$) b1 = 0. b2_grad = value_of_attribute(ele,'B2_GRADIENT') if (b2_grad .eq. real_garbage$) b2_grad = 0. b2 = value_of_attribute(ele,'B2') if (b2 .eq. real_garbage$) b2 = 0. ! 'bnl' notation here follows Bengtsson, where n starts at 1 rather than zero-- i.e, b2l(Bengtsson) = b1l(Bmad) b2l = (b1_grad + b1) * value_of_attribute(ele,'L') b3l = (b2_grad + b2) * value_of_attribute(ele,'L')/2. xix = xix + (-1/(4.*pi))*(b2l - 2.*b3l*ele%a%eta)*ele%a%beta xiy = xiy + (1/(4.*pi))*(b2l - 2.*b3l*ele%a%eta)*ele%b%beta endif enddo ring%a%chrom = xix ring%b%chrom = xiy end subroutine calc_xi !============================================================================= !============================================================================= ! subroutine for calculating dbeta_ddelta and deta_ddelta: subroutine calc_dbeta_ddelta(ring, etax2, betax1, betay1) type(lat_struct) :: ring type(ele_struct) :: ele, ele2 real(rp) :: b2l=0.,b3l=0., b2l2=0.,b3l2=0., b1_grad, b1, b2_grad, b2 integer ix, jx real(rp), allocatable, intent(out) :: etax2(:), betax1(:), betay1(:) if (.not. allocated(etax2)) allocate(etax2(1:ring%n_ele_track)) if (.not. allocated(betax1)) allocate(betax1(1:ring%n_ele_track)) if (.not. allocated(betay1)) allocate(betay1(1:ring%n_ele_track)) ! initialize to zero: etax2(:) = 0. betax1(:) = 0. betay1(:) = 0. if (ring%ele(0)%ix_pointer .ne. -1) then ! have not filtered by k1/k2/b1/b2 yet. do it now. call filter_eles(ring) endif loop1: do ix=1, ring%n_ele_track if (value_of_attribute(ring%ele(ix), 'L') .eq. 0.) cycle loop1 ! throw away markers ! use k1/k2/b1/b2 filter only for speed purposes. This assumes somewhat uniform distribution of elements with these attributes, ! in order to compute a ringwide average etax2, betax1, betay1 if ((ring%ele(ix)%ix_pointer .eq. 1) .or. (ring%ele(ix)%ix_pointer .eq. 2)) then ele = ring%ele(ix) call twiss_at_element (ring%ele(ix), average = ele) etax2(ix) = -ele%a%eta ! initialize loop2: do jx=1, ring%n_ele_track if ((ring%ele(jx)%ix_pointer .eq. 1) .or. (ring%ele(jx)%ix_pointer .eq. 2)) then ! filter by k1/k2/b1/b2 if (match_reg(ring%ele(jx)%name, "SK")) cycle loop2 ele2 = ring%ele(jx) call twiss_at_element (ring%ele(jx), average = ele2) b1_grad = value_of_attribute(ele2,'B1_GRADIENT') if (b1_grad .eq. real_garbage$) b1_grad = 0. b1 = value_of_attribute(ele2,'B1') if (b1 .eq. real_garbage$) b1 = 0. b2_grad = value_of_attribute(ele2,'B2_GRADIENT') if (b2_grad .eq. real_garbage$) b2_grad = 0. b2 = value_of_attribute(ele2,'B2') if (b2 .eq. real_garbage$) b2 = 0. ! 'bnl' notation here follows Bengtsson, where n starts at 1 rather than zero-- i.e, b2l(Bengtsson) = b1l(Bmad) b2l2 = (b1_grad + b1) * value_of_attribute(ele2,'L') b3l2 = (b2_grad + b2) * value_of_attribute(ele2,'L')/2. etax2(ix) = etax2(ix) + (sqrt(ele%a%beta)/(2.*sin(ring%a%tune/2.)))*(b2l2 - b3l2*ele2%a%eta)*& ele2%a%eta * sqrt(ele2%a%beta) * cos(abs(ele%a%phi - ele2%a%phi)-ring%a%tune/2.) betax1(ix) = betax1(ix) + (ele%a%beta/(2.*sin(ring%a%tune)))*(b2l2 - 2.*b3l2*ele2%a%eta)*& (ele2%a%beta) * cos(2*abs(ele%a%phi - ele2%a%phi)-ring%a%tune) betay1(ix) = betay1(ix) + (-ele%b%beta/(2.*sin(ring%b%tune)))*(b2l2 - 2.*b3l2*ele2%a%eta)*& (ele2%b%beta) * cos(2*abs(ele%b%phi - ele2%b%phi)-ring%b%tune) endif enddo loop2 endif enddo loop1 end subroutine calc_dbeta_ddelta !============================================================================= !============================================================================= ! subroutine for calculating dnu_dJ: subroutine calc_dnu_dJ(ring, dnux_dJx, dnux_dJy, dnuy_dJy) type(lat_struct) :: ring type(ele_struct) :: ele, ele2 real(rp), intent(out) :: dnux_dJx, dnux_dJy, dnuy_dJy real(rp) :: b2l=0.,b3l=0., b2l2=0.,b3l2=0., b2_grad, b2 integer ix, jx dnux_dJx=0. dnux_dJy=0. dnuy_dJy=0. if (ring%ele(0)%ix_pointer .ne. -1) then ! have not filtered by k1/k2/b1/b2 yet. do it now. call filter_eles(ring) endif loop1: do ix=1, ring%n_ele_track if (ring%ele(ix)%ix_pointer .eq. 2) then ! filter down to elements with k2 or b2 only if (match_reg(ring%ele(ix)%name, "SK")) cycle loop1 ele = ring%ele(ix) call twiss_at_element (ring%ele(ix), average = ele) b2_grad = value_of_attribute(ele,'B2_GRADIENT') if (b2_grad .eq. real_garbage$) b2_grad = 0. b2 = value_of_attribute(ele,'B2') if (b2 .eq. real_garbage$) b2 = 0. !see this location in calc_h for comment b3l = (b2_grad + b2) * value_of_attribute(ele,'L')/2. loop2: do jx=1, ring%n_ele_track if (ring%ele(jx)%ix_pointer .eq. 2) then ! filter down to elements with k2 or b2 only if (match_reg(ring%ele(ix)%name, "SK")) cycle loop2 ele2 = ring%ele(jx) call twiss_at_element (ring%ele(jx), average = ele2) b2_grad = value_of_attribute(ele2,'B2_GRADIENT') if (b2_grad .eq. real_garbage$) b2_grad = 0. b2 = value_of_attribute(ele2,'B2') if (b2 .eq. real_garbage$) b2 = 0. b3l2 = (b2_grad + b2) * value_of_attribute(ele2,'L')/2. dnux_dJx = dnux_dJx + (-1/(16.*pi))*b3l*b3l2*(ele%a%beta)**(3/2.)*(ele2%a%beta)**(3/2.)*& ((3.*cos(abs(ele%a%phi-ele2%a%phi)-ring%a%tune/2.))/sin(ring%a%tune/2.) + & (cos(3.*abs(ele%a%phi-ele2%a%phi)-3.*ring%a%tune/2.))/sin(3.*ring%a%tune/2.)) dnux_dJy = dnux_dJy + (1/(8.*pi))*b3l*b3l2*sqrt(ele%a%beta*ele2%a%beta)*& ele%b%beta * ((2.*ele2%a%beta*cos(abs(ele%a%phi-ele2%a%phi)-ring%a%tune/2.))/& sin(ring%a%tune/2.) + (-ele2%b%beta*cos(abs((ele%a%phi-ele2%a%phi) + & 2.*(ele%b%phi-ele2%b%phi))-(ring%a%tune + 2.*ring%b%tune)/2.))/sin((ring%a%tune + & 2.*ring%b%tune)/2.) + (ele2%b%beta*cos(abs((ele%a%phi-ele2%a%phi) - 2.*(ele%b%phi-ele2%b%phi))& -(ring%a%tune - 2.*ring%b%tune)/2.))/sin((ring%a%tune - 2.*ring%b%tune)/2.)) dnuy_dJy = dnuy_dJy + (-1/(16.*pi))*b3l*b3l2*sqrt(ele%a%beta*ele2%a%beta)*& ele%b%beta * ele2%b%beta * ((4.*cos(abs(ele%a%phi-ele2%a%phi)-ring%a%tune/2.))/& sin(ring%a%tune/2.) + (cos(abs((ele2%a%phi-ele%a%phi) + 2.*(ele2%b%phi-ele%b%phi))-& (ring%a%tune + 2.*ring%b%tune)/2.))/sin((ring%a%tune + 2.*ring%b%tune)/2.) + & (cos(abs((ele2%a%phi-ele%a%phi) - 2.*(ele2%b%phi-ele%b%phi))-(ring%a%tune - & 2.*ring%b%tune)/2.))/sin((ring%a%tune - 2.*ring%b%tune)/2.)) endif enddo loop2 endif enddo loop1 end subroutine calc_dnu_dJ !============================================================================= !============================================================================= subroutine calc_xi2(ring, xix2, xiy2) type(lat_struct), intent(in) :: ring type(ele_struct) :: ele, ele2 real(rp), allocatable :: etax2(:), betax1(:), betay1(:) real(rp) :: b2l=0.,b3l=0., b1_grad, b1, b2_grad, b2 real(rp), intent(out) :: xix2,xiy2 integer ix, jx call calc_dbeta_ddelta(ring, etax2, betax1, betay1) ! initialize: xix2 = (-1/2.) * ring%a%chrom ! constant offset xiy2 = (-1/2.) * ring%b%chrom ! constant offset do ix=0, ring%n_ele_track if ((ring%ele(ix)%ix_pointer .eq. 1) .or. (ring%ele(ix)%ix_pointer .eq. 2)) then ! passes k1/k2/b1/b2 filter done in calc_dbeta_ddelta if (match_reg(ring%ele(ix)%name, "SK")) cycle ele = ring%ele(ix) call twiss_at_element (ring%ele(ix), average = ele) b1_grad = value_of_attribute(ele,'B1_GRADIENT') if (b1_grad .eq. real_garbage$) b1_grad = 0. b1 = value_of_attribute(ele,'B1') if (b1 .eq. real_garbage$) b1 = 0. b2_grad = value_of_attribute(ele,'B2_GRADIENT') if (b2_grad .eq. real_garbage$) b2_grad = 0. b2 = value_of_attribute(ele,'B2') if (b2 .eq. real_garbage$) b2 = 0. b2l = (b1_grad + b1) * value_of_attribute(ele,'L') b3l = (b2_grad + b2) * value_of_attribute(ele,'L')/2. xix2 = xix2 + (1/(8.*pi))*(2.*b3l*etax2(ix)*ele%a%beta - (b2l - 2.*b3l * ele%a%eta)*betax1(ix)) xiy2 = xiy2 + (-1/(8.*pi))*(2.*b3l*etax2(ix)*ele%b%beta + (b2l - 2.*b3l * ele%a%eta)*betay1(ix)) endif enddo end subroutine calc_xi2 !============================================================================= !============================================================================= subroutine meas_xi2(ring, delta, xix2, xiy2) type(lat_struct), intent(in) :: ring type(lat_struct) ring2(-1:1) type (coord_struct), allocatable :: orb(:) real(rp), intent(out) :: xix2,xiy2 real(rp), intent(in) :: delta real(rp) :: xix(-1:1)=0., xiy(-1:1)=0. integer ix integer :: i_dim = 4 allocate(orb(0:ring%n_ele_track)) ring2 = ring do ix = -1,1, 2 orb(0)%vec(6) = delta *ix call closed_orbit_calc(ring2(ix), orb, i_dim) call lat_make_mat6(ring2(ix), -1, orb) call twiss_at_start(ring2(ix)) call twiss_propagate_all(ring2(ix)) call chrom_calc(ring2(ix), delta, xix(ix), xiy(ix)) end do xix2 = 0.5*(xix(1) - xix(-1))/(2.*delta) xiy2 = 0.5*(xiy(1) - xiy(-1))/(2.*delta) end subroutine meas_xi2 !============================================================================= !============================================================================= subroutine meas_dbeta_ddelta(ring, delta, etax2, betax1, betay1) type(lat_struct), intent(in) :: ring type(lat_struct) ring2(-1:1) type (coord_struct), allocatable :: orb(:) real(rp), allocatable, intent(out) :: etax2(:), betax1(:), betay1(:) real(rp), intent(in) :: delta real(rp) :: xix(-1:1)=0., xiy(-1:1)=0. real(rp), allocatable :: etax(:,:), betax(:,:), betay(:,:) integer ix, jx integer :: i_dim = 4 allocate(orb(0:ring%n_ele_track),etax2(0:ring%n_ele_track),betax1(0:ring%n_ele_track),& betay1(0:ring%n_ele_track)) allocate(etax(-1:1,0:ring%n_ele_track), betax(-1:1,0:ring%n_ele_track), betay(-1:1,0:ring%n_ele_track)) ring2(-1) = ring ring2(1) = ring do ix = -1,1,2 orb(0)%vec(6) = delta *ix call closed_orbit_calc(ring2(ix), orb,i_dim) call lat_make_mat6(ring2(ix), -1, orb) call twiss_at_start(ring2(ix)) call twiss_propagate_all(ring2(ix)) betax(ix,:) = ring2(ix)%ele(0:ring2(ix)%n_ele_track)%a%beta betay(ix,:) = ring2(ix)%ele(0:ring2(ix)%n_ele_track)%b%beta etax(ix,:) = ring2(ix)%ele(0:ring2(ix)%n_ele_track)%a%eta enddo betax1(:) = (betax(1,:) - betax(-1,:))/(2.*delta) betay1(:) = (betay(1,:) - betay(-1,:))/(2.*delta) etax2(:) = 0.5*( etax(1,:) - etax(-1,:))/(2.*delta) ! etax2 == (d_etax/d_delta)/2 end subroutine meas_dbeta_ddelta !============================================================================= !============================================================================= subroutine filter_eles(ring) type(lat_struct), intent(inout) :: ring integer ix do ix=1, ring%n_ele_track if (value_of_attribute(ring%ele(ix), 'L') .eq. 0.) cycle ! flag k2/b2 first, since that takes precidence when both k1/b1 and k2/b2 are present. if((value_of_attribute(ring%ele(ix), 'K2') .ne. 0 .and. value_of_attribute(ring%ele(ix), 'K2') .ne. real_garbage$) & .or. (value_of_attribute(ring%ele(ix), 'B2') .ne. 0 .and. value_of_attribute(ring%ele(ix), 'B2') .ne. real_garbage$)) then ring%ele(ix)%ix_pointer = 2 elseif ( (value_of_attribute(ring%ele(ix), 'K1') .ne. 0 .and. value_of_attribute(ring%ele(ix), 'K1') .ne. real_garbage$) & .or. (value_of_attribute(ring%ele(ix), 'B1') .ne. 0 .and. value_of_attribute(ring%ele(ix), 'B1') .ne. real_garbage$)) then ring%ele(ix)%ix_pointer = 1 endif enddo ring%ele(0)%ix_pointer = -1 ! set this to flag that we've already done the filtering end subroutine filter_eles end module calc_res_mod