program spin_momentum_correlation 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 real(rp) revolution/149.2/ real(rp) angle_correction real(rp) xp_limit2, xp_limit, xpinf2_avg, sigma_xp real(rp) sigma_x 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) xmin/0.068/,xmax/0.086/, dx, pmin, pmax, px real(rp) x0/0.077/,dxp_dx real(rp) alpha_i,beta_i real(rp), allocatable :: dp_max_turn(:),dp_min_turn(:) integer kick_points, pulse_points, sum_multiple_turn_kick/0/ integer n_momenta !number of momenta to write_amp_phase integer p_steps, x_steps, i_dx 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./ logical Bpeak_only/.false./ character*290 pulse_trace, kick_trace namelist/spin_momentum_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, p_steps, x_steps, dxp_dx, sigma_x, Bpeak_only,alpha_i,beta_i sum_multiple_turn_kick=0 lun=lunget() open(unit=lun, file='spin_momentum_input.dat', STATUS='old', IOSTAT=ios) READ(lun,NML=spin_momentum_input, IOSTAT=ios) write(6, nml=spin_momentum_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 dxp_dx = -alpha_i/beta_i 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, '(2(a,es12.4))')' alpha_i=',alpha_i,'; beta_i=',beta_i write(6, '(a,es12.4)')' dxp_dx= ', dxp_dx 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(6, '(a,es12.4)')' sigma_x= ', sigma_x 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 write(lun1, '(a,es12.4)')' sigma_x= ', sigma_x write(lun1, '(a,es12.4)')' index= ', index write(lun1, '(2(a,es12.4))')' alpha_i=',alpha_i,'; beta_i=',beta_i write(lun1, '(a,es12.4)')' dxp_dx= ', dxp_dx if(Bpeak_only)write(lun1, '(a,a)')' Bpeak_only =', '"T"' if(.not. Bpeak_only)write(lun1, '(a,a)')' Bpeak_only =', '"F"' ! set up spacing of xinf dx = (xmax-xmin)/x_steps pmin = -Aperture/eta pmax = Aperture/eta px = (2*Aperture/eta)/p_steps ! 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)) 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 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! cycle through xinf do i_dx = 1, x_steps xinf = i_dx * dx + xmin xpinf = dxp_dx * (xinf - x0) ! 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 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) dp_max = dp_max_turn(0) dp_min=dp_min_turn(0) if(sum_multiple_turn_kick>0 .and. dp_min < dp_max)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 endif if(dp_max > dp_min )then ! print '(a,4es12.4)','debug2',dp_min, dp_max, B, Bpeak if((Bpeak_only .and. B >= Bpeak*0.95) .or. .not. Bpeak_only)then ! print '(i10,4es12.4)', (k,dp_min_turn(k), dp_max_turn(k), dp_min, dp_max,k=0,sum_multiple_turn_kick) ! print '(4es12.4)', xinf,dp_min, dp_max, dp_max-dp_min write(19, '(a,i5,8es12.4)')' i= ',i, kick_and_pulse(i)%kick_time, delta_time,pulse, kick_tot,xinf,xpinf, dp_min, dp_max call histogram_height_momentum_spin(Aperture, eta, pulse*time_step_pulse, sigma_x, dx, px, x_steps, p_steps, xmin, pmin, xinf, dp_min, dp_max, .false.) endif endif end do !cycle through xinf ! 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.)cycle 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 !!!!!!! end do !cycle through times call histogram_height_momentum_spin(Aperture, eta, pulse, sigma_x, dx, px, x_steps, p_steps, xmin, pmin, xinf, dp_min, dp_max, .true.) 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