module ring_aperture_mod use bmad_struct !use quick_plot !use random_mod !$ use omp_lib !use ptc_layout_mod use beam_mod use quick_plot implicit none type amplitude_info_struct type (coord_struct) :: orb real(rp) :: Qa0 = 0.0_rp real(rp) :: Qb0 = 0.0_rp real(rp) :: Qa1 = 0.0_rp real(rp) :: Qb1 = 0.0_rp integer :: n_turns end type type amplitude_block_struct real(rp) :: a_min = 0.0_rp real(rp) :: b_min = 0.0_rp real(rp) :: pz_min = 0.0_rp real(rp) :: a_max = 0.001_rp real(rp) :: b_max = 0.001_rp real(rp) :: pz_max = 0.001_rp integer :: Na = 10 integer :: Nb = 1 integer :: Npz = 10 type (amplitude_info_struct), allocatable :: info(:,:,:) end type contains !-------------------------------------------------------------------- !- ! ! !- !-------------------------------------------------------------------- subroutine write_amplitude_block(iu, amplitude_block) implicit none type (amplitude_block_struct) :: amplitude_block integer :: iu, ia, ib, ipz do ia=1, size(amplitude_block%info, 1) do ib=1, size(amplitude_block%info, 2) do ipz=1, size(amplitude_block%info, 3) !write (iu, '(6es15.7)') amplitude_block%info(ia, ib, ipz)%orb%vec call write_amplitude_info(iu, amplitude_block%info(ia, ib, ipz)) end do end do end do end subroutine subroutine write_amplitude_info(iu, amplitude_info) implicit none type (amplitude_info_struct) amplitude_info integer :: iu write (iu, '(10es15.7)') amplitude_info%orb%vec, & amplitude_info%Qa0, & amplitude_info%Qa1, & amplitude_info%Qb0, & amplitude_info%Qb1 end subroutine !-------------------------------------------------------------------- !- ! ! !- !-------------------------------------------------------------------- function tune_diffusion(Qa0, Qb0, Qa1, Qb1, n_turns) result (d) implicit none real(rp) :: d, Qa0, Qb0, Qa1, Qb1 integer :: n_turns d = log( hypot(Qa0-Qa1, Qb0-Qb1)/n_turns) end function !-------------------------------------------------------------------- !- ! ! !- !-------------------------------------------------------------------- subroutine frequency_drift(lat, amplitude_block, n_turns, plot_on) implicit none type (lat_struct) :: lat type (coord_struct) :: orb0 type (amplitude_block_struct) :: amplitude_block real(rp) :: Qa0, Qb0, Qa1, Qb1 real(rp) :: diffusion_rate integer :: iu, ia, ib, ipz, n_turns logical :: plot_on, err ! do ib=1, size(amplitude_block%info, 2) do ipz=1, size(amplitude_block%info, 3) do ia=1, size(amplitude_block%info, 1) !write (*, '(6es15.7)') amplitude_block%info(ia, ib, ipz)%orb%vec orb0=amplitude_block%info(ia, ib, ipz)%orb call tunes_from_tracking(lat, orb0, n_turns/2, Qa0, Qb0) call tunes_from_tracking(lat, orb0, n_turns/2, Qa1, Qb1, err) diffusion_rate = tune_diffusion(Qa0, Qb0, Qa1, Qb1, n_turns/2) write (*, '(6es15.7, 5f13.7)') amplitude_block%info(ia, ib, ipz)%orb%vec, & Qa0, Qb0, abs(Qa0-Qa1), abs(Qb0-Qb1), diffusion_rate amplitude_block%info(ia, ib, ipz)%Qa0 = Qa0 amplitude_block%info(ia, ib, ipz)%Qb0 = Qb0 amplitude_block%info(ia, ib, ipz)%Qa1 = Qa1 amplitude_block%info(ia, ib, ipz)%Qb1 = Qb1 if (plot_on .and. (.not. err)) then call qp_set_box(1, 1, 2, 1) call qp_calc_and_set_axis ('X', 0.0_rp, 0.01_rp, 4, 8, 'GENERAL') call qp_calc_and_set_axis ('Y', 0.0_rp, 0.01_rp, 4, 8, 'GENERAL') call qp_draw_symbol( & amplitude_block%info(ia, ib, ipz)%orb%vec(1), & amplitude_block%info(ia, ib, ipz)%orb%vec(6), & color=qp_continuous_color(-diffusion_rate/20)) call qp_set_box(2, 1, 2, 1) call qp_calc_and_set_axis ('X', 0.0_rp, 0.5_rp, 4, 8, 'GENERAL') call qp_calc_and_set_axis ('Y', 0.0_rp, 0.5_rp, 4, 8, 'GENERAL') call qp_draw_symbol(Qa0, Qb0, color=green$) call qp_draw_symbol(Qa1, Qb1, color=qp_continuous_color(-diffusion_rate/20)) endif end do end do end do end subroutine !-------------------------------------------------------------------- !- ! ! !- !-------------------------------------------------------------------- subroutine amplitude_block_init(lat, amplitude_block) use bmad implicit none type(lat_struct) :: lat type(amplitude_block_struct) :: amplitude_block type(ele_struct), pointer :: ele0 type(coord_struct) :: orb_mode real(rp) :: a0, b0, pz0, da, db, dpz integer :: Na, Nb, Npz integer :: i, ia, ib, ipz logical :: err character(30), parameter :: r_name='amplitude_block_init' ! ele0 => lat%ele(0) Na = amplitude_block%Na Nb = amplitude_block%Nb Npz = amplitude_block%Npz allocate(amplitude_block%info(Na, Nb, Npz)) ! Step sizes da = (amplitude_block%a_max-amplitude_block%a_min)/max(Na-1, 1) db = (amplitude_block%b_max-amplitude_block%b_min)/max(Nb-1, 1) dpz = (amplitude_block%pz_max-amplitude_block%pz_min)/max(Npz-1, 1) call init_coord(orb_mode, ele=ele0, element_end=upstream_end$) do ia = 1, Na do ib = 1, Nb do ipz = 1, Npz a0 = amplitude_block%a_min + (ia-1)*da b0 = amplitude_block%b_min + (ib-1)*db pz0 = amplitude_block%pz_min + (ipz-1)*dpz orb_mode%vec=[a0, 0.0_rp, b0 , 0.0_rp, 0.0_rp, pz0] !write (*, '(6es15.7)') orb_mode%vec call convert_coords('MODE', orb_mode, ele0, 'LAB', & amplitude_block%info(ia, ib, ipz)%orb, err) if (err) print *, 'err in '//r_name end do end do end do end subroutine !-------------------------------------------------------------------- !- ! ! !- !-------------------------------------------------------------------- subroutine amplitude_scan(lat, orb0, normal_mode, plot_on, max_turns, & n_sigma_a_max, n_sigma_b_max, n_sigma_z_max, & step_a, step_b, step_z) implicit none type(lat_struct) :: lat type (normal_modes_struct) :: normal_mode type(coord_struct), allocatable :: orbit(:) type(coord_struct) ::orb0, orb_J_phi type(ele_struct), pointer :: ele0 !type (aperture_struct) aperture real(rp) :: sigma_x, sigma_y, sigma_z, sigma_px, sigma_py, sigma_pz,emit_a, emit_b real(rp) :: Ja, Jb, Jz, dJa, dJb, dJz, max_Ja, max_Jb, max_Jz real(rp) :: n_sigma_a_max, n_sigma_b_max, n_sigma_z_max, step_a, step_b, step_z real(rp), allocatable :: x_arr(:), y_arr(:) real(rp) :: vec0(6) integer :: n, na, nb, nz, nJa, nJb, nJz, max_turns, n_lost integer :: track_state, ix_branch integer :: color integer :: outfile logical :: plot_on, err ! print *, 'amplitude_scan' !TEMP outfile = lunget() open(outfile, file = 'amplitude_scan.dat') write(outfile, '(9a15)') 'sqrt(Ja/emit_a)', 'Jz/sigma_pz', 'n_lost', & 'vec(1)', 'vec(2)', 'vec(3)', 'vec(4)', 'vec(5)', 'vec(6)' ix_branch = 0 call reallocate_coord(orbit, lat) ele0 => lat%ele(0) !max_turns = 1000 emit_a = normal_mode%a%emittance !emit_b = normal_mode%b%emittance emit_b = 0.01*emit_a sigma_z = normal_mode%sig_z sigma_pz = normal_mode%sigE_E !dJa = emit_a !dJb = emit_b !dJz = sigma_pz max_Ja = n_sigma_a_max**2*emit_a max_Jb = n_sigma_b_max**2*emit_b max_Jz = n_sigma_z_max**2*sigma_pz nJa = nint(n_sigma_a_max/step_a) nJb = nint(n_sigma_b_max/step_b) nJz = nint(n_sigma_z_max/step_z) print *, 'nJa, nJb, nJz: ', nJa, nJb, nJz print *, 'n_sigma_a_max, n_sigma_b_max, n_sigma_z_max', n_sigma_a_max, n_sigma_b_max, n_sigma_z_max if (plot_on) then call qp_clear_page () call qp_set_box(1, 1, 2, 1) call qp_calc_and_set_axis ('X', 0.0_rp, n_sigma_a_max, 4, 8, 'GENERAL') call qp_calc_and_set_axis ('Y', 0.0_rp, n_sigma_z_max, 4, 8, 'GENERAL') call qp_draw_axes ('sqrt(J\da\u/\ge\da\u)', '\gd\d0\u / \gs\d\gd\u (1)' )!'J\db\u/\ge\db\u') allocate(x_arr(max_turns), y_arr(max_turns)) endif !nz = 0 do nb = 0, 0 do nz = 0, nJz do na = 0, nJa ! Set initial orbit Ja=(na**2)*emit_a * step_a**2 Jb=(nb**2)*emit_b * step_b**2 Jz=nz*sigma_pz * step_z call convert_coords('LAB', orb0, lat%ele(0),'ACTION-ANGLE', orb_J_phi, err) orb_J_phi%vec = orb_J_phi%vec + [Ja, 0.0_rp, Jb , 0.0_rp, 0.0_rp, 0.0_rp] call convert_coords('ACTION-ANGLE', orb_J_phi, lat%ele(0), 'LAB', orbit(0), err) orbit(0)%vec(6) = Jz vec0 = orbit(0)%vec write(*,'(6es15.7, 2i8)') orbit(0)%vec, na, nz !print *, orbit(0)%vec n_lost = -1 !not lost do n = 1, max_turns if(plot_on) then call convert_coords('LAB', orbit(0), lat%ele(0),'NORMALIZED', orb_J_phi, err) x_arr(n) = orb_J_phi%vec(1)!sqrt(2*orb_J_phi%vec(1))*cos(orb_J_phi%vec(2)) y_arr(n) = orb_J_phi%vec(2)!-sqrt(2*orb_J_phi%vec(1))*sin(orb_J_phi%vec(2)) endif call track_all (lat, orbit, ix_branch, track_state, err) if (err) then n_lost = n exit endif ! loop back to beginning orbit(0) = orbit(lat%n_ele_track) end do write(*,'(6es15.7, 2i8)') orbit(0)%vec, n, nz !print *, orbit(0)%vec if (plot_on) then call qp_set_box(2, 1, 2, 1) call qp_clear_box() call qp_calc_and_set_axis ('X', -2*sqrt(2*max_Ja), 2*sqrt(2*max_Ja), 4, 8, 'GENERAL') call qp_calc_and_set_axis ('Y', -2*sqrt(2*max_Ja), 2*sqrt(2*max_Ja), 4, 8, 'GENERAL') call qp_draw_axes ('a_bar (sqrt(m))', "a_bar' (sqrt(m))") call qp_draw_symbols (x_arr(1:n), y_arr(1:n), type = dot_sym$) endif if (plot_on) then call qp_set_box(1, 1, 2, 1) call qp_calc_and_set_axis ('X', 0.0_rp, n_sigma_a_max, 4, 8, 'GENERAL') call qp_calc_and_set_axis ('Y', 0.0_rp, n_sigma_z_max, 4, 8, 'GENERAL') if (err) then color = red$ else color = green$ endif call qp_draw_symbol(sqrt(Ja/emit_a), Jz/sigma_pz, color=color) write(outfile, '(2es15.7, i8, 6es15.7)') sqrt(Ja/emit_a), Jz/sigma_pz, n_lost, vec0 endif end do end do end do close(outfile) end subroutine amplitude_scan !--------------------------------------------------------------------------- !+ ! ! !- !--------------------------------------------------------------------------- subroutine tunes_from_tracking(lat, orb0, n_turns, Qa, Qb, error) use fourier_mod implicit none type(lat_struct) :: lat type(coord_struct), allocatable :: orbit(:) type(coord_struct) ::orb0, orb_J_phi type(ele_struct), pointer :: ele0 real(rp) :: Qa, Qb real(rp), allocatable :: vec_set(:,:) integer :: n, n_turns integer :: ix_branch, track_state logical, optional :: error logical :: err character(30), parameter :: r_name = 'tunes_from_tracking' ! err = .false. if (present(error)) error = .false. ix_branch = 0 call reallocate_coord(orbit, lat) orbit(0) = orb0 ele0 => lat%ele(0) ix_branch = 0 allocate(vec_set(n_turns, 6)) !write (*, '(6es15.7)') orbit(0)%vec do n = 1, n_turns call track_all (lat, orbit, ix_branch, track_state, err) if (err) then print *, 'particle lost at n: ', n, track_state if (present(error)) error = err exit endif ! loop back to beginning orbit(0) = orbit(lat%n_ele_track) vec_set(n, :) = orbit(0)%vec end do if (.not. err) n = n_turns ! Fourier analysis (assumes n_turns is a power of 2) ! Remove the average Qa = fine_frequency_estimate(vec_set(1:n,1) - sum(vec_set(1:n,1))/n) Qb = fine_frequency_estimate(vec_set(1:n,3) - sum(vec_set(1:n,3))/n) ! Return last orb orb0 = orbit(0) ! Cleanup deallocate(vec_set, orbit) end subroutine !--------------------------------------------------------------------------- !+ ! Subroutine plot_particle_phase_space(particles, dim1, dim2, scale1, scale1) ! ! ! ! !- !--------------------------------------------------------------------------- subroutine plot_particle_phase_space(particles, dim1, dim2, scale1, scale2, label1, label2) implicit none type (coord_struct) :: particles(:) real(rp) :: scale1, scale2 integer :: dim1, dim2 character(*) :: label1, label2 ! !call qp_clear_box() call qp_calc_and_set_axis ('X', -scale1, scale1, 4, 8, 'GENERAL') call qp_calc_and_set_axis ('Y', -scale2, scale2, 4, 8, 'GENERAL') call qp_draw_axes ( trim(label1), trim(label2) ) call qp_draw_symbols (particles(:)%vec(dim1), particles(:)%vec(dim2), type = dot_sym$) end subroutine end module ring_aperture_mod