!!! Header ! Add comments about module here module tune_shift_mod use bmad type wake_values_struct real(rp) :: z=0., w=0., w_in=0. end type wake_values_struct type tuneshift_struct real(rp) :: x=0., y=0., z_pos=0., z_disp=0. end type tuneshift_struct type wake_ele_struct character(60) :: name = '', material = '', shape = '' character(200) :: filepath = '' logical :: use = .true. real(rp) :: sig_z_calc = 0., x_disp = 0., y_disp = 0.00500, a = 0., b = 0., dz = 0., rho = 0., d_y = 0. , wake_scale = (5.0/4.6)**(3.0/2.0) type(wake_values_struct) :: wake_values(10000) integer :: n_step = 0 end type wake_ele_struct type ele_tuneshift_struct character(60) :: name = '' real(rp) :: s0 = 0., s1 = 0., beta_x_avg = 0., beta_y_avg = 0. type(tuneshift_struct) :: tuneshift(1000) logical :: use = .true. end type ele_tuneshift_struct type d_func_struct real(rp) :: q=0., f0=0., Fx=0., Fy=0., G0=0., G1x=0., G1y=0. end type d_func_struct type Q_values_struct real(rp) :: z_hat=0., del_Q_min=2., del_Q_max=0., del_Q_avg=0., del_Q_rms=0. end type Q_values_struct contains !===================================================================== !--- subroutine parse_wakes ! Reads in wake and resistive wall definitions from wake_list namelist subroutine parse_wakes(init_file, wakes, sig_z) use fgsl implicit none character(200) :: init_file, line='' type(wake_ele_struct) :: wakes(100) integer :: ix = 0, iz = 0, j = 1, i = 1, iost = 0, ios = 0 real(rp) :: sig_z, d_y = 0., c, Z0, rho, i_eqn, i_eqn_neg, k_eqn namelist /wake_list/ wakes open(1, file=init_file) read(1, nml=wake_list) close(1) write(*,*) "sig_z = ", sig_z do ix=1,size(wakes) ! if no data, or if intended to skip, skip it if (trim(wakes(ix)%name) == '') cycle if (wakes(ix)%use .eqv. .false.) cycle ! Two scenarios: ! 1) if filepath is defined, it's a T3P output, therefore discontinuity. write(*,*) '***', ix, wakes(ix)%name if (wakes(ix)%filepath /= '') then ! read in z and w values from file open(2, file=wakes(ix)%filepath) ! formatted read of individual lines in wakes(ix)%file, populating wakes(ix)%wake_values%z,w_in j=0 ios=0 do while (ios==0) read(2,'(a)',iostat=ios)line if (line(:1)=='#') cycle j=j+1 read(line,'(F24.18,1x,F25.18)') wakes(ix)%wake_values(j)%z, wakes(ix)%wake_values(j)%w_in if (ix == 1) then write(90 ,*) wakes(ix)%wake_values(j)%z, wakes(ix)%wake_values(j)%w_in end if ! offset all z entries by -5sig_z wakes(ix)%wake_values(j)%z = wakes(ix)%wake_values(j)%z - (5*wakes(ix)%sig_z_calc) enddo ! set dz and n_step wakes(ix)%dz = wakes(ix)%wake_values(2)%z - wakes(ix)%wake_values(1)%z wakes(ix)%n_step = j close(2) ! 2) Resistive wall, computed analytically: else wakes(ix)%n_step = int(10*sig_z/wakes(ix)%dz) ! translate material name to resistivity: rho = material_to_rho(wakes(ix)%material) Z0 = 1/(c_light*eps_0_vac) ! set d_y (geometric prefactor) d_y = evaluate_d_y(wakes(ix)%shape, wakes(ix)%a, wakes(ix)%b) ! generate wake curve, evaluated over -5sigma to +5sigma in steps of dz do iz=1, wakes(ix)%n_step wakes(ix)%wake_values(iz)%z = -5*sig_z + wakes(ix)%dz*iz ! evaluate bessel functions !write(23,'(i5,3es14.4)') ix, wakes(ix)%wake_values(iz)%z, wakes(ix)%wake_values(iz)%z**2/(4*sig_z**2), fgsl_sf_bessel_knu(dble(0.25), wakes(ix)%wake_values(iz)%z**2/(4*sig_z**2)) if (wakes(ix)%wake_values(iz)%z == 0) then ! knu blows up at z=0 ! so set k_eqn = (k_eqn[iz-1]+k_eqn[iz+1])/2 !k_eqn = (fgsl_sf_bessel_knu(dble(0.25),wakes(ix)%wake_values(iz-1)%z**2/(4*sig_z**2)) + & ! fgsl_sf_bessel_knu(dble(0.25),wakes(ix)%wake_values(iz+1)%z**2/(4*sig_z**2)))/2 !i_eqn = (fgsl_sf_bessel_Inu(dble(0.25), wakes(ix)%wake_values(iz-1)%z**2/(4*sig_z**2)) + & ! fgsl_sf_bessel_Inu(dble(0.25), wakes(ix)%wake_values(iz+1)%z**2/(4*sig_z**2)))/2 !i_eqn_neg = (2*SIN(pi/4)/pi)*k_eqn + i_eqn k_eqn = fgsl_sf_bessel_knu(dble(0.25),wakes(ix)%wake_values(iz-1)%z**2/(4*sig_z**2)) i_eqn = fgsl_sf_bessel_Inu(dble(0.25), wakes(ix)%wake_values(iz-1)%z**2/(4*sig_z**2)) i_eqn_neg = (2*SIN(pi/4)/pi)*k_eqn + i_eqn else k_eqn = fgsl_sf_bessel_knu(dble(0.25), wakes(ix)%wake_values(iz)%z**2/(4*sig_z**2)) i_eqn = fgsl_sf_bessel_Inu(dble(0.25), wakes(ix)%wake_values(iz)%z**2/(4*sig_z**2)) i_eqn_neg = (2*SIN(pi/4)/pi)*k_eqn + i_eqn endif ! compute wake; requires separate calculation for z < 0 and z > 0 (Changed w's to w_in) if (wakes(ix)%wake_values(iz)%z > 0) then wakes(ix)%wake_values(iz)%w_in = d_y*(c_light/(pi**2*wakes(ix)%b**3))*SQRT(2*Z0*rho/sig_z)*(pi/4)*& exp(-(wakes(ix)%wake_values(iz)%z**2)/(4*sig_z**2))*SQRT(wakes(ix)%wake_values(iz)%z/sig_z)*(i_eqn_neg + i_eqn)*1E-12 else if (wakes(ix)%wake_values(iz)%z < 0) then wakes(ix)%wake_values(iz)%w_in = d_y*(c_light/(pi**2*wakes(ix)%b**3))*SQRT(2*Z0*rho/sig_z)*(1/(8**(0.5)))*& exp(-((wakes(ix)%wake_values(iz)%z**2)/(4*sig_z**2)))*SQRT(-wakes(ix)%wake_values(iz)%z/sig_z)*k_eqn*1E-12 end if write(88,*) wakes(ix)%wake_values(iz)%z, wakes(ix)%wake_values(iz)%w_in, wakes(ix)%name, wakes(ix)%material, rho write(89,*) wakes(ix)%wake_values(iz)%z, wakes(ix)%wake_values(iz)%w_in !if (ix == 2) then !write(12,*) wakes(ix)%wake_values(iz)%z !end if enddo endif enddo end subroutine parse_wakes !===================================================================== !--- subroutine scale_wakes ! reconvolves T3P wakes to a new bunch length subroutine scale_wakes(wakes, sig_z) implicit none type(wake_ele_struct) :: wakes(100) real(rp) :: sig_z, sig_c, const, wjx = 0., z0 = 0., z1 = 0., z2 = 0., dz = 0. integer :: i = 0, ix = 0, jx = 0, n_z_step = 0 print *, "sig_z = ", sig_z print *, "wakes(1)%sig_z_calc = ", wakes(1)%sig_z_calc do i=1, size(wakes) if (trim(wakes(i)%name) == '') cycle if (wakes(i)%use .eqv. .false.) cycle ! check that we are reconvolving to a longer bunch length if (sig_z <= wakes(i)%sig_z_calc) then print *, "*** ERROR (Rtn scale_wakes): requested reconvolved bunch length not larger than calculated" print *, " iwake = ", i print *, " sig_z = ", sig_z print *, " sig_z_calc = ", wakes(i)%sig_z_calc exit end if ! set reconvolve Gaussian sigma sig_c = sqrt(sig_z**2 - wakes(i)%sig_z_calc**2) ! Set const const = 1/SQRT(2*pi*sig_c**2) ! set constants z0, dz, and n_z_step z0 = wakes(i)%wake_values(1)%z dz = wakes(i)%wake_values(2)%z - wakes(i)%wake_values(1)%z n_z_step = 5*sig_c/dz ! +/- 5 sigma ! Scale wakes if (wakes(i)%filepath /= '') then wakes(i)%wake_values(:)%w = 0. do ix=1, wakes(i)%n_step do jx=-n_z_step, n_z_step ! Gaussian is centered at z1 z1 = wakes(i)%wake_values(ix)%z z2 = z1 + jx*dz ! value of wake at this jx if (ix+jx < 1) then ! 0 if before head of bunch wjx = 0.0 else if (ix+jx > wakes(i)%n_step) then ! last data point value if past valid data in tail wjx = wakes(i)%wake_values(wakes(i)%n_step)%w_in else wjx = wakes(i)%wake_values(ix+jx)%w_in end if ! convolution wakes(i)%wake_values(ix)%w = wakes(i)%wake_values(ix)%w & + wjx * exp(-(z2-z1)**2/(2*sig_c**2))/(sqrt(2*pi)*sig_c) * dz enddo enddo end if enddo end subroutine scale_wakes !===================================================================== !--- subroutine calc_tuneshift ! computes the contribution to tuneshift from each element subroutine calc_tuneshift(init_file, wakes, ele_tuneshifts, ring, sig_z, current) implicit none character(200) :: init_file real(rp) :: betay_avg = 0., energy = 0., wake = 0., w_val = 0., z_val = 2., z_val2 = 2., sig_z, value1 = 0., wake_val = 0., del_Q_theta = 0., del_theta = 0., del_Q_obs = 0., delta = 0., Q_tot = 0., Q_rms_tot = 0., current real(rp) :: wake_scale = (5.0/4.5)**(3.0/2.0) type(wake_ele_struct) :: wakes(100) type(ele_tuneshift_struct) :: ele_tuneshifts(200) type(Q_values_struct) :: Q_values(100) type(lat_struct) :: ring integer :: ix=0, jx=0, kx=0, j=0, z0 = 1, z1 = 1, k=1, n=1 namelist /tuneshift_list/ ele_tuneshifts open(1, file=init_file) read(1, nml=tuneshift_list) close(1) ! save energy for later use: energy = value_of_attribute(ring%ele(0), 'E_TOT') open(1,file='tuneshift_output_resistivewall.txt', status='replace') open(2,file='tuneshift_output_wake.txt', status='replace') ! iterate over elements from ele_tuneshifts array: do kx=1, size(ele_tuneshifts) if (trim(ele_tuneshifts(kx)%name) == '') cycle !if ((trim(ele_tuneshifts(kx)%name) /= 'ccu wake')) cycle !if ((trim(ele_tuneshifts(kx)%name) /= 'ccu resistive')) cycle do ix=1, size(wakes) ! only proceed if name of this template (in "wakes") matches requested name in ele_tuneshifts if (wakes(ix)%name == ele_tuneshifts(kx)%name) then !write(59,*) wakes(ix)%name ! reset z_val = 2 z_val2 = 2 ele_tuneshifts(kx)%use = wakes(ix)%use if (wakes(ix)%use .eqv. .false.) cycle ! Calculate average beta function betay_avg = average_beta(ring, ele_tuneshifts(kx)%s0, ele_tuneshifts(kx)%s1, wakes(ix)%filepath) ! Find values of z closest to 0 and set z0,z1 indices call find_z(wakes(ix)%wake_values(:)%z, wakes(ix)%n_step, 0._rp, z0, z1) do j=1, wakes(ix)%n_step w_val = wakes(ix)%wake_values(j)%w_in ! Again there are two cases: ! 1) if filepath is specified, it's a T3P wake from discontinuity if (wakes(ix)%filepath /= '') then ele_tuneshifts(kx)%tuneshift(j)%y = current*(-1)*(1/energy)*2.56E3*w_val*betay_avg/(4*pi) if(wakes(ix)%name == 'ccu wake') then ! hack - t3p run for 5mm, need 4.6mm aperture ele_tuneshifts(kx)%tuneshift(j)%y = wake_scale*ele_tuneshifts(kx)%tuneshift(j)%y end if ele_tuneshifts(kx)%tuneshift(j)%z_pos = wakes(ix)%wake_values(j)%z write(99,*) ele_tuneshifts(kx)%tuneshift(j)%y, j write(45,*) ele_tuneshifts(kx)%tuneshift(j)%y*390.10*1000 ! 2) if filepath is NOT specified, it's resistive wall; need length of element too else ele_tuneshifts(kx)%tuneshift(j)%y = (-1)*(1/energy)*current*(2.56E3*w_val)*betay_avg*(ele_tuneshifts(kx)%s1 - ele_tuneshifts(kx)%s0)/(4*pi) ele_tuneshifts(kx)%tuneshift(j)%z_pos = wakes(ix)%wake_values(j)%z write(44,*) ele_tuneshifts(kx)%tuneshift(j)%y*390.10*1000, w_val endif enddo ! Find tune shift at z = 0 wake_val = -((ele_tuneshifts(kx)%tuneshift(z1)%y - ele_tuneshifts(kx)%tuneshift(z0)%y)/(wakes(ix)%wake_values(z1)%z - wakes(ix)%wake_values(z0)%z))*wakes(ix)%wake_values(z0)%z + ele_tuneshifts(kx)%tuneshift(z0)%y ! Print out the value of the tuneshift at z = 0 for each element write(40,*)'wake_val', ele_tuneshifts(kx)%name, wake_val, wakes(ix)%filepath value1 = value1 + wake_val ! Print out the summed total of the tuneshifts at z=0 for all elements so far write(41,*)'value', ele_tuneshifts(kx)%name, value1, wakes(ix)%filepath end if enddo enddo close(1) close(2) !--------- ! calculate tuneshift vs. synchrotron osciallation amplitude, in 26 steps of z-amplitude do k=0, 25 ! Set maximum amplitude of position oscilations as z_hat Q_values(k+1)%z_hat = 5*sig_z*k/25 Q_values(k+1)%del_Q_max = -999999 Q_values(k+1)%del_Q_rms = 0 !----- ! evaluate in 10-degree increments around one synchrotron oscillation do n=0, 36 do kx=1, size(ele_tuneshifts) ! Do no calculations for elements not included if (ele_tuneshifts(kx)%use .eqv. .false.) cycle ! Do no calculations for elements without any information if (trim(ele_tuneshifts(kx)%name) == '') cycle !if ((trim(ele_tuneshifts(kx)%name) /= 'ccu wake')) cycle !if ((trim(ele_tuneshifts(kx)%name) /= 'ccu resistive')) cycle if (Q_values(k+1)%z_hat == 0) write(47,*) ele_tuneshifts(kx)%tuneshift(Q_values(k+1)%z_hat)%y*390.10*1000 ! reset z_val = 2 z_val2 = 2 ! Find delta Theta for each value n from 0 to 36 (making 10 degeree steps) del_theta = 2*pi*n/36 ! Find the two z positions closest to z_hat do jx=1, size(ele_tuneshifts(kx)%tuneshift) if (ele_tuneshifts(kx)%tuneshift(jx)%z_pos == 0) cycle if (ABS(ele_tuneshifts(kx)%tuneshift(jx)%z_pos - Q_values(k+1)%z_hat*sin(del_theta)) < z_val) then z1 = z0 z0 = jx z_val = ABS(ele_tuneshifts(kx)%tuneshift(z0)%z_pos - Q_values(k+1)%z_hat*sin(del_theta)) z_val2 = ABS(ele_tuneshifts(kx)%tuneshift(z1)%z_pos - Q_values(k+1)%z_hat*sin(del_theta)) end if if (ABS(ele_tuneshifts(kx)%tuneshift(jx)%z_pos - Q_values(k+1)%z_hat*sin(del_theta)) < z_val2) then if (jx /= z0) then z1 = jx z_val2 = ABS(ele_tuneshifts(kx)%tuneshift(z1)%z_pos - Q_values(k+1)%z_hat*sin(del_theta)) end if end if if (z1 == z0) z1 = z1+1 end do ! over jx ! Interpolate using the two closest z positons to find tuneshift at each z_hat del_Q_theta = del_Q_theta + ((ele_tuneshifts(kx)%tuneshift(z1)%y - ele_tuneshifts(kx)%tuneshift(z0)%y)/(ele_tuneshifts(kx)%tuneshift(z1)%z_pos - ele_tuneshifts(kx)%tuneshift(z0)%z_pos))*(Q_values(k+1)%z_hat*sin(del_theta) - ele_tuneshifts(kx)%tuneshift(z0)%z_pos) + ele_tuneshifts(kx)%tuneshift(z0)%y write(43,*) del_Q_theta*390.10*1000, n end do ! over kx ! Set max or min if del_Q_theta is either if (del_Q_theta > Q_values(k+1)%del_Q_max) Q_values(k+1)%del_Q_max = del_Q_theta if (del_Q_theta < Q_values(k+1)%del_Q_min) Q_values(k+1)%del_Q_min = del_Q_theta ! Sum together all values of del_Q_theta to start solving for del_Q_avg if (n /= 36) Q_values(k+1)%del_Q_avg = Q_values(k+1)%del_Q_avg + del_Q_theta del_Q_theta = 0 end do ! over n ! Divide summed del_Q_thetas by 36 to find del_Q_avg Q_values(k+1)%del_Q_avg = Q_values(k+1)%del_Q_avg / 36 !--------- ! Repeat steps from above to cycle through del_Q_theta to calculate del_Q_rms do n=0, 36 ! Find delta Theta for each value n from 0 to 36 (making 10 degeree steps) del_theta = 2*pi*n/36 do kx=1, size(ele_tuneshifts) if (ele_tuneshifts(kx)%use .eqv. .false.) cycle if (trim(ele_tuneshifts(kx)%name) == '') cycle !if ((trim(ele_tuneshifts(kx)%name) /= 'ccu wake')) cycle !if ((trim(ele_tuneshifts(kx)%name) /= 'ccu resistive')) cycle ! reset z_val = 2 z_val2 = 2 ! Find the two z positions closest to z_hat do jx=1, wakes(1)%n_step if (ele_tuneshifts(kx)%tuneshift(jx)%z_pos == 0) cycle if (ABS(ele_tuneshifts(kx)%tuneshift(jx)%z_pos - Q_values(k+1)%z_hat*sin(del_theta)) < z_val) then z1 = z0 z0 = jx z_val = ABS(ele_tuneshifts(kx)%tuneshift(z0)%z_pos - Q_values(k+1)%z_hat*sin(del_theta)) z_val2 = ABS(ele_tuneshifts(kx)%tuneshift(z1)%z_pos - Q_values(k+1)%z_hat*sin(del_theta)) end if if (ABS(ele_tuneshifts(kx)%tuneshift(jx)%z_pos - Q_values(k+1)%z_hat*sin(del_theta)) < z_val2) then if (jx /= z0) then z1 = jx z_val2 = ABS(ele_tuneshifts(kx)%tuneshift(z1)%z_pos - Q_values(k+1)%z_hat*sin(del_theta)) end if end if if (z1 == z0) z1 = z1+1 end do ! over jx ! Interpolate using the two closest z positons to find tuneshift at each z_hat del_Q_theta = del_Q_theta + ((ele_tuneshifts(kx)%tuneshift(z1)%y - ele_tuneshifts(kx)%tuneshift(z0)%y)/(ele_tuneshifts(kx)%tuneshift(z1)%z_pos - & ele_tuneshifts(kx)%tuneshift(z0)%z_pos))*Q_values(k+1)%z_hat*sin(del_theta) - ((ele_tuneshifts(kx)%tuneshift(z1)%y -& ele_tuneshifts(kx)%tuneshift(z0)%y)/(ele_tuneshifts(kx)%tuneshift(z1)%z_pos - ele_tuneshifts(kx)%tuneshift(z0)%z_pos))*ele_tuneshifts(kx)%tuneshift(z0)%z_pos + & ele_tuneshifts(kx)%tuneshift(z0)%y end do ! over kx ! Sum together all values of (del_Q_theta - del_Q_avg) to start solving for del_Q_rms if (n /= 36) Q_values(k+1)%del_Q_rms = Q_values(k+1)%del_Q_rms + (del_Q_theta - Q_values(k+1)%del_Q_avg)**2 write(93,*) Q_values(k+1)%del_Q_rms, k del_Q_theta = 0 end do ! over n ! Divide summed del_Q_thetas by 36 and squareroot to find del_Q_rms Q_values(k+1)%del_Q_rms = SQRT(Q_values(k+1)%del_Q_rms / 36) if (k /= 25) Q_tot = Q_tot + k*(((5*sig_z)/25.0)**2)*EXP((((k*5*sig_z)/25.0)**2)/(-2*(sig_z**2)))*(Q_values(k+1)%del_Q_avg) if (k /= 25) Q_rms_tot = Q_rms_tot + k*(((5*sig_z)/25)**2)*EXP((((k*5*sig_z)/25)**2)/(-2*(sig_z)**2))*(Q_values(k+1)%del_Q_rms**2) end do ! over k Q_tot = 390.10*1000*(Q_tot)/(sig_z**2) Q_rms_tot = 390.10*1000*(SQRT(Q_rms_tot))/sig_z write(65,*) Q_tot, Q_rms_tot write(*,*) '' write(*,'(a,f8.3,a)') ' I_bunch = ', current, ' mA' write(*,'(a,f8.3,a)') 'Vertical tune shift = ', Q_tot, ' Hz' write(*,'(a,f8.3,a)') 'Vertical tune spread = ', Q_rms_tot, ' Hz' ! Write del_Q values to a file open(1,file='Q_values.txt', status='old') write(1,'(a20, 6a23)')'k+1', 'z_hat', 'del_Q_min', 'del_Q_max', 'del_Q_avg', 'del_Q_rms', 'del_Q_obs' do k=0,25 del_Q_obs = (5*k/25.0)*exp(((5*k/25.0)**2)/(-2.0))*Q_values(k+1)%del_Q_avg write(1,'(I20, 6F23.15)')k+1, Q_values(k+1)%z_hat, & Q_values(k+1)%del_Q_min, Q_values(k+1)%del_Q_max, Q_values(k+1)%del_Q_avg, Q_values(k+1)%del_Q_rms, del_Q_obs delta = del_Q_obs + delta !write(11,*) Q_values(k+1)%z_hat !write(89,*) Q_values(k+1)%del_Q_min !write(13,*) Q_values(k+1)%del_Q_max !write(14,*) Q_values(k+1)%del_Q_avg !write(15,*) Q_values(k+1)%del_Q_rms end do close(1) delta = delta/26 end subroutine calc_tuneshift !===================================================================== !--- subroutine find_z ! locates z indices closest to z_target subroutine find_z(z, n_step, z_target, z0, z1) implicit none real(rp), intent(in) :: z(:), z_target real(rp) z_val, z_val2 integer, intent(out) :: z0, z1 integer jx integer, intent(in) :: n_step z0 = 0 z1 = 0 z_val = 1000. z_val2 = 1000. do jx=1, n_step if (ABS(z(jx)-z_target) < z_val) then z1 = z0 z0 = jx z_val = ABS(z(z0)-z_target) z_val2 = ABS(z(z1)-z_target) endif if (ABS(z(jx)-z_target) < z_val2) then if (jx /= z0) z1 = jx end if enddo end subroutine find_z !===================================================================== !===================================================================== !=== FUNCTIONS ======================================================= !===================================================================== !--- function material_to_rho ! translates a material name (str) into a resistivity function material_to_rho(material) implicit none character(20) :: material real(rp) :: material_to_rho, stainless_steel = 6.897E-7, aluminum = 2.65E-8 Select Case (material) Case ('stainless steel') material_to_rho = stainless_steel Case ('aluminum') material_to_rho = aluminum End Select end function material_to_rho !===================================================================== !--- function evaluate_d_y ! evaluates the geometric prefactor D_y for resistive wall function evaluate_d_y(shape, a, b) implicit none character(20) :: shape character(200) :: line='' type(d_func_struct) :: funcs(150) real(rp) :: evaluate_d_y, a, b, q=0., x=2., y = 2. integer :: j=0, ix=0, ij=1, ios=0 namelist /func_list/ funcs ! open table containing F and G values open(1, file='/home/shanksj/bmad_dist/bsim_cesr/tuneshift_from_impedance/f_and_g.txt', iostat=ios) j=0 do while (ios==0) ! enter tabel into array read(1,'(a)',iostat=ios)line line=adjustl(line) if (line(:1)=='!') cycle j=j+1 read(line,'(F6.4,10x,F6.4,10x,F6.4,10x,F6.4,10x,F6.4,10x,F6.4,10x,F6.4)')funcs(j)%q, funcs(j)%f0, & funcs(j)%Fx, funcs(j)%Fy, funcs(j)%G0, funcs(j)%G1x, funcs(j)%G1y end do close(1) ! find q q=(1-b/a)/(1+b/a) ! find closest two points in array to q do j=1,size(funcs) if (ABS(funcs(j)%q - q) < x) then ix=ij ij=j x = ABS(funcs(ij)%q - q) y = ABS(funcs(ix)%q - q) end if if (ABS(funcs(j)%q - q) < y) then if (j /= ij) ix = j end if end do ! evaluate F or G depending on shape of beam pipe if (shape=='rectangular')then evaluate_d_y = (funcs(ij)%Fy-funcs(ix)%Fy)/(funcs(ij)%q-funcs(ix)%q)*q-(funcs(ij)%Fy-& funcs(ix)%Fy)/(funcs(ij)%q-funcs(ix)%q)*funcs(ix)%q+funcs(ix)%Fy else evaluate_d_y = (funcs(ij)%G1y-funcs(ix)%G1y)/(funcs(ij)%q-funcs(ix)%q)*q-(funcs(ij)%G1y-& funcs(ix)%G1y)/(funcs(ij)%q-funcs(ix)%q)*funcs(ix)%q+funcs(ix)%G1y end if end function evaluate_d_y !===================================================================== !--- function average_beta ! averages beta functions over the length of an element function average_beta(ring, s0, s1, filepath) implicit none type(lat_struct) :: ring type(ele_struct) :: twiss0, twiss1, twissb character(200) :: filepath real(rp) :: average_beta, s0, s1, sb=0., sum_beta=0., delta_l = 0. integer :: n=0, i=0 ! Determine which method to use to calculate the average beta sum_beta = 0. if (filepath == '') then if (s1-s0 > 0.1) then ! Evaluate beta at every 10cm n = nint((s1-s0)/0.1) delta_l = (s1-s0)/n call twiss_and_track_at_s(ring, s0, twiss0) sum_beta = sum_beta + twiss0%b%beta/2 call twiss_and_track_at_s(ring, s1, twiss1) sum_beta = sum_beta + twiss1%b%beta/2 do i=1,n-1 sb=s0+i*delta_l ! write(*,*)'sb =', sb call twiss_and_track_at_s(ring, sb, twissb) sum_beta = sum_beta + twissb%b%beta ! write(*,*)'average_beta =', average_beta end do average_beta = sum_beta/n else ! Calculate average beta call twiss_and_track_at_s(ring, s0, twiss0) call twiss_and_track_at_s(ring, s1, twiss1) ! Return average average_beta = (twiss0%b%beta + twiss1%b%beta)/2 end if else ! Calculate average beta call twiss_and_track_at_s(ring, s0, twiss0) call twiss_and_track_at_s(ring, s1, twiss1) ! Return average average_beta = (twiss0%b%beta + twiss1%b%beta)/2 end if end function average_beta end module tune_shift_mod