program captured_momenta use precision_def use sim_utils use muon_mod use muon_interface implicit none type(kick_and_pulse_struct), allocatable :: kick_and_pulse(:) integer, parameter :: npoints = 10000 integer ix,i,k,j integer jx integer lun, ios integer lun1 integer kicks_per_section, tot_individual_kicks real(rp) Bpeak, xinf,xpinf, Aperture,momentum/3.095/, index/0.11/, R0/7.112/, length_kick/3.81/, Aperture_first_pass real(rp) eta, beta, xtune, k_min, B_to_kick, Bmin, Bmax, Bmin_first_pass real(rp) k_max real(rp) B real(rp) kick_to_gauss, pulse_height_normalization real(rp), allocatable :: delta_central(:) real(rp) pulse_height_sum, delta_sum, delta2, delta_ave, variance, sigma real(rp), allocatable :: slice_ave(:), slice_mean_square(:), Nxp(:) real(rp) sigma_d, gauss_diff, dp_min, dp_max, delta_p real(rp) Ngauss real(rp) root2, pulse_delay, kick_delay, delay real(rp) kick_tot, delta_time real(rp) kick_tot_cos, kick_tot_sin, kick_tot_sin_to_angle_beta, kick_tot1, kick_tot_cos1 real(rp) revolution/149.2/ real(rp) angle_correction real(rp) xp_limit2, xp_limit, xpinf2_avg, sigma_xp real(rp) y_min, y_max real(rp) start_time, end_time real(rp) dt, time_step_pulse, time_step_kick, pulse, pulse_sum, time real(rp) dxpinf ! change in the effective injection angle due to phase kicker at twopi xtune/4 real(rp) bxpinf !xpinf + dxpinf real(rp) kicker_ends(6)/8.4315,9.7015,10.292,11.562,12.152,13.422/ real(rp) kicker_start_ends(3)/8.4315,10.292,12.152/ real(rp) kicker_length/1.27/ real(rp), allocatable :: applied_kick(:), sub_kick(:) real(rp) deltak real(rp) kick_delay0 real(rp) b_2,b_3 ! quad sextupole and octupole moments in units of field index real(rp), allocatable :: dp_max_turn(:),dp_min_turn(:) real(rp) sub_time integer kick_points, pulse_points, sum_multiple_turn_kick/0/ integer n_momenta !number of momenta to write_amp_phase logical gaussian_momentum_distribution/.false./ logical first/.true./ logical itexists logical, allocatable :: mask(:) logical distributed_kick logical only_once/.true./, all_done/.false./ logical collimation_always/.true./ character*290 pulse_trace, kick_trace namelist/captured_momenta_input/pulse_trace, kick_trace, pulse_delay, kick_delay,xinf,xpinf, Aperture, Bpeak,index,sigma_d, sigma_xp, sum_multiple_turn_kick ,kicks_per_section, n_momenta, b_2, b_3 namelist/kick_delay_input/kick_delay sum_multiple_turn_kick=0 lun=lunget() open(unit=lun, file='captured_momenta_input.dat', STATUS='old', IOSTAT=ios) READ(lun,NML=captured_momenta_input, IOSTAT=ios) write(6, nml=captured_momenta_input) close(lun) kick_delay0=kick_delay inquire (file='kick_delay_input.dat', exist=itexists) if(itexists)then lun=lunget() open(unit=lun, file='kick_delay_input.dat', STATUS='old', IOSTAT=ios) READ(lun,NML=kick_delay_input, IOSTAT=ios) write(6, nml=kick_delay_input) close(lun) endif beta = R0/sqrt(1.-index) eta = R0/(1.-index) xtune = sqrt(1.-index) k_min = (xinf - Aperture+abs(xpinf*beta))/beta k_min = (xinf - sqrt(Aperture**2-(xpinf*beta)**2))/beta k_max = (xinf + sqrt(Aperture**2-(xpinf*beta)**2))/beta if(sigma_xp >0.)k_min= (xinf-Aperture)/beta B_to_kick = length_kick *0.3/momentum/1.e4*beta ! kick[G] X (B_to_kick) = kick[rad]*beta Bmin = k_min*beta/B_to_kick Aperture_first_pass = 0.050 Bmin_first_pass = (xinf - sqrt(Aperture_first_pass**2-(xpinf*beta)**2))/B_to_kick if(sigma_xp >0)Bmin_first_pass = (xinf - Aperture_first_pass)/B_to_kick ! Bmax = (xinf + Aperture-abs(xpinf*beta))/B_to_kick Bmax=k_max*beta/B_to_kick if(sigma_xp > 0)Bmax = (xinf + Aperture)/B_to_kick print '(a,es12.4)',' beta [m] = ', beta print '(a,es12.4)',' eta [m] = ', eta print '(a,es12.4)',' xtune = ', xtune print '(a,es12.4)',' k_min[rad] = ', k_min print '(a,es12.4,a,es12.4)',' Bmin [G] = ', Bmin,' Bmax [G] = ', Bmax print '(a,es12.4)',' Bpeak [G] = ', Bpeak print '(a,es12.4)',' Bmin_first_pass = ', Bmin_first_pass print '(a,es12.4)',' perfect angle = ', -xinf*cos(twopi*xtune/4.)/sin(twopi*xtune/4.)/beta write(6, '(a,es12.4)')' xinf= ', xinf*1000 write(6, '(a,es12.4)')' xpinf= ', xpinf*1000 write(6, '(a,es12.4)')' sigma_d= ', sigma_d write(6, '(a,es12.4)')' sigma_xp= ', sigma_xp lun1=lunget() open(lun1, file='fort.16') write(lun1, '(a,es12.4)') ' beta = ', beta write(lun1, '(a,es12.4)')' eta = ', eta write(lun1, '(a,es12.4)')' xtune = ', xtune write(lun1, '(a,es12.4)')' k_min = ', k_min write(lun1, '(a,es12.4)')' Bmin = ', Bmin write(lun1, '(a,es12.4)')' Bmin_first_pass = ', Bmin_first_pass write(lun1, '(a,es12.4)')' Bmax= ', Bmax write(lun1, '(a,es12.4)')' Bpeak= ',Bpeak write(lun1, '(a,es12.4)')' xinf= ', xinf*1000 write(lun1, '(a,es12.4)')' xpinf= ', xpinf*1000 write(lun1, '(a,es12.4)')' sigma_d= ', sigma_d write(lun1, '(a,es12.4)')' sigma_xp= ', sigma_xp if(kicks_per_section>0)then ! place kicks kicker_length/(kicks_per_section+1) of the length of the kicker from ends allocate(applied_kick(1:3*kicks_per_section), sub_kick(1:3*kicks_per_section)) k=0 do i = 1,3 !starting end do k= 1,kicks_per_section applied_kick(k+(i-1)*kicks_per_section) = kicker_start_ends(i) + k * kicker_length/(kicks_per_section+1) print '(i10,2es12.4)',k+(i-1)*kicks_per_section , kicker_start_ends(i),applied_kick(k+(i-1)*kicks_per_section) end do end do tot_individual_kicks = 3 * kicks_per_section print '(a,i10)',' tot_individual_kicks =',tot_individual_kicks endif if(sigma_d > 0) gaussian_momentum_distribution = .true. if(gaussian_momentum_distribution)then print '(a,es12.4)',' sigma_d = ',sigma_d Ngauss = 1./sigma_d/sqrt(twopi) ! Ngauss * exp(-((x-x0)/(2*sigma))**2) normalization if -infinity < x < infinity ! if limits are finite then define Ngauss so that int_DLow^DHigh (Ngauss *exp(-(delta/sigma)**2/2) d delta)=1 ! => 1./Ngauss = int _Dlow^Dhigh exp(-(delta/sigma)**2/2) d delta = erf(dp_max/root2/sigma) - erf(dp_min/root2/sigma) root2 = sqrt(2.) ! Ngauss = 1./(erf(dp_max/root2/sigma_d) - erf(dp_min/root2/sigma_d)) ! = -sigma_p**2 * Ngauss *(exp(-(dp_max/sigma_d)**2/2)-exp(-(dp_min/sigma_d)**2/2) endif ! B to central value of accepted momentum offset delta ! delta_central[%] = (xinf - B * B_to_kick)/eta/2. * 100 allocate(kick_and_pulse(0:npoints), delta_central(1:npoints)) allocate(slice_ave(1:npoints), slice_mean_square(1:npoints), Nxp(1:npoints)) allocate(mask(0:npoints)) allocate(dp_max_turn(0:sum_multiple_turn_kick),dp_min_turn(0:sum_multiple_turn_kick)) call combine_kick_and_pulse(kick_and_pulse, pulse_trace, kick_trace, pulse_delay, kick_delay, pulse_points,kick_points) ! kick_and_pulse(i)%kick_time - time ! kick_and_pulse(i)%kick - kick - arbitrary units ! kick_and_pulse(i)%pulse_height - intensity ! scale kick to Gauss ix = maxloc(kick_and_pulse(1:npoints)%kick,1) kick_to_gauss = Bpeak/kick_and_pulse(ix)%kick ! print '(a,5es12.4)','ix-2:ix+2',kick_and_pulse(ix-2)%kick,kick_and_pulse(ix-1)%kick,kick_and_pulse(ix)%kick,kick_and_pulse(ix+1)%kick,kick_and_pulse(ix+2)%kick do i=1,kick_points write(20,'(2es12.4)')kick_and_pulse(i)%kick_time,kick_to_gauss * kick_and_pulse(i)%kick end do ! normalize pulse height ix = maxloc(kick_and_pulse(1:npoints)%pulse_height,1) pulse_sum=0 print '(3a12)','pulse_sum','pulse_height','Delta t' do i=2,pulse_points !integrate pulse_height(t) dt pulse_sum = pulse_sum + kick_and_pulse(i)%pulse_height * (kick_and_pulse(i)%pulse_time - kick_and_pulse(i-1)%pulse_time) ! print '(3es12.4)',pulse_sum, kick_and_pulse(i)%pulse_height, (kick_and_pulse(i)%pulse_time - kick_and_pulse(i-1)%pulse_time) end do print '(a, es12.4)',' pulse sum =', pulse_sum ! n=0 pulse_height_sum =0 delta_sum = 0 delta2=0 start_time=0. print *, 'kick_points = ', kick_points kick_and_pulse(0)%pulse_time = kick_and_pulse(1)%pulse_time do i=1,kick_points write(20,'(2es12.4)')kick_and_pulse(i)%kick_time, kick_and_pulse(i)%kick * kick_to_gauss pulse=0 ! Get pulse height at the time of this slice time = kick_and_pulse(i)%kick_time call get_pulse_height_at_time(kick_and_pulse, time, pulse,ix) if(ix == 1)cycle do k=1,tot_individual_kicks sub_time = time + (applied_kick(k) - twopi*R0/4.)/c_light*1.e9 !ns call get_kick_at_time(kick_and_pulse,i, sub_time,sub_kick(k)) end do ! Add up all the kick at this time modula revolution period kick_tot = kick_and_pulse(i)%kick * sin(twopi*xtune/4.) ! first pass kick_tot_cos= kick_and_pulse(i)%kick * cos(twopi*xtune/4.) ! first pass kick_tot_sin= kick_and_pulse(i)%kick * sin(twopi*xtune/4.) ! first pass if(kicks_per_section > 0)then kick_tot = sum(sub_kick(1:tot_individual_kicks) * sin(xtune*applied_kick(1:tot_individual_kicks)/R0))/tot_individual_kicks ! first pass kick_tot_cos= sum(sub_kick(1:tot_individual_kicks) * cos(xtune*applied_kick(1:tot_individual_kicks)/R0))/tot_individual_kicks ! first pass kick_tot_sin= sum(sub_kick(1:tot_individual_kicks) * sin(xtune*applied_kick(1:tot_individual_kicks)/R0)) /tot_individual_kicks ! first pass kick_tot1 = kick_and_pulse(i)%kick/tot_individual_kicks * sum(sin(xtune*applied_kick(1:tot_individual_kicks)/R0)) ! first pass kick_tot_cos1= kick_and_pulse(i)%kick/tot_individual_kicks * sum(cos(xtune*applied_kick(1:tot_individual_kicks)/R0)) ! first pass ! kick_tot_sin1= kick_and_pulse(i)%kick/tot_individual_kicks * sum(sin(xtune*applied_kick(1:tot_individual_kicks)/R0)) ! first pass write(21,'(5es12.4)')time, kick_tot, kick_tot1, kick_tot_cos, kick_tot_cos1 endif ! max and min magnetic field to capture dxpinf=kick_tot_cos * kick_to_gauss * B_to_kick/beta ! call get_max_min_b(xinf, Aperture, dxpinf, xpinf, beta, Bmin, Bmax) bxpinf=xpinf+dxpinf k_min = (xinf - sqrt(Aperture**2-(bxpinf*beta)**2))/beta !added 3/11/22 k_min = (xinf - Aperture+abs(xpinf*beta))/beta k_max = (xinf + sqrt(Aperture**2-(bxpinf*beta)**2))/beta ! Bmin=k_min*beta/B_to_kick ! Bmax=k_max*beta/B_to_kick! 3/11/22 call get_momentum_range(kick_tot, kick_to_gauss, B_to_kick, eta, beta, kick_tot_cos, kick_tot_sin, xinf, xpinf, & Aperture, Bmin, Bmax, sigma_xp,delta_p, dp_max_turn(0), dp_min_turn(0), angle_correction, & delta_central(i), Nxp(i),B,kick_tot_sin_to_angle_beta, bxpinf, xpinf2_avg, xp_limit) if(angle_correction >=1.)cycle dp_max = dp_max_turn(0) dp_min=dp_min_turn(0) if(sum_multiple_turn_kick>0 .and. dp_min < dp_max)then ! if(sum_multiple_turn_kick >0 .and. pulse > 0)then jx = minloc(abs(kick_and_pulse(:)%kick_time -kick_and_pulse(i)%kick_time),1) mask(:) = .true. mask(1:jx)=.false. do k = 1,sum_multiple_turn_kick j = minloc(abs(kick_and_pulse(:)%kick_time -kick_and_pulse(i)%kick_time -k*revolution),1, mask) delta_time = kick_and_pulse(j)%kick_time - kick_and_pulse(i)%kick_time ! print '(4i10,4es12.4)',i,k,j,jx,kick_and_pulse(j)%kick_time, delta_time, kick_and_pulse(j)%kick_time -kick_and_pulse(i)%kick_time-k*revolution, & ! kick_and_pulse(j+1)%kick_time -kick_and_pulse(i)%kick_time-k*revolution if(k*revolution - 10 < delta_time .and. delta_time < k*revolution + 10)then if(kicks_per_section>0)then deltak = kick_and_pulse(j)%kick/tot_individual_kicks * sum(sin(twopi*xtune*(k+applied_kick(1:tot_individual_kicks)/twopi/R0))) kick_tot = kick_tot + kick_and_pulse(j)%kick/tot_individual_kicks * sum(sin(twopi*xtune*(k+applied_kick(1:tot_individual_kicks)/twopi/R0))) kick_tot_cos = kick_tot_cos + kick_and_pulse(j)%kick/tot_individual_kicks * sum(cos(twopi*xtune*(k+applied_kick(1:tot_individual_kicks)/twopi/R0))) kick_tot_sin = kick_tot_sin + kick_and_pulse(j)%kick/tot_individual_kicks * sum(sin(twopi*xtune*(k+applied_kick(1:tot_individual_kicks)/twopi/R0))) else deltak = kick_and_pulse(j)%kick * sin(twopi*xtune*(k+1/4.)) kick_tot = kick_tot + kick_and_pulse(j)%kick * sin(twopi*xtune*(k+1/4.)) kick_tot_cos = kick_tot_cos + kick_and_pulse(j)%kick * cos(twopi*xtune*(k+1/4.)) kick_tot_sin = kick_tot_sin + kick_and_pulse(j)%kick * sin(twopi*xtune*(k+1/4.)) endif write(18, '(a,2i5,8es12.4)')' k,i= ',k,i, kick_and_pulse(i)%kick_time, delta_time, kick_tot, & kick_tot * kick_to_gauss,kick_and_pulse(j)%kick * sin(twopi*xtune*(k+1/4.)) , deltak*kick_to_gauss, kick_tot_cos, kick_tot_sin endif mask(j)=.false. call get_momentum_range(kick_tot, kick_to_gauss, B_to_kick, eta, beta, kick_tot_cos, kick_tot_sin, xinf, xpinf, & Aperture, Bmin, Bmax, sigma_xp,delta_p, dp_max_turn(k), dp_min_turn(k), angle_correction, & delta_central(i), Nxp(i),B,kick_tot_sin_to_angle_beta, bxpinf, xpinf2_avg, xp_limit) if(collimation_always)then !collimation after each pass through kicker dp_min = max(dp_min,dp_min_turn(k)) dp_max = min(dp_max,dp_max_turn(k)) else !collimation only after all passes through kicker dp_min = dp_min_turn(k) dp_max = dp_max_turn(k) endif end do write(19, '(a,i5,4es12.4)')' i= ',i, kick_and_pulse(i)%kick_time, delta_time, kick_tot, kick_tot * kick_to_gauss endif ! B=kick_and_pulse(i)%kick* kick_to_gauss if(B>Bmin .and. B < Bmax)then if(start_time == 0.)start_time = kick_and_pulse(i)%kick_time end_time = kick_and_pulse(i)%kick_time print '(i10,5es12.4)',i,dxpinf,bxpinf, dp_max, dp_min, B if(dp_min == 0. .and. dp_max == 0.)then print '(a,2es12.4)','dp_min,dp_max',dp_min,dp_max cycle endif if(gaussian_momentum_distribution)then !compute average and rms momenta for each slice Ngauss = 2./(erf(dp_max/root2/sigma_d) - erf(dp_min/root2/sigma_d)) gauss_diff = exp(-(dp_max/sigma_d)**2/2.)-exp(-(dp_min/sigma_d)**2/2.) slice_ave(i) = -Ngauss*sigma_d/sqrt(twopi) * gauss_diff slice_mean_square(i) = -Ngauss*sigma_d/sqrt(twopi)*(dp_max*exp(-(dp_max/sigma_d)**2/2.)-dp_min*exp(-(dp_min/sigma_d)**2/2.)) + sigma_d**2 print '(i10,a,2es12.4)',i,'slice_ave,slice_mean_square', slice_ave(i),slice_mean_square(i) else slice_ave(i) =(dp_max+dp_min)/2. slice_mean_square(i) = 1/3.*(dp_max**3-dp_min**3)/(dp_max-dp_min) print '(a,i10,3es12.4)','slice_ave(i) ',i,slice_ave(i),xinf,B*B_to_kick endif time_step_kick = kick_and_pulse(i)%kick_time - kick_and_pulse(i-1)%kick_time if(gaussian_momentum_distribution)then pulse_height_sum = pulse_height_sum + pulse * time_step_kick /Nxp(i)/Ngauss !total stored in all slices delta_sum = delta_sum + slice_ave(i) * pulse * time_step_kick /Nxp(i)/Ngauss delta2 = delta2 + slice_mean_square(i) * pulse * time_step_kick / Nxp(i)/Ngauss print '(i10,a,5es12.4)',i,' pulse_height_sum,delta_sum,delta2,slice_ave(i),slice_mean_square(i)',pulse_height_sum,delta_sum,delta2,slice_ave(i), slice_mean_square(i) else pulse_height_sum = pulse_height_sum + pulse * time_step_kick delta_sum = delta_sum + slice_ave(i) * pulse * time_step_kick delta2 = delta2 + slice_mean_square(i) * pulse * time_step_kick endif endif if(first)then lun=lunget() open(lun,file='captured_momenta_out.dat') write(lun,'(14a14)')'kick_time', 'B[G]', 'mid range', ' pulse_height', 'slice_ave', 'slice_rms', 'dp_max','dp_min',' ',' xp_limit ',' Nxp(i) ',& 'xp_correction',' Ngauss' first=.false. endif if(pulse == 0.)cycle write(lun,'(13es14.4)')kick_and_pulse(i)%kick_time, B, delta_central(i), pulse, slice_ave(i), slice_mean_square(i), dp_max,dp_min, xpinf2_avg, xp_limit,& Nxp(i), angle_correction, Ngauss !!!!!!! if(B >= Bpeak*0.96 .and. only_once)then ! if(slice_ave(i-1) > 0 )then if(slice_ave(i) == 0)all_done=.true. print *,' call write_amp_phase_momentum at i= ',i print '(a,2es12.4)','slice_ave(i), (dp_max-dp_min)/2',slice_ave(i), (dp_max-dp_min)/2 call write_amp_phase_momentum(npoints,i,index, eta,beta, revolution, b_2, b_3, pulse, kick_and_pulse(i)%kick_time, & time_step_kick, slice_ave(i), (dp_max-dp_min)/2,xinf, bxpinf,kick_tot_sin_to_angle_beta , n_momenta, all_done) only_once=.false. all_done = .true. endif !!!!!! end do close(unit=lun) delta_ave = delta_sum / pulse_height_sum variance = delta2 / pulse_height_sum - delta_ave**2 sigma = sqrt(variance) print '(3(a,es12.4))',' Ngauss = ', Ngauss,' dp_max = ', dp_max,' dp_min = ', dp_min print '(a,es12.4,a,es12.4)',' = ', delta_ave, ' [m] = ', delta_ave*eta print '(a,es12.4, a, es12.4)',' delta2/pulse_height_sum = ', delta2/pulse_height_sum, ' sqrt(delta2/pulse_height_sum)*eta = ',sqrt(delta2/pulse_height_sum)*eta print '(a,es12.4,a,es12.4,a,es12.4)',' variance = ', variance, ' sigma_delta = ', sigma,' sigma_r [m] = ', sigma*eta print '(a,es12.4)',' pulse_height_sum = ', pulse_height_sum print '(a,es12.4)',' fraction_stored = ', pulse_height_sum/pulse_sum print '(a,i10)',' multi_kick_turns = ', sum_multiple_turn_kick write(6,'(a,es12.4)')'pulse length [ns] = ', end_time - start_time print '(a,i10)',' kicks_per_section = ', kicks_per_section write(lun1,'(a,es12.4)')'r_ave =', delta_ave*eta*1000 write(lun1,'(a,es12.4)')'sigma_r =', sigma*eta*1000 write(lun1,'(a,es12.4)')'pulse_length = ', end_time - start_time write(lun1,'(a,es12.4)')' pulse_height_sum = ', pulse_height_sum write(lun1,'(a,es12.4)')' fraction_stored = ', pulse_height_sum/pulse_sum write(lun1,'(a,i10)')' multi_kick_turns = ', sum_multiple_turn_kick write(lun1,'(a,i10)')' kicks_per_section = ',kicks_per_section close(lun1) inquire (file='capture_vs_timing.dat', exist=itexists) lun=lunget() if(.not. itexists)then open(unit=lun, file ='capture_vs_timing.dat') write(lun,'(9a12,a12,a12)')'pulse_delay','kick_delay','xinf','xpinf','Aperture','Bpeak','index','sigma_d','sigma_xp','multi-turn','fraction' close(lun) endif if(kick_delay /= kick_delay0)then print '(a)',':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::' print '(a)',' kick delay from kick_delay_input.dat does not agree with kick_delay from captured_momenta_input.dat' print '(a)',':::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::' endif open(unit=lun, file='capture_vs_timing.dat',status="old", access="append", action="write") write(lun,'(9es12.4,6x,i1,5x,es12.4)')pulse_delay, kick_delay,xinf,xpinf, Aperture, Bpeak,index,sigma_d, sigma_xp, sum_multiple_turn_kick, pulse_height_sum/pulse_sum close(lun) end program subroutine combine_kick_and_pulse(kick_and_pulse, pulse_trace, kick_trace, pulse_delay, kick_delay, pulse_points, kick_points) use precision_def use sim_utils use muon_mod use muon_interface, dummy => combine_kick_and_pulse implicit none type (kick_and_pulse_struct), allocatable :: kick_and_pulse(:) character*290 new_string, pulse_trace, kick_trace logical itexists logical end_flag/.false./ integer, parameter :: npoints=10000 integer reason integer lun,i integer nargs integer j integer kick_points, pulse_points character*20 delay_word character*1 junk real(rp) pulse_delay, kick_delay, delay real(rp) kick_time(0:npoints), b1(0:npoints),b2(0:npoints), b3(0:npoints) real(rp), allocatable:: pulse_time(:), pulse_height(:) real(rp) pulse_height_int ! nargs = command_argument_count() ! if (nargs >= 1)then ! call get_command_argument(1, kick_trace) ! print '(a,a)','kick trace from file: ',kick_trace ! elseif(nargs >=2)then ! call get_command_argument(2, pulse_trace) ! print '(a,a)','pulse trace from file: ',pulse_trace ! elseif(nargs ==3)then ! call get_command_argument(3,delay_word) ! read(delay_word,*)delay ! endif print '(a,a)','kick trace from file: ',kick_trace print '(a,a)','pulse trace from file: ',pulse_trace print '(a,es12.4)',' pulse delay [ns] =', pulse_delay print '(a,es12.4)',' kick delay [ns] =', kick_delay lun=lunget() open(unit=lun, file=kick_trace,status='old') i=0 do while(.true.) read(lun, '(a)', IOSTAT =reason)new_string ! print '(a)', new_string if(reason < 0) exit if(index(new_string,'#')/=0)cycle i=i+1 if(i>npoints)then print '(a,i10)', 'Number of points in kick_trace > upper bound ',npoints exit endif read(new_string, *, IOSTAT=reason) kick_time(i), b1(i), b2(i), b3(i) kick_time(i) = kick_time(i)-kick_delay kick_points = i kick_and_pulse(i)%kick_time = kick_time(i) kick_and_pulse(i)%kick = b1(i) end do close(unit=lun) allocate(pulse_time(0:npoints), pulse_height(0:npoints)) lun=lunget() inquire (file=trim(pulse_trace), exist=itexists) if (.not.itexists) print *, trim(pulse_trace) ,' file not found' if (itexists) print *, trim(pulse_trace) ,' file found' open(unit=lun, file=trim(pulse_trace),status='old') i=0 reason=0 do while(.true.) read(lun, '(a)', IOSTAT =reason)new_string ! print '(a)',trim(new_string) ! if(new_string(1:2) == ' ')cycle if(index(new_string,'#')/=0)cycle if(index(new_string,'NaN')/=0)cycle if(reason < 0) exit i=i+1 if(i>npoints)then print '(a,i10)', 'Number of points in pulse_trace > upper bound ',npoints goto 99 endif read(new_string, *) pulse_time(i), pulse_height(i), junk pulse_time(i) = pulse_time(i)-pulse_delay kick_and_pulse(i)%pulse_time = pulse_time(i) kick_and_pulse(i)%pulse_height = pulse_height(i) ! print '(i10,2es12.4,a)',i,pulse_time(i),pulse_height(i),junk pulse_points = i end do 99 close(unit=lun) do i=1,kick_points write(11,'(i10,2es12.4)')i, kick_and_pulse(i)%kick_time, kick_and_pulse(i)%kick end do do i=1,pulse_points write(12,'(i10,2es12.4)')i, kick_and_pulse(i)%pulse_time, kick_and_pulse(i)%pulse_height end do !do i=1,kick_points !if(i > npoints)then ! print '(a,i10)', 'i is greater than the number of points ', npoints ! goto 199 !endif !call kick_to_pulse_time(kick_time(i), pulse_time, delay, j) !pulse_height_int = (pulse_height(j+1)*(kick_time(i)-(pulse_time(j)+pulse_delay)) + pulse_height(j)*(pulse_time(j+1)+pulse_delay-kick_time(i)))/(pulse_time(j+1)-pulse_time(j)) !pulse_height_int = (pulse_height(j+1)*(kick_time(i)-(pulse_time(j))) + pulse_height(j)*(pulse_time(j+1)-kick_time(i)))/(pulse_time(j+1)-pulse_time(j)) !write(13,'(i10,3es12.4,i10,2es12.4, 2es12.4)')i,kick_time(i),b1(i), pulse_height_int, j, pulse_height(j), pulse_height(j+1), pulse_time(j), pulse_time(j+1) !write(14,'(i10,3es12.4)')i,kick_time(i),b1(i), pulse_height_int !end do !199 continue end subroutine subroutine kick_to_pulse_time (kick_time, pulse_time, delay, j) use precision_def implicit none real(rp) kick_time, delay real(rp), allocatable :: pulse_time(:) integer j,i !for kick time kick_time(i) find corresponding pulse_time(j) !print '(a,i10)','size(pulse_time)= ', size(pulse_time) do i=1,size(pulse_time)-2 ! print '(i10, 2es12.4)',i,kick_time, pulse_time(i)+delay ! if(pulse_time(i)+delay == kick_time)then if(pulse_time(i) == kick_time)then j=i exit endif ! if(pulse_time(i)+delay < kick_time .and. pulse_time(i+1)+delay > kick_time)then if(pulse_time(i) < kick_time .and. pulse_time(i+1) > kick_time)then j=i exit endif ! if(pulse_time(i)+delay > kick_time)then if(pulse_time(i) > kick_time)then j=0 exit endif end do end subroutine subroutine write_amp_phase_momentum(npoints,ix_slice,index, eta, beta, revolution, b_2, b_3, pulse, kick_time, time_step_kick, slice_ave, width, xinf, bxpinf, kick_tot_sin_to_angle_beta, n_momenta, all_done) use precision_def use sim_utils implicit none real(rp) index, eta, slice_ave, width, xinf, bxpinf, kick_tot_sin_to_angle_beta, beta, pulse,time_step_kick, kick_time real(rp) revolution, xp0 real(rp), save :: omega_c real(rp), allocatable, save :: omega_b(:), amp(:), delta(:), x0(:) real(rp) b_2, b_3 ! sextupole and quadrupole terms in units of the field index real(rp) Qx, Qx_0 real(rp) momentum_step real(rp), allocatable, save :: t(:), x(:,:) real(rp), allocatable, save :: xsum(:) real(rp) pulse_sum/0./ integer n_momenta, i,j, ix_slice, npoints,n_slice/0/ integer lun, lun1 integer, parameter :: ntimes = 4000 logical all_done, first_time/.true./ character*100, save :: format_string1,format_string2,format_string3 print *,'n_momenta=',n_momenta if(first_time)then allocate(t(1:ntimes)) allocate(xsum(1:ntimes)) xsum(1:ntimes)=0 allocate(x(1:ntimes,1:npoints)) allocate(omega_b(0:n_momenta), amp(0:n_momenta), delta(0:n_momenta), x0(0:n_momenta)) write(format_string1,'(a11,i2,a12)')'(a12,i2,a5,',n_momenta+1,'(es12.4,a1))' write(format_string2,'(a11,i2,a12)')'(a9,i2,a5,',n_momenta+1,'(es12.4,a1))' write(format_string3,'(a11,i2,a12)')'(a14,i2,a5,',n_momenta+1,'(es12.4,a1))' print *, 'format_string1 = ', format_string1 print *, 'format_string2 = ', format_string2 print *, 'format_string3 = ', format_string3 ! lun = lunget() ! open(unit=lun, file= 'momenta_amp_phase.dat') omega_c = twopi/(revolution/1000.) !revolution is time in units of ns first_time = .false. else ! lun = lunget() ! open(unit=lun, file= 'momenta_amp_phase.dat', status="old",access="append",action="write") endif pulse_sum = pulse_sum + pulse if(.not. all_done)then n_slice = n_slice + 1 xp0 = bxpinf momentum_step = 2*width/n_momenta print '(a,4es12.4)',' momentum_step,width,omega_c, kick_tot_sin_to_angle_beta ', momentum_step, width, omega_c,kick_tot_sin_to_angle_beta do i=0,n_momenta delta(i) = slice_ave-width + i * momentum_step ! print '(a,3es12.4)','xinf, delta(i),kick_tot_sin_to_angle_beta', xinf,delta(i),kick_tot_sin_to_angle_beta x0(i) = xinf - eta* delta(i) - kick_tot_sin_to_angle_beta amp(i) = sqrt((beta*bxpinf)**2 + x0(i)**2) Qx_0 = sqrt(1- (index+2*b_2*(eta*delta(i))**2+3*b_3*(eta*delta(i))**3)) Qx = (3*b_3/8/Qx_0 - 5*b_2**2/12/Qx_0**3)*amp(i)**2 +(1-Qx_0**2)/2/Qx_0 * delta(i) + Qx_0 omega_b(i) = Qx*omega_c/(1+delta(i)) ! print '(a,7es12.4)','delta,x0,xp0,amp,Qx_0,Qx,omega_b',delta(i),x0(i),xp0, amp(i), Qx_0, Qx, omega_b(i) ! write(lun, '(6es12.4)')pulse,time_step_kick,delta(i), x0(i), xp0, omega_b(i) end do ! write(lun, '(a1,/)')'#' ! write(lun,'(a17,i2)')'n_momenta = ', n_momenta ! write(lun,'(a17,es12.4)')'pulse = ', pulse ! write(lun,'(a17,es12.4)')'time_step_kick = ', time_step_kick ! write(lun,'(a17,es12.4)')'xp0 = ',xp0 ! write (lun ,format_string1)'array delta[',n_momenta+1,'] = [', (delta(i),',' , i=0,n_momenta-1),delta(n_momenta),']' ! write (lun ,format_string2)'array x0[',n_momenta+1,'] = [', (x0(i),',' , i=0,n_momenta-1),x0(n_momenta),']' ! write (lun ,format_string3)'array omega_b[',n_momenta+1,'] = [', (omega_b(i),',' , i=0,n_momenta-1),omega_b(n_momenta),']' ! write (lun ,format_string3)'array amp[',n_momenta+1,'] = [', (amp(i),',' , i=0,n_momenta-1),amp(n_momenta),']' ! close(lun) x(1:ntimes,1:npoints)=0 do j = 1, ntimes !micro s t(j)=float(j)*0.1 ! t(j) in units of 100 ns. kick_time is in ns. ! do i=0, n_momenta i= 10 x(j,n_slice) = x(j,n_slice) + xp0*sin(omega_b(i)*(t(j)+kick_Time/100.)) + x0(i)*cos(omega_b(i)*(t(j)+kick_time/100.)) + delta(i) * eta ! print *,'i,n_slice,x(j,n_slice)',i,n_slice,x(j,n_slice) ! end do ! x(j)=x(j)*pulse/(n_momenta+1) x(j,n_slice)=x(j,n_slice)/(n_momenta+1) * pulse * time_step_kick xsum(j) = xsum(j)+ x(j,n_slice) ! print *,'j,x,xsum',j,x(j,n_slice),xsum(j) end do write(21, '(a1,/)')'#' write(21,'(2es12.4)')(t(j),x(n_slice,i), j=1,ntimes) else ! (all_done) open(unit=lun1, file='cbo_vs_time.dat') write(lun1,'(2es12.4)')(t(j),xsum(j)/pulse_sum, j=1,ntimes) endif return end subroutine write_amp_phase_momentum !function delta_for_kick(xinf,B_to_kick,eta, Bfield) !use precision_def !implicit none !real(rp) xinf, B_to_kick, eta ! xinf-Bfield * B_to_kick)/eta*100/2 subroutine get_momentum_range(kick_tot, kick_to_gauss, B_to_kick, eta,beta, kick_tot_cos, kick_tot_sin, xinf, xpinf, & Aperture,Bmin, Bmax, sigma_xp,delta_p, dp_max, dp_min, angle_correction, delta_central, Nxp,B,kick_tot_sin_to_angle_beta,bxpinf, xpinf2_avg, xp_limit) use precision_def use sim_utils implicit none real(rp) kick_tot, Aperture, kick_to_gauss, B_to_kick, beta, kick_tot_cos, kick_tot_sin, xinf, xpinf, Bmin, Bmax, sigma_xp real(rp) delta_p, dp_max, dp_min, angle_correction real(rp) Nxp, delta_central real(rp) B, y_min, y_max, dxpinf, kick_tot_sin_to_angle_beta real(rp) bxpinf, eta, xpinf2_avg, xp_limit, xp_limit2 real(rp) root2 real(rp) ang_cor_plus, ang_cor_minus real(rp) k_min, k_max root2=sqrt(2.) B=kick_tot * kick_to_gauss dxpinf = kick_tot_cos * kick_to_gauss * B_to_kick/beta kick_tot_sin_to_angle_beta = kick_tot_sin * kick_to_gauss * B_to_kick delta_central=0 dp_max=0. dp_min=0. bxpinf=xpinf+dxpinf xpinf2_avg=bxpinf**2 angle_correction=0. ! Nxp(i)=1 Nxp=1 if(B>Bmin .and. B < Bmax)then ! if(start_time == 0.)start_time = kick_and_pulse(i)%kick_time ! end_time = kick_and_pulse(i)%kick_time ! n=n+1 xp_limit2 = (Aperture**2-(xinf-B*B_to_kick)**2) if(sigma_xp > 0)then !compute limits on angular acceptance and the average angular variance for each slice xp_limit = sqrt(Aperture**2-(xinf-B*B_to_kick)**2)/beta y_min = -xp_limit - bxpinf y_max = xp_limit - bxpinf Nxp= 2./(erf(y_max/root2/sigma_xp) - erf(y_min/root2/sigma_xp)) xpinf2_avg = -Nxp*sigma_xp/sqrt(twopi) * (y_max*exp(-(y_max/sigma_xp)**2/2.) - y_min*exp(-(y_min/sigma_xp)**2/2.)) + sigma_xp**2 & -2*Nxp*bxpinf*sigma_xp/sqrt(twopi)*(exp(-(y_max/sigma_xp)**2/2.)-exp(-(y_min/sigma_xp)**2/2.)) +bxpinf**2 print '(a,4es12.4)','xp_limit, y_min,y_max,sqrt(xpinf2_avg) ',xp_limit, y_min, y_max, sqrt(xpinf2_avg) endif angle_correction=xpinf2_avg*(beta)**2/xp_limit2 ! if(xpinf2_avg == 0)angle_correction=0 if(angle_correction >=1)return delta_central = (xinf - B*B_to_kick)/eta/2. * (1- angle_correction) ![%] delta_p = delta_central ! slice_ave(i) = delta_central(i) ! slice_mean_square(i) = (xinf/eta/2.)**2/12. dp_max = delta_p + Aperture/eta/2.*(1- angle_correction) dp_min = delta_p - Aperture/eta/2.*(1-angle_correction) ang_cor_plus=xpinf2_avg*beta**2/(Aperture-(xinf-B*B_to_kick)) ang_cor_minus=xpinf2_avg*beta**2/(Aperture+(xinf-B*B_to_kick)) dp_max=0.5*(Aperture+(xinf-B*B_to_kick)- xpinf2_avg*beta**2/(Aperture-(xinf-B*B_to_kick))) dp_min=-0.5*(Aperture-(xinf-B*B_to_kick)- xpinf2_avg*beta**2/(Aperture+(xinf-B*B_to_kick))) dp_max = dp_max/eta dp_min = dp_min/eta endif return end subroutine get_momentum_range subroutine get_pulse_height_at_time(kick_and_pulse, time, pulse,ix) use precision_def use sim_utils use muon_mod use muon_interface, dummy => get_pulse_height_at_time implicit none type (kick_and_pulse_struct), allocatable :: kick_and_pulse(:) real(rp) time, pulse real(rp) dt, time_step_pulse integer ix ix = minloc( abs (kick_and_pulse(:)%pulse_time - time),1) if(ix==1)return ix=ix-1 dt =kick_and_pulse(ix)%pulse_time -time if(abs(dt)<10.)then if(kick_and_pulse(ix-1)%pulse_time <= time .and. time < kick_and_pulse(ix)%pulse_time)then time_step_pulse = kick_and_pulse(ix)%pulse_time -kick_and_pulse(ix-1)%pulse_time pulse = kick_and_pulse(ix)%pulse_height - (kick_and_pulse(ix)%pulse_height - kick_and_pulse(ix-1)%pulse_height)*dt/time_step_pulse ! print '(a20,2i10,6es12.4)','if first',i,ix, kick_and_pulse(i)%kick_time,kick_and_pulse(ix-2)%pulse_time ,kick_and_pulse(ix-1)%pulse_time, & ! kick_and_pulse(ix)%pulse_time, kick_and_pulse(ix+1)%pulse_time,kick_and_pulse(ix+2)%pulse_time elseif(kick_and_pulse(ix)%pulse_time <= time .and. time < kick_and_pulse(ix+1)%pulse_time)then time_step_pulse = kick_and_pulse(ix+1)%pulse_time -kick_and_pulse(ix)%pulse_time pulse = kick_and_pulse(ix)%pulse_height - (kick_and_pulse(ix+1)%pulse_height - kick_and_pulse(ix)%pulse_height)*dt/time_step_pulse ! print '(a20,2i10,6es12.4)','if second',i,ix, kick_and_pulse(i)%kick_time,kick_and_pulse(ix-2)%pulse_time ,kick_and_pulse(ix-1)%pulse_time, & ! kick_and_pulse(ix)%pulse_time, kick_and_pulse(ix+1)%pulse_time,kick_and_pulse(ix+2)%pulse_time else print '(a20,i10,6es12.4)','check it out',ix, time,kick_and_pulse(ix-2)%pulse_time ,kick_and_pulse(ix-1)%pulse_time, & kick_and_pulse(ix)%pulse_time, kick_and_pulse(ix+1)%pulse_time,kick_and_pulse(ix+2)%pulse_time print *,'stop' ! stop endif write(17, '(i10,6es12.4)') ix,dt,time_step_pulse, kick_and_pulse(ix+1)%pulse_height, kick_and_pulse(ix)%pulse_height, pulse, time endif ! pulse height at slice time is pulse return end subroutine subroutine get_kick_at_time(kick_and_pulse,ix, time,kick) use precision_def use sim_utils use muon_mod use muon_interface, dummy => get_kick_at_time implicit none type (kick_and_pulse_struct), allocatable :: kick_and_pulse(:) real(rp) time, pulse real(rp) dt, time_step_pulse, kick, dgrid integer ix ! central index integer jx, jxx ! first guess is kick_and_pulse(ix)%kick kick=0 if(kick_and_pulse(ix-1)%kick == 0 .and. kick_and_pulse(ix)%kick==0 .and. kick_and_pulse(ix+1)%kick == 0)return jx=minloc(abs(kick_and_pulse(:)%kick_time-time),1) dt = kick_and_pulse(jx)%kick_time - time if(kick_and_pulse(jx)%kick_time-time >0)then do jxx= jx-1,1, -1 if(kick_and_pulse(jxx)%kick_time-time >0 )cycle dgrid=kick_and_pulse(jx)%kick_time - kick_and_pulse(jxx)%kick_time kick = (dt*kick_and_pulse(jxx)%kick +(dgrid-dt) * kick_and_pulse(jx)%kick)/dgrid exit end do endif if(kick_and_pulse(jx)%kick_time-time < 0)then do jxx = jx+1,size(kick_and_pulse(:)%kick_time)-1 if(kick_and_pulse(jxx)%kick_time-time <0)cycle dgrid=kick_and_pulse(jxx)%kick_time - kick_and_pulse(jx)%kick_time kick = ((dgrid +dt) *kick_and_pulse(jx)%kick + dt * kick_and_pulse(jxx)%kick)/dgrid exit end do endif return end subroutine