program osc_track use bmad use beam_mod use mode3_mod use random_mod use bsim_cesr_interface use bsim_interface use nr use sim_utils use bmadz_interface use z_tune_mod use transfer_map_mod use sim_utils_interface use histogram_mod use und_energy_trans_mod use ele_noise_mod implicit none type (lat_struct) lat, lat_ideal type (ele_struct) ele type (coord_struct), allocatable :: orbit(:), co(:) type (beam_init_struct) beam_init type (bunch_struct) bunch type (normal_modes_struct) mode type (rad_int_all_ele_struct) rad_int type (bunch_params_struct) bunch_params, bunch_params2 type (ele_pointer_struct), allocatable :: eles(:) type (taylor_struct) t_map(6) integer number, m, i, ix, nargs,j,i_turn integer ix_cache/0/ integer lun, ios integer n_turns, n_particle integer k,l integer track_clock integer n_arg integer ix_branch integer n_obs_pts/1.0/, turn_step/50000/ integer, allocatable :: obs_ele_idx(:) integer i_steps/1/, i_ix integer track_specie, track_direction integer idx real(rp) bunch_charge, bunch_current real(rp) density/1.e12/, sigx/0.01/, sigy/0.002/, sigz/0.01/ real(rp) initial_offset(6), harvest(6) real(rp) track_time real(rp) i_0/0.125/ real(rp) actual_time, delta_time real(rp) x_emit, y_emit, z_emit, ratio real(rp) Qx_hz/0.0/, Qy_hz/0.0/, Qz_hz/0.0/, frev/390.136/, target_tunes(3) real(rp), allocatable :: dk1(:) character*120 lat_file, init_dist_file character*120 file_name character*20 :: obs_ele_names(1:10) character*5 i_char logical error/.false./, excitation_on/.true./, damping_on/.true./ logical random_offset/.false./ logical at_obs_pt/.false./ logical set_tunes/.false./, ok logical use_init_dist logical output_last_turn ! Feedback definitions logical feedback_on/.false./ character*20 feed_pick, feed_kick integer feed_pick_ele, feed_kick_ele real(rp) feedback_gain/1.0/ real(rp) sz_slice/1E-3/ ! longitudinal kick real(rp) m51, m52, m56, mu0/2.405/, m56_t, delta_s, delta_s0, delta_ds real(rp), allocatable :: pick_vec(:,:) real(rp) k_light, lambda, d_s_curv, d_s_lin, d_s, d_phi, zp_kick real(rp) gamma, n_und, lambda_und, K_und real(rp) trans_mat(1:6,1:6) logical use_am_factor, display_d_s ! cooling range parameters calc real(rp) beta_p, alpha_p, gamma_p, eta_p, etap_p, phi_p, alphaz_p/0.0/, gammaz_p/0.1/ real(rp) beta_k, alpha_k, gamma_k, eta_k, etap_k, phi_k real(rp) emit_max, sig_p_max ! cooling boundary real(rp) lambda_x, lambda_s ! cooling rate real(rp) delta_emit, d_e_factor, d_f real(rp) a_s ! sample lengtherning real(rp) theta_xt, theta_zt, theta_xc real(rp) e_x, ax, az, xi_max real(rp) e_x_factor, ax_factor, az_factor ! Emittance calc definitions real(rp) sigma_mat(1:6,1:6) real(rp) normal(1:3) real(rp), ALLOCATABLE :: centroid_position(:,:) ! Radiation and fluctuation definitions real(rp) rannum(3) logical turn_on_rf/.true./ integer track_start, track_end, bypass_first_ele, bypass_last_ele ! incoherent kick integer, allocatable :: z_idx(:) integer i_sort, inco_count real(rp) inco_kick/0.0/, d_psi/0.0/ real(rp) raninco, sigma_inco, sigma_cons, sigma0 logical add_incoherent_kick/.false./ logical output_co_inco_kick/.false./, save_init_distr/.false./ logical output_tbt_vec/.false./, output_emit/.true./ logical use_dist_for_inco/.false./, z_gauss_dist/.true./ real(rp) binsize, slice_part, n_tot_part/1.0E9/ real(rp), allocatable :: z(:) integer, allocatable :: h(:) real(rp)s_model/0.0/ logical output_ds_vec/.false./ ! taylor map elements and theta real(rp) t511, t521, t522, t533, t543, t544 real(rp) t555, t565, t566 real(rp) ds_x, ds_y, ds_z real(rp) theta2_x, theta2_y, f1, f2, theta2_z real(rp) phi0_a, phi0_b logical calc_nonlin_ds/.true./ ! use real undulator radiation real(rp) delta_s_t, light_amplification_ratio/1.0/, energy_kick/0.0/, part_energy type (ele_struct), pointer :: source_ele, kick_ele type (ele_struct), pointer :: ele1, ele2 type (coord_struct), allocatable :: source_start(:), kick_start(:) logical und_track_energy_transfer/.false./ integer source_idx_start, kick_idx_start ! check particles action real(rp) jx logical output_action/.false./ real(rp) emitx_cutoff type (twiss_struct) twiss_obs ! longitudinal kick files character*120 long_kick_file real(rp) kickbuffer(1:2200,1:3) real(rp), allocatable :: longk_ds(:), longk_zp(:) logical use_kick_file/.false./, use_real_kick/.false./ ! transverse Ex file character*120 tran_kick_file integer ex_m/101/, ex_n/101/ real(rp) exbuffer(1:100000,1:3) real(rp), allocatable :: ex_x(:), ex_y(:), ex(:,:), ex_buffer(:,:) real(rp) x_light, y_light, ex_ratio ! apply magnet noise type (ele_noise_struct) :: ele_noise(n_ele_noise_max) logical apply_mag_noise character*20 bend1_name, bend2_name integer bend1_idx, bend2_idx integer n_apply_noise/2/ ! at which turn to apply noise delta_s_t = 6.695854354104528E-012 * c_light !~70meV obs_ele_names(1)='BEGINNING' track_specie = 1 track_direction = 1 use_init_dist = .false. initial_offset(1:6) = 0.0 x_emit = -1 y_emit = -1 z_emit = -1 long_kick_file ="/home/sw565/sw565/test_with_bmad/osc_und/und/dz_kick_helical.txt" tran_kick_file ="/home/sw565/sw565/test_with_bmad/osc_track/pycode/E_array_new.txt" feed_pick = "p_mark" feed_kick = "k_mark" lambda = 0.8E-6 k_light = 2*pi/lambda n_und = 4 ! =4 for cesr =7 for iota lattice use_am_factor = .false. display_d_s = .false. output_last_turn = .false. sigma_cons=251 ! for n=127661 particles within [-n_u*lambda, 0] with 1E9 particles in a bunch, n_u=4, labmda=0.8um ! iota lattice: 12.2 for n_slice=300 partices with 1E7 particles in a bunch, n_u=7, lambda=2.2E-6m sigma0 = 10.0E-3 ! bunch length at which the sigma_cons was calculated bypass_first_ele = 0 binsize=0.5E-3 bend1_name = "b12w" bend2_name = "b46w" ! input list namelist /input/ n_turns, n_particle, i_0, excitation_on, damping_on, turn_on_rf, feedback_on, sz_slice, & initial_offset, random_offset, lat_file, turn_step, x_emit, y_emit, z_emit, n_obs_pts, obs_ele_names, & set_tunes, Qx_hz, Qy_hz, Qz_hz, track_specie, feedback_gain, use_init_dist, & init_dist_file, feed_pick, feed_kick, output_last_turn, output_co_inco_kick, use_am_factor, display_d_s, & add_incoherent_kick, output_tbt_vec, output_emit, use_dist_for_inco, sigma_cons, output_ds_vec, z_gauss_dist, & binsize,n_tot_part, lambda, n_und, sigma0, und_track_energy_transfer, save_init_distr, output_action, & use_kick_file, long_kick_file, tran_kick_file, use_real_kick, calc_nonlin_ds, apply_mag_noise, ele_noise, & bend1_name, bend2_name, n_apply_noise ! read in the paramter file n_arg = cesr_iargc() if (n_arg > 1) then print *, 'Usage: osc_track ' print *, 'Default: = osc_track.in' stop endif if (n_arg == 1) call cesr_getarg(1, file_name) if (n_arg == 0) then print '(a, $)', ' Input command file : ' read (*, '(a)') file_name endif if (file_name == '') file_name = 'osc_track.in' print *, 'Opening: ', trim(file_name) open (unit= 1, file = file_name, status = 'old') read (1, nml = input) close (unit = 1) write(6, nml = input) ! read longitudinal kick file if (use_kick_file) then open(unit=2, file=long_kick_file, status = 'old', iostat = ios) if(ios .ne. 0) then write(*,*) 'Error: cannot open file: ', trim(long_kick_file) stop else i = 1 do read(2, *, iostat=ios) kickbuffer(i,:) if(ios > 0) then write(*,*) "error: cannot read in from file: ", trim(long_kick_file) stop elseif(ios < 0) then write(*,*) "end of longitudinal kick file reached" exit endif i = i + 1 enddo end if close(2) allocate(longk_ds(1:i-1)) allocate(longk_zp(1:i-1)) longk_ds(:) = kickbuffer(1:i-1,2) longk_zp(:) = kickbuffer(1:i-1,3) print *, "max ds: ", maxval(longk_ds), " min ds: ", minval(longk_ds) print *, "max energy transfer (eV): ", maxval(longk_zp) open(unit=2, file=tran_kick_file, status = 'old', iostat = ios) if(ios .ne. 0) then write(*,*) 'Error: cannot open file: ', trim(tran_kick_file) stop else i = 1 do read(2, *, iostat=ios) exbuffer(i,:) if(ios > 0) then write(*,*) "error: cannot read in from file: ", trim(tran_kick_file) stop elseif(ios < 0) then write(*,*) "end of transverse Ex file reached" exit endif i = i + 1 enddo end if close(2) allocate(ex_buffer(1:i-1,1:3)) allocate(ex_x(1:ex_m)) allocate(ex_y(1:ex_n)) ex_buffer(:,1)=exbuffer(1:i-1,1) ex_buffer(:,2)=exbuffer(1:i-1,2) ex_buffer(:,3)=exbuffer(1:i-1,3) call extract_2darrays_to_arrays(ex_buffer, ex_m, ex_n, ex_x, ex_y, ex) end if k_light = 2*pi/lambda ! record start time call tick(track_clock) ! prepare lattice bmad_com%radiation_damping_on = .false. bmad_com%radiation_fluctuations_on = .false. if (excitation_on) then bmad_com%radiation_fluctuations_on = .true. end if if (damping_on) then bmad_com%radiation_damping_on = .true. end if bmad_com%max_aperture_limit = .045 bmad_com%csr_and_space_charge_on = .true. call bmad_parser (lat_file, lat) print *,' bmad_parser done' allocate(centroid_position(lat%n_ele_track,1:6)) allocate(pick_vec(n_particle,6)) allocate(z_idx(n_particle)) allocate(source_start(n_particle)) allocate(kick_start(n_particle)) ! Check tracking species if ( track_specie .eq. 1 ) then lat%param%particle = positron$ !set the ref charge to be positron, default tracking charge should be the same. track_direction = 1 else lat%param%particle = electron$ track_direction = -1 endif ! locate observation points allocate(obs_ele_idx(1:n_obs_pts)) write (*,'(a)') 'Find observation element index:' write (*,'(3a12)') 'obs', 'Index', 'Name' ! for real undulator energy transfer calc source_ele => pointer_to_ele(lat, 'osc_pickup') call find_element_ends(source_ele, ele1, ele2) source_idx_start = ele1%ix_ele print *, 'source_idx_start: ', source_idx_start kick_ele => pointer_to_ele(lat, 'osc_kicker') call find_element_ends(kick_ele, ele1, ele2) kick_idx_start = ele1%ix_ele print *, 'kick_idx_start: ', kick_idx_start do i=1, n_obs_pts call lat_ele_locator(obs_ele_names(i),lat,eles,l,error) obs_ele_idx(i) = eles(1)%ele%ix_ele write(*,'(2i12, a24)')i, obs_ele_idx(i), obs_ele_names(i) end do ! locate the OSC feedback pickup and kicker element (undulators) call lat_ele_locator(feed_pick, lat, eles, l, error) feed_pick_ele=eles(1)%ele%ix_ele call lat_ele_locator(feed_kick, lat, eles, l, error) feed_kick_ele=eles(1)%ele%ix_ele ! if pick and kick elements are same if (feed_pick_ele .eq. feed_kick_ele) then feed_pick_ele = feed_pick_ele -1 end if write(*,'(2i12, a24)') i, feed_pick_ele, feed_pick write(*,'(2i12, a24)') i, feed_kick_ele, feed_kick ! find bend1 and bend2 indices for checking field call lat_ele_locator(bend1_name, lat, eles, l, error) bend1_idx=eles(1)%ele%ix_ele call lat_ele_locator(bend2_name, lat, eles, l, error) bend2_idx=eles(1)%ele%ix_ele i = index(file_name,'.') ! open(unit=28,file=file_name(1:i-1)//'.out') if (apply_mag_noise) then open(unit=28, file='bfield_'//file_name(1:i-1)//'.out') end if bunch_current=i_0 if (output_tbt_vec) then open(unit=24,file=file_name(1:i-1)//'vec_tbt.out') end if if (output_action) then open(unit=101,file=file_name(1:i-1)//'_action.out') end if ! output tracking status after finish open(unit=25,file='particle_status_'//file_name(1:i-1)//'.out') write(25,'(5a12)') 'Turn', 'Total N', 'Live N','Lost ele','Live Q' bunch_charge = bunch_current*1E-3*lat%param%total_length/c_light frev = c_light/lat%param%total_length/1E3 ! in kHz call set_on_off(rfcavity$, lat, on$) call reallocate_coord (orbit, lat%n_ele_track) call closed_orbit_calc(lat,orbit,6) call track_all(lat,orbit) call lat_make_mat6(lat, -1, ref_orb = orbit) call twiss_at_start(lat) call twiss_propagate_all(lat) call radiation_integrals (lat, orbit, mode, ix_cache, 0, rad_int) call calc_z_tune(lat) if (Qx_hz .eq. 0.0) then Qx_hz = lat%a%tune/twopi*frev endif if (Qy_hz .eq. 0.0) then Qy_hz = lat%b%tune/twopi*frev end if if (Qz_hz .eq. 0.0) then Qz_hz = lat%z%tune/twopi*frev end if ! set tunes if (set_tunes) then write (*,'(a33,f12.6,a9,f12.6,a9,f12.6)') 'Current tunes (kHz): Qx = ',lat%a%tune/twopi*frev, & ' Qy = ',lat%b%tune/twopi*frev,' Qz = ', lat%z%tune/twopi*frev write (*,'(a28,f12.6,a9,f12.6,a9,f12.6)') 'Target tunes : Qx = ',Qx_hz, ' Qy = ',Qy_hz,' Qz = ', Qz_hz target_tunes(1)=(Qx_hz/frev+int(lat%ele(lat%n_ele_track)%a%phi/twopi))*twopi ! integer+fractional tunes target_tunes(2)=(Qy_hz/frev+int(lat%ele(lat%n_ele_track)%b%phi/twopi))*twopi target_tunes(3)=Qz_hz/frev*twopi call choose_quads(lat, dk1) call set_tune(target_tunes(1), target_tunes(2), dk1, lat, orbit, ok) call set_z_tune(lat, target_tunes(3)) else print *, 'Use lattice design tunes!' print *, lat%a%tune/twopi, lat%b%tune/twopi,lat%z%tune/twopi end if call reallocate_coord (co, lat%n_ele_track) co = orbit call longitudinal_beta(lat, mode) call twiss3_at_start(lat, error) call twiss3_propagate_all(lat) call twiss_at_start(lat) call twiss_propagate_all(lat) ! turn on or off RF if (.not. turn_on_rf) then call set_on_off(rfcavity$, lat, off$) end if track_start=1 track_end=lat%n_ele_track part_energy = lat%ele(1)%value(e_tot$) if (und_track_energy_transfer) light_amplification_ratio = feedback_gain * 1E9 /0.07 if (use_kick_file) then light_amplification_ratio = feedback_gain * 1E9 / maxval(longk_zp) if (use_real_kick) then feedback_gain = maxval(longk_zp)/part_energy light_amplification_ratio = 1.0 end if end if ! output twiss parameters with wakes off open(13,file='twiss.out') write(13, '(13a12)') "# ele", "s", "phi_x", "beta_x", "alpha_x", "phi_y", "beta_y", "alpha_y", "eta_x", "etap_x", "eta_y", "etap_y", "beta_z" do i=1, lat%n_ele_track write(13,'(i12,12es12.4)') i, lat%ele(i)%s, lat%ele(i)%a%phi, lat%ele(i)%a%beta, lat%ele(i)%a%alpha, lat%ele(i)%b%phi, lat%ele(i)%b%beta, & lat%ele(i)%b%alpha, lat%ele(i)%x%eta, lat%ele(i)%x%etap, lat%ele(i)%y%eta, lat%ele(i)%y%etap, lat%ele(i)%z%beta enddo close(13) open(14,file='twiss3.out') write(14, '(14a12)') "# ele", "s", "phi_x", "beta_x", "alpha_x", "phi_y", "beta_y", "alpha_y", "eta_x", "etap_x", "eta_y", "etap_y", "beta_c", "gamma_c" do i=1, lat%n_ele_track write(14,'(i12,13es12.4)') i, lat%ele(i)%s, lat%ele(i)%mode3%a%phi, lat%ele(i)%mode3%a%beta, lat%ele(i)%mode3%a%alpha, & lat%ele(i)%mode3%b%phi, lat%ele(i)%mode3%b%beta, lat%ele(i)%mode3%b%alpha, & lat%ele(i)%mode3%x%eta, lat%ele(i)%mode3%x%etap, lat%ele(i)%mode3%y%eta, lat%ele(i)%mode3%y%etap, lat%ele(i)%mode3%c%beta, lat%ele(i)%mode3%c%gamma enddo close(14) ! Output closed orbit open(12,file='closed_orbit.out') do i = 1, lat%n_ele_track write(12, '(i12,7es12.4)') i, orbit(i)%vec(:), lat%ele(i)%s end do close(12) print '(3(a,f12.6), a6)','horizontal tune = ', lat%a%tune/twopi*frev,' vertical tune = ', & lat%b%tune/twopi*frev,' synch tune = ',lat%z%tune/twopi*frev, ' kHz' print '(2(a23,es12.4))', 'horizontal emittance = ',mode%a%emittance,'vertical emittance = ', mode%b%emittance ! calculate the trasfer matrix from pickup to kicker i=feed_pick_ele+1 trans_mat(:,:) = lat%ele(i)%mat6(:,:) if (feed_kick_ele .lt. feed_pick_ele) then print *, 'Kicker element index smaller then pickup element index!!!' stop else if (feed_kick_ele - 1 .eq. feed_pick_ele) then trans_mat(:,:) = lat%ele(feed_kick_ele)%mat6(:,:) else do i=feed_pick_ele+2, feed_kick_ele trans_mat(:,:) = matmul(lat%ele(i)%mat6(:,:), trans_mat(:,:)) end do end if m51=trans_mat(5,1) m52=trans_mat(5,2) m56=trans_mat(5,6) print *,'-----' print *,'pickup to kicker transfer matrix: ' do i=1,6 print '(6es12.4)', trans_mat(i,:) end do delta_s0 = lat%ele(feed_kick_ele)%s-lat%ele(feed_pick_ele)%s -abs(lat%ele(feed_kick_ele)%floor%r(3)-lat%ele(feed_pick_ele)%floor%r(3)) s_model= lat%ele(feed_kick_ele)%s - lat%ele(feed_pick_ele)%s_start m56_t = m51*lat%ele(feed_pick_ele)%a%eta + m52*lat%ele(feed_pick_ele)%a%etap + m56 print *,'-----' print '(a,es12.4)', 'Bypass delay: ', delta_s0 print '(4(a,es12.4))', 'm51: ', m51, ' m52: ', m52, ' m56: ', m56, ' m56_t: ', m56_t print *,'-----' ! Calculate the nonlinear bypass delay using T elements from Taylor maps if (calc_nonlin_ds) then call transfer_map_calc(lat, t_map, error, ix1 = feed_pick_ele, ix2 = feed_kick_ele, concat_if_possible = .true.) t511 = taylor_coef(t_map(5),taylor_expn([1,1])) t521 = taylor_coef(t_map(5),taylor_expn([2,1])) t522 = taylor_coef(t_map(5),taylor_expn([2,2])) t533 = taylor_coef(t_map(5),taylor_expn([3,3])) t543 = taylor_coef(t_map(5),taylor_expn([4,3])) t544 = taylor_coef(t_map(5),taylor_expn([4,4])) t555 = taylor_coef(t_map(5),taylor_expn([5,5])) t565 = taylor_coef(t_map(5),taylor_expn([6,5])) t566 = taylor_coef(t_map(5),taylor_expn([6,6])) print '(3(a,es12.4))', 't511: ', t511, ' t521: ', t521, ' t522: ', t522 print '(3(a,es12.4))', 't533: ', t533, ' t543: ', t543, ' t544: ', t544 print '(3(a,es12.4))', 't555: ', t555, ' t565: ', t565, ' t566: ', t566 ds_x = t511*lat%ele(feed_pick_ele)%a%beta - t521*lat%ele(feed_pick_ele)%a%alpha+ t522*lat%ele(feed_pick_ele)%a%gamma ds_y = t533*lat%ele(feed_pick_ele)%b%beta - t543*lat%ele(feed_pick_ele)%b%alpha+ t544*lat%ele(feed_pick_ele)%b%gamma ds_z = t555*lat%ele(feed_pick_ele)%mode3%c%beta - t565*lat%ele(feed_pick_ele)%mode3%c%alpha+ t566*lat%ele(feed_pick_ele)%mode3%c%gamma print *,'-----' print '(3(a,es12.4))', 'nonlinear_ds_x: ', ds_x, ' nonlinear_ds_y: ', ds_y, ' nonlinear_ds_z: ', ds_z print *,'-----' end if theta2_x = 0.0 theta2_y = 0.0 theta2_z = 0.0 phi0_a=lat%ele(feed_pick_ele)%a%phi phi0_b=lat%ele(feed_pick_ele)%b%phi do i=feed_pick_ele+1,feed_kick_ele f1 = ( lat%ele(i-1)%a%alpha**2 + 1 )/lat%ele(i-1)%a%beta f2 = ( lat%ele(i )%a%alpha**2 + 1 )/lat%ele(i )%a%beta theta2_x = theta2_x + lat%ele(i)%value(l$)*(f1+f2)/2.0 ! print '(i10, 4es12.4)', i, f1, f2, lat%ele(i)%value(l$)*(f1+f2)/2.0, theta2_x f1 = ( lat%ele(i-1)%b%alpha**2 + 1 )/lat%ele(i-1)%b%beta f2 = ( lat%ele(i )%b%alpha**2 + 1 )/lat%ele(i )%b%beta theta2_y = theta2_y + lat%ele(i)%value(l$)*(f1+f2)/2.0 ! print '(i10, 4es12.4)', i, f1, f2, lat%ele(i)%value(l$)*(f1+f2)/2.0, theta2_y f1 = lat%ele(i-1)%mode3%c%gamma f2 = lat%ele(i )%mode3%c%gamma theta2_z = theta2_z + lat%ele(i)%value(l$)*(f1+f2)/2.0 end do print '(3(a,es12.4))', 'theta2_x: ', theta2_x, ' theta2_y: ', theta2_y, ' theta2_z: ', theta2_z print *,'-----' ! Twiss parameters at the first observation point twiss_obs%beta=lat%ele(obs_ele_idx(1))%a%beta twiss_obs%alpha=lat%ele(obs_ele_idx(1))%a%alpha twiss_obs%gamma=lat%ele(obs_ele_idx(1))%a%gamma twiss_obs%eta=lat%ele(obs_ele_idx(1))%a%eta twiss_obs%etap=lat%ele(obs_ele_idx(1))%a%etap ! Twiss parameters at the OSC undulators beta_p=lat%ele(feed_pick_ele)%a%beta alpha_p=lat%ele(feed_pick_ele)%a%alpha gamma_p=lat%ele(feed_pick_ele)%a%gamma eta_p=lat%ele(feed_pick_ele)%a%eta etap_p=lat%ele(feed_pick_ele)%a%etap phi_p=lat%ele(feed_pick_ele)%a%phi ! alphaz_p=lat%ele(feed_pick_ele)%z%alpha alphaz_p=lat%ele(feed_pick_ele)%mode3%c%alpha gammaz_p=lat%ele(feed_pick_ele)%mode3%c%gamma beta_k=lat%ele(feed_kick_ele)%a%beta alpha_k=lat%ele(feed_kick_ele)%a%alpha gamma_k=lat%ele(feed_kick_ele)%a%gamma eta_k=lat%ele(feed_kick_ele)%a%eta etap_k=lat%ele(feed_kick_ele)%a%etap phi_k=lat%ele(feed_kick_ele)%a%phi d_f=phi_k-phi_p theta_xt=atan2(m52*alpha_p-m51*beta_p, m52) theta_zt=atan2(alphaz_p, 1.0_rp) theta_xc= - atan2(eta_k, beta_k*etap_k+alpha_k*eta_k) e_x_factor = eta_k**2*gamma_k+beta_k*etap_k**2+2*alpha_k*eta_k*etap_k ax_factor= sqrt(m51**2*beta_p+m52**2*gamma_p-2*alpha_p*m51*m52) az_factor= m56_t ! az_factor= sqrt(gammaz_p)*m56_t ! seems wrong write(*,'(3(a12, es12.4))') 'ex_factor:', e_x_factor, 'ax_factor:',ax_factor,'az_factor:',az_factor write(*,'(3(a12, es12.4))') 'theta_xt:', theta_xt, 'theta_zt:', theta_zt, 'theta_xc:',theta_xc print *,'-----' ! print the essential quantities print *,'Pickup undulator Twiss parameters: ' write(*,'(3(a10,es12.4))') 'beta_x: ', beta_p , 'alpha_x: ', alpha_p, 'gamma_x:', gamma_p write(*,'(3(a10,es12.4))') 'eta_x: ', eta_p, 'etap_x: ', etap_p, 'phi_x:', phi_p write(*,'(3(a10,es12.4))') 'alpha_z: ',alphaz_p, 'gamma_z: ', gammaz_p print *,'-----' print *,'Kicker undulator Twiss parameters: ' write(*,'(3(a10,es12.4))') 'beta_x: ', beta_k , 'alpha_x: ', alpha_k, 'gamma_x:', gamma_k write(*,'(3(a10,es12.4))') 'eta_x: ', eta_k, 'etap_x: ', etap_k,'phi_x:', phi_k print *,'-----' ! sample lengtherning bypass line: a_s = sqrt(mode%a%emittance*ax_factor**2)+m56_t*mode%sige_e print *,'sample lengtherning amplitude: ' write(*,'(a, es12.4)') 'a_s: ', a_s print *,'-----' if (ax_factor .ne. 0.0) then emit_max = (mu0/k_light)**2/ax_factor**2 end if ! sig_p_max = mu0/k_light/m56_t sig_p_max = mu0/k_light/az_factor print *,'Cooling range: ' write(*,'((a10, es12.4, a15, es12.4))') 'emit_max: ', emit_max, 'sig_p_max: ', sig_p_max write(*,'((a10, es12.4, a15, es12.4))') 'n_x: ', sqrt(emit_max/mode%a%emittance), 'n_z: ', sig_p_max/mode%sigE_E print *,'-----' lambda_x = k_light/2*(m56-m56_t)*feedback_gain lambda_s = k_light/2*m56_t*feedback_gain print *,'Cooling rates: ' write(*,'((a10, es12.4, a15, es12.4))') 'lambda_x:',lambda_x , 'lambda_s: ', lambda_s write(*,'(a,es12.4)') 'Cooling ratio x/z: ', lambda_x/lambda_s print *,'-----' write(*,'(a, f12.4)') 'phase advance: ', d_f ! change of emittance every turn: d_e_factor= m51*sqrt(beta_p/beta_k)*( eta_k*(cos(d_f)-alpha_k*sin(d_f))-etap_k*beta_k*sin(d_f) ) + & m52*eta_k/sqrt(beta_k*beta_p)*( (alpha_k-alpha_p)*cos(d_f) + (1+alpha_k*alpha_p)*sin(d_f) ) + & m52*etap_k*sqrt(beta_k/beta_p)*( cos(d_f)+alpha_p*sin(d_f) ) delta_emit = 2*pi*feedback_gain*k_light*d_e_factor*x_emit write(*,'(2(a,es12.4))') 'emit_factor: ',d_e_factor, ' delta_emit: ', delta_emit ! setup bunch parameters if( x_emit .gt. 0.0 ) then mode%a%emittance = x_emit endif if( y_emit .gt. 0.0 ) then mode%b%emittance = y_emit endif if( z_emit .gt. 0.0 ) then ratio = mode%sig_z / mode%sige_e mode%sige_e = sqrt(z_emit / ratio) mode%sig_z = sqrt(ratio * z_emit) endif beam_init%n_particle= n_particle beam_init%bunch_charge = bunch_charge beam_init%a_emit = mode%a%emittance beam_init%b_emit = mode%b%emittance beam_init%sig_z = mode%sig_z beam_init%sig_e = mode%sigE_E beam_init%full_6D_coupling_calc = .false. call init_bunch_distribution (lat%ele(bypass_first_ele), lat%param, beam_init, 0, bunch) if (apply_mag_noise) then call init_ele_noise(lat, ele_noise) lat_ideal = lat ! Save the ideal no-noise lattice end if ! read initial distribution if true if (use_init_dist) then open(unit=101, file=init_dist_file) do j=1, n_particle read(101, '(i14, 6es14.4)') idx, bunch%particle(j)%vec(:) end do close(101) end if print '(2(a23,es12.4))', 'horizontal emittance = ',beam_init%a_emit,'vertical emittance = ', beam_init%b_emit print '(2(a23,es12.4))', 'bunch length = ',beam_init%sig_z,'sig E/E = ', beam_init%sig_e print '(a23,es12.4)','bunch charge =', beam_init%bunch_charge print '(1(a23,i10))','number of bunches =', beam_init%n_bunch print '(1(a23,i10))', 'particles/bunch =', beam_init%n_particle print '(1(a23,i10))', 'number of turns =', n_turns print '(a,6es12.4)','initial_offset =',initial_offset(1:6) print '(a,es12.4)','Energy: ', lat%ele(1)%value(e_tot$) write(*,'(a,i10)') 'Number of track elements: ', lat%n_ele_track write(*,'(a,i10)') 'Number of max elements: ', lat%n_ele_max write(*,'(a, es12.4)') 'BMAD max aperture limit', bmad_com%max_aperture_limit write(*,'(a, i5)') 'Tracking reference particle: ', lat%param%particle write(*,'(a, i5)') 'Tracking particle charge: ', lat%param%default_tracking_species write(*,'(a, a10)') 'Beam init struct specie: ', beam_init%species write(*, '(1(a23,L1))') 'radiation_damping_on = ', bmad_com%radiation_damping_on write(*, '(1(a28,L1))') 'radiation_fluctuations_on = ', bmad_com%radiation_fluctuations_on if (use_kick_file .and. use_real_kick) then write(*,'(2(a28,es12.4))') 'Passive gain: ', maxval(longk_zp)/part_energy, 'Feedback gain: ', feedback_gain end if ! first turn change: e_x=sqrt(beam_init%a_emit)*e_x_factor ax= sqrt(beam_init%a_emit)*ax_factor az= sqrt(beam_init%sig_z*beam_init%sig_e)*az_factor d_e_factor = 2*e_x*bessel_jn(0,k_light*az)*sqrt(2.0)*sin(theta_zt+pi/4)*& bessel_jn(1,k_light*ax)*cos(d_f+theta_xc-theta_xt) delta_emit = feedback_gain*d_e_factor xi_max = -beam_init%a_emit/d_e_factor*bessel_jn(0,k_light*az) write(*,'(a)') '--------------------' write(*,'(2(a,es12.4))') 'd_e_factor: ', d_e_factor, ' maxiume xi: ', xi_max write(*,'(1(a,es12.4))') 'First turn emit change: ', delta_emit ! add random bunch jitter ! call ran_seed_put(101) ! debug if (.not. use_init_dist) then if (random_offset) then call ran_uniform(harvest) ! generate random numbers between 0 and 1 harvest=harvest*2-1 ! generate random numbers between -1 and 1 forall(i=1:6) bunch%particle(:)%vec(i) = orbit(0)%vec(i)+bunch%particle(:)%vec(i) + initial_offset(i)*harvest(i) else forall(i=1:6) bunch%particle(:)%vec(i) = orbit(0)%vec(i) + bunch%particle(:)%vec(i) + initial_offset(i) endif end if j = index(file_name,'.') ! Write out initial particle distribution if (save_init_distr) then open(unit=23,file='init_distr_'//file_name(1:j-1)//'.out') do i=1,n_particle write(23,'(i14,6es14.4)') i, bunch%particle(i)%vec(:) enddo close(23) end if ! Write out last initial particle distribution if (output_last_turn) then open(unit=27,file='last_distr_'//file_name(1:j-1)//'.out') end if if(output_co_inco_kick) then open(unit=29,file='co_inco_kick_'//file_name(1:j-1)//'.out') end if if(output_ds_vec) then open(unit=30, file='ds_vec_'//file_name(1:j-1)//'.out') end if ! Open observation point files do i=1, n_obs_pts j = index(file_name,'.') open(unit=(10000 + obs_ele_idx(i)), file=file_name(1:j-1)//'.out') open(unit=(20000 + obs_ele_idx(i)), file=file_name(1:j-1)//'_2.out') write(10000 + obs_ele_idx(i), '(9a12)') "# turn", "emit_a", "emit_b", "emit_c", "sqrt()", "sqrt()", & "sqrt()", "sqrt()", "delta_ex" write(20000 + obs_ele_idx(i), '(8a12)') "# turn", "emit_a", "emit_b", "emit_c", "sqrt()", "sqrt()", & "sqrt()", "sqrt()" enddo ! track multiple particles in a single bunch write(*, '(2(a15,i6))') 'start ele idx: ', track_start, 'end ele idx: ', track_end if (display_d_s) then print '(2a8,8a12)', 'i','j','ds_lin', 'ds', 'ds/nlambda', 'd_phi', 'z_pick', 'z_kick', 'zp_kick','sigma_inco' end if ! open(unit=111,file=file_name(1:i-1)//'_emit.out') ! open(unit=111,file=file_name(1:i-1)//'_action.out') do i_turn=1,n_turns if (apply_mag_noise) then if (i_turn .ge. n_apply_noise ) then call apply_ele_noise(lat, lat_ideal, ele_noise, i_turn) if (turn_on_rf) then call twiss3_at_start(lat, error) call twiss3_propagate_all(lat) else call twiss_at_start(lat) call twiss_propagate_all(lat) end if end if delta_s = lat%ele(feed_kick_ele)%s-lat%ele(feed_pick_ele)%s -abs(lat%ele(feed_kick_ele)%floor%r(3)-lat%ele(feed_pick_ele)%floor%r(3)) delta_ds = delta_s - delta_s0 write(28,'(i5, 5es20.10)') i_turn, lat%ele(bend1_idx)%value(g_err$), lat_ideal%ele(bend1_idx)%value(g$), & lat%ele(bend2_idx)%value(g_err$), lat_ideal%ele(bend2_idx)%value(g$), delta_ds end if do m=track_start, track_end call track1_bunch (bunch, lat, lat%ele(m), bunch, error, orbit, track_direction) if (n_particle .gt. 10) then call calc_bunch_params(bunch, bunch_params, error, .true.) centroid_position(m,:) = bunch_params%centroid%vec(:) sigma_mat = bunch_params%sigma normal(1) = bunch_params%a%emit normal(2) = bunch_params%b%emit normal(3) = bunch_params%c%emit end if ! write(111, '(2i12,3es12.4)') i_turn,m, normal(1), normal(2), normal(3) ! gamma_p=lat%ele(m)%a%gamma ! alpha_p=lat%ele(m)%a%alpha ! beta_p =lat%ele(m)%a%beta ! jx = gamma_p*bunch%particle(1)%vec(1)**2 + 2*alpha_p*bunch%particle(1)%vec(2)*bunch%particle(1)%vec(2) + beta_p*bunch%particle(1)%vec(2)**2 ! write(111,'(2i14,1es14.4)') i_turn, m, jx if (m .eq. source_idx_start) source_start(:)=bunch%particle(:) if (m .eq. kick_idx_start) kick_start(:)=bunch%particle(:) ! if at the undulator pickup if (m .eq. feed_pick_ele) then z_idx(:)=0 do j=1, n_particle pick_vec(j,:)=bunch%particle(j)%vec(:) if (output_tbt_vec) then write(24,'(2i14,6es14.4)') i_turn, j, bunch%particle(j)%vec(:) end if if (output_action) then jx = gamma_p*pick_vec(j,1)**2 + 2*alpha_p*pick_vec(j,1)*pick_vec(j,2) + beta_p*pick_vec(j,2)**2 write(101,'(2i14,1es14.4)') i_turn, j, jx end if end do ! sort and find the ascending indices call indexx(pick_vec(:,5),z_idx) if (add_incoherent_kick) then call hist_real(pick_vec(:,5),z_idx,binsize,h,z) end if end if ! check element at feedback pickup ! check if we are at an observation point at_obs_pt = .false. do j=1, n_obs_pts if (m .eq. obs_ele_idx(j)) then at_obs_pt = .true. endif enddo if( at_obs_pt ) then e_x=sqrt(normal(1))*e_x_factor ax= sqrt(normal(1))*ax_factor az= sqrt(normal(3))*az_factor d_e_factor = 2*e_x*bessel_jn(0,k_light*az)*sqrt(2.0)*sin(theta_zt+pi/4)*& bessel_jn(1,k_light*ax)*cos(d_f+theta_xc-theta_xt) delta_emit = feedback_gain*d_e_factor if (output_emit) then write((10000 + m), '(i12,8es12.4)') i_turn, normal(1), normal(2), normal(3), sqrt(sigma_mat(1,1)), & sqrt(sigma_mat(3,3)), sqrt(sigma_mat(5,5)), sqrt(sigma_mat(6,6)), delta_emit end if emitx_cutoff = emit_max call calc_bunch_params_cutoff (bunch, bunch_params2, twiss_obs, emitx_cutoff, error) sigma_mat = bunch_params2%sigma normal(1) = bunch_params2%a%emit normal(2) = bunch_params2%b%emit normal(3) = bunch_params2%c%emit if (output_emit) then write((20000 + m), '(i12,7es12.4)') i_turn, normal(1), normal(2), normal(3), sqrt(sigma_mat(1,1)), & sqrt(sigma_mat(3,3)), sqrt(sigma_mat(5,5)), sqrt(sigma_mat(6,6)) end if end if ! if at the undulator kicker if (m .eq. feed_kick_ele) then do i_ix = n_particle, 1 , -1 ! start from the largest z particle (sorted) j=z_idx(i_ix) ! original particle index if (bunch%particle(j)%state .eq. alive$) then inco_kick = 0 d_s_lin=m51*pick_vec(j,1) + m52*pick_vec(j,2) + m56*pick_vec(j,6) ! linear map (first order taylor map) d_s= bunch%particle(j)%vec(5) - pick_vec(j,5) + delta_ds d_phi=k_light*d_s ! d_phi=k_light*d_s_lin ! test if (und_track_energy_transfer) then call und_plane_kick (lat, source_start(j), kick_start(j), source_ele, kick_ele, delta_s_t, d_s, light_amplification_ratio, energy_kick) end if !calculate the incoherent kick for this particle if (add_incoherent_kick) then if (use_dist_for_inco) then if (z_gauss_dist) then sigma_inco = sigma_cons*exp(-((pick_vec(j,5)-orbit(m)%vec(5))/mode%sig_z)**2/4)*sqrt(sigma0/mode%sig_z) ! proportional to sqrt(N) else call interp1d(z,h,pick_vec(j,5),slice_part) slice_part = slice_part*n_tot_part/n_particle*n_und*lambda/binsize ! normalize sigma_inco = sigma_cons*sqrt(slice_part/127661) ! sigma_cons=250.8 for N=127661, lambda=0.8um end if call ran_gauss(raninco) inco_kick = -sign(1.0, m56)*feedback_gain*raninco*sigma_inco else i_sort = i_ix - 1 inco_count = 0 if (i_sort .gt. 0) then do while ( pick_vec(z_idx(i_sort),5)-pick_vec(z_idx(i_ix),5) .ge. -sz_slice ) d_psi=2*pi*(pick_vec(z_idx(i_ix),5)-pick_vec(z_idx(i_sort),5))/sz_slice inco_kick = inco_kick - sign(1.0, m56)*feedback_gain*sin(d_phi + d_psi) i_sort = i_sort - 1 if (i_sort .eq. 0) then exit end if inco_count = inco_count + 1 end do end if ! i_sort end if ! use distribution to calcalute inco or not end if ! add_incoherent_kick if (use_kick_file) then ! call interp1d_real(longk_ds, longk_zp, d_s, zp_kick) call interp1d_real_linear(longk_ds, longk_zp, d_s, zp_kick) call convert_pc_to((1+bunch%particle(j)%vec(6))*bunch%particle(j)%p0c,lat%param%particle, part_energy) zp_kick = zp_kick * light_amplification_ratio /part_energy x_light = bunch%particle(j)%vec(1) + pick_vec(j,1) y_light = bunch%particle(j)%vec(3) + pick_vec(j,3) call interp2d_real(ex_x, ex_y, ex, x_light, y_light, ex_ratio) zp_kick = zp_kick * ex_ratio else if (use_am_factor) then zp_kick = - sign(1.0, m56)*feedback_gain*(1-d_s/n_und/lambda)*sin(d_phi) else zp_kick = - sign(1.0, m56)*feedback_gain*sin(d_phi) end if end if if (und_track_energy_transfer) then call convert_pc_to((1+bunch%particle(j)%vec(6))*bunch%particle(j)%p0c,lat%param%particle, part_energy) zp_kick = energy_kick/part_energy end if if (add_incoherent_kick) then zp_kick = zp_kick + inco_kick end if if (display_d_s) then print '(2i8,8es12.4)', i_turn,j, d_s_lin, d_s, d_s/n_und/lambda, d_phi, pick_vec(j,5), bunch%particle(j)%vec(5), zp_kick, sigma_inco end if if (output_co_inco_kick) then write(29,'(3i8,3es15.5)') i_ix, j, inco_count, inco_kick, pick_vec(z_idx(i_ix),5), -sign(1.0, m56)*feedback_gain*sin(d_phi) end if if (output_ds_vec) then if ((mod(i_turn,turn_step) .eq. 0) .or. (i_turn .eq. 1)) then write(30,'(2i8,14es12.4)') i_turn,j,d_s_lin, d_s, pick_vec(j,:), bunch%particle(j)%vec(:) end if end if if (feedback_on .eq. .true.) then bunch%particle(j)%vec(6) = bunch%particle(j)%vec(6) + zp_kick end if end if ! check whether the particle is alive end do ! i_ix particle number in sorted order end if ! check whether at kicker end do ! m element ! Output last turn vec if (output_last_turn) then if (i_turn .eq. n_turns) then do j=1, n_particle write(27,'(i14, 6es14.4)') j, bunch%particle(j)%vec(:) enddo end if end if end do !i turn call calc_bunch_params(bunch, bunch_params, error, .true.) write(25,'(4i12,es12.4)') i_turn, bunch_params%n_particle_tot, bunch_params%n_particle_live, & bunch_params%n_particle_lost_in_ele, bunch_params%charge_live close(25) ! Write out centroid orbit for last turn ! open(unit=26,file='cent_lturn_'//i_char//'.out') ! do i=1, lat%n_ele_track ! write(26,'(i12,12es12.4)') i, centroid_position(i,1:6) ! enddo ! close(26) ! Close observation point files do i=1, n_obs_pts close((10000 + obs_ele_idx(i))) close((20000 + obs_ele_idx(i))) enddo if (output_tbt_vec) close(unit=24) if (output_last_turn) close(unit=27) if (output_co_inco_kick) close(unit=29) if (output_ds_vec) close(unit=30) if (output_action) close(unit=101) close(unit=28) ! close(unit=111) deallocate(centroid_position) deallocate(obs_ele_idx) deallocate(pick_vec) deallocate(z_idx) deallocate(source_start) deallocate(kick_start) if (use_kick_file) then deallocate(longk_ds) deallocate(longk_zp) deallocate(ex_x) deallocate(ex_y) deallocate(ex) deallocate(ex_buffer) end if ! record end time track_time = tock(track_clock) print *, 'Run time = ', track_time contains subroutine tick(t) integer, intent(out) :: t call system_clock(t) end subroutine tick ! returns time in seconds from now to time described by t real function tock(t) integer, intent(in) :: t integer :: now, clock_rate call system_clock(now,clock_rate) tock = real(now - t)/real(clock_rate) end function tock end program osc_track