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 n,lun, ios 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 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 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(:) real(rp) deltak real(rp) kick_delay0 real(rp) b_2,b_3 ! quad sextupole and octupole moments in units of field index 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./ 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 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 write(16, '(a,es12.4)') ' beta = ', beta write(16, '(a,es12.4)')' eta = ', eta write(16, '(a,es12.4)')' xtune = ', xtune write(16, '(a,es12.4)')' k_min = ', k_min write(16, '(a,es12.4)')' Bmin = ', Bmin write(16, '(a,es12.4)')' Bmin_first_pass = ', Bmin_first_pass write(16, '(a,es12.4)')' Bmax= ', Bmax write(16, '(a,es12.4)')' Bpeak= ',Bpeak write(16, '(a,es12.4)')' xinf= ', xinf*1000 write(16, '(a,es12.4)')' xpinf= ', xpinf*1000 write(16, '(a,es12.4)')' sigma_d= ', sigma_d write(16, '(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)) 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)) 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 ix = minloc( abs (kick_and_pulse(:)%pulse_time - kick_and_pulse(i)%kick_time),1) if(ix==1)cycle ix=ix-1 dt =kick_and_pulse(ix)%pulse_time -kick_and_pulse(i)%kick_time if(abs(dt)<10.)then if(kick_and_pulse(ix-1)%pulse_time <= kick_and_pulse(i)%kick_time .and. kick_and_pulse(i)%kick_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 <= kick_and_pulse(i)%kick_time .and. kick_and_pulse(i)%kick_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,2i10,6es12.4)','check it out',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 print *,'stop' ! stop endif write(17, '(2i10,7es12.4)') i,ix,dt,time_step_pulse, kick_and_pulse(ix+1)%pulse_height, kick_and_pulse(ix)%pulse_height, pulse, kick_and_pulse(i)%kick_time, kick_and_pulse(i)%kick endif ! pulse height at slice time is pulse ! 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 = kick_and_pulse(i)%kick/tot_individual_kicks * sum(sin(xtune*applied_kick(1:tot_individual_kicks)/R0)) ! first pass kick_tot_cos= kick_and_pulse(i)%kick/tot_individual_kicks * sum(cos(xtune*applied_kick(1:tot_individual_kicks)/R0)) ! first pass kick_tot_sin= kick_and_pulse(i)%kick/tot_individual_kicks * sum(sin(xtune*applied_kick(1:tot_individual_kicks)/R0)) ! first pass endif if(sum_multiple_turn_kick>0 .and. kick_tot * kick_to_gauss > Bmin_first_pass)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,6es12.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 endif mask(j)=.false. 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 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(i)=0 slice_ave(i)=0. slice_mean_square(i)=0. dp_max=0. dp_min=0. bxpinf=xpinf+dxpinf xpinf2_avg=bxpinf**2 angle_correction=0. Nxp(i)=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(i)= 2./(erf(y_max/root2/sigma_xp) - erf(y_min/root2/sigma_xp)) xpinf2_avg = -Nxp(i)*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(i)*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 delta_central(i) = (xinf - B*B_to_kick)/eta/2. * (1- angle_correction) ![%] delta_p = delta_central(i) 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) ! print '(i10,7es12.4)',i,dxpinf,bxpinf, dp_max, dp_min, y_min, y_max, B 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 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 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,'(12a14)')'kick_time', 'B[G]', 'mid range', ' pulse_height', 'slice_ave', 'slice_rms', 'half_width',' ',' xp_limit ',' Nxp(i) ',& 'xp_correction',' Ngauss' first=.false. endif if(pulse == 0.)cycle write(lun,'(12es14.4)')kick_and_pulse(i)%kick_time, B, delta_central(i), pulse, slice_ave(i), slice_mean_square(i), (dp_max-dp_min)/2, xpinf2_avg, xp_limit,& Nxp(i), angle_correction, Ngauss !!!!!!! if(B >= Bpeak*0.96 .and. only_once)then ! if(slice_ave(i-1) > 0 .and. .not. all_done)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) ! if(all_done)only_once=.false. 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(16,'(a,es12.4)')'r_ave =', delta_ave*eta*1000 write(16,'(a,es12.4)')'sigma_r =', sigma*eta*1000 write(16,'(a,es12.4)')'pulse_length = ', end_time - start_time write(16,'(a,es12.4)')' pulse_height_sum = ', pulse_height_sum write(16,'(a,es12.4)')' fraction_stored = ', pulse_height_sum/pulse_sum write(16,'(a,i10)')' multi_kick_turns = ', sum_multiple_turn_kick write(16,'(a,i10)')' kicks_per_section = ',kicks_per_section 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), save :: pulse_sum, pulse_sumx integer n_momenta, i,j, ix_slice, npoints,n_slice/0/ integer lun, lun1 integer, parameter :: ntimes = 40000 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') pulse_sum=0. pulse_sumx=0. 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 * time_step_kick write(150,'(3es12.4)')kick_time,pulse, time_step_kick 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.01 ! t(j) in units of 100 ns. kick_time is in ns. do i=0, n_momenta x(j,n_slice) = x(j,n_slice) + xp0*sin(omega_b(i)*(t(j)+kick_Time/10.)) + x0(i)*cos(omega_b(i)*(t(j)+kick_time/10.)) + 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) ! xsum(j) = xsum(j)+ x(j,n_slice) ! print *,'j,x,xsum',j,x(j,n_slice),xsum(j) end do write(21, '(a1,/)')'#' if((n_slice/10)*10 == n_slice)then pulse_sumx = pulse_sumx +pulse*time_step_kick write(100+n_slice/10,'(i10,2es12.4)')(n_slice,t(j),x(j,n_slice)*pulse*time_step_kick, j=1,ntimes) x(1:ntimes,n_slice)=x(1:ntimes,n_slice) * pulse * time_step_kick xsum(1:ntimes)=xsum(1:ntimes)+x(1:ntimes,n_slice) write(130+n_slice/10,'(i10,2es12.4)')(n_slice,t(j),xsum(j)/pulse_sumx, j=1,ntimes) endif print '(a,es12.4)','pulse_sumx =',pulse_sumx else ! (all_done) xsum(1:ntimes) = xsum(1:ntimes)/pulse_sumx lun1=lunget() open(unit=lun1, file='cbo_vs_time.dat') write(lun1,'(2es12.4)')(t(j),xsum(j), 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