! get a list of directories and combine content of Time_Dep_Moments_at_END.dat into single file program compile_Time_dep use precision_def use nr use bmad use sim_utils use compile_interface !use muon_mod implicit none type g2moment_struct real(rp) sigma(6,6), ave(6), third(6) real(rp) max_sigma(6,6), min_sigma(6,6), max_ave(6), min_ave(6) end type TYPE (g2moment_struct), allocatable :: moments(:) !moment structure for each time bin, 1000 turns => an array with 1000*149e-9/bin_width TYPE (g2moment_struct), allocatable :: moments_sum(:) !moment structure for each time bin, 1000 turns => an array with 1000*149e-9/bin_width type running_moment_struct real(rp) product(1:6,1:6), sum(1:6), third(6) integer n end type type (running_moment_struct), allocatable :: running_moment(:) type (running_moment_struct), allocatable :: running_moments_sum(:) !moment structure for each time bin, 1000 turns => an array with 1000*149e-9/bin_width integer, parameter:: n_time_bins=700000 integer tbin, tbin_max/0/, lun integer i, n, lines/0/ integer isign/1/ integer m integer k, ix integer j integer nbins_max/590000/ integer reason integer lun1 integer nargs integer n_files/0/ integer lun_0 integer m_rev integer peak_bin real(rp) two/2./, zero/0./ complex(rp), allocatable ::data_fft_n(:), data_fft_x(:), data_fft_y(:), data_fft_xx(:), data_fft_yy(:), data_fft_xxx(:), data_fft_yyy(:) real(rp) time_bin/1.e-8/ real(rp), allocatable :: time(:) real(rp) freq real(rp) total_hits real(rp) time_off real(rp) qy,qx, betax, eta real(rp) y_last, yp_last, betay,Cp, x_last real(rp) x_average,pz,Ce real(rp) s real(rp) f_rev, Q_vert character*290 string, new_string character*100 dir, dir_file character*12000 file_string/' '/ character*20 time_off_char character*200 cwd logical itexists logical end_flag/.false./, efield_end_flag/.false./ allocate(moments(1:nbins_max), running_moment(1:nbins_max)) allocate(running_moments_sum(1:nbins_max), moments_sum(1:nbins_max)) allocate(time(1:nbins_max)) nargs = cesr_iargc() if (nargs == 1)then call cesr_getarg(1, time_off_char) print *,time_off_char read(time_off_char,*)time_off print *, 'Using ', ' time_off = ', time_off else time_off = 35.e-6 print '(a,es12.4)',' time_off = ', time_off endif running_moments_sum(1:nbins_max)%n = 0 forall(i=1:6)moments_sum(1:nbins_max)%ave(i)=0 call getcwd(cwd) string='ls > out.dat' call execute_command_line(string) lun=lunget() open(unit=lun,file='out.dat', status='old') do while(.true.) read(lun, '(a)', IOSTAT = reason)new_string if(reason < 0)exit if( (index(new_string,'2018')==0) .and. (index(new_string,'2019')==0) .and. len(trim(new_string)) /= 7)cycle dir=trim(new_string) dir_file = trim(dir)//'/co_tunes_scraping_off.dat' print '(a)',dir_file inquire (file=trim(dir_file), exist=itexists) if (.not.itexists) cycle call get_tunes(dir,Qx,Qy,betax,betay,eta) exit end do close(lun) open(unit=lun,file='out.dat', status='old') do while(.true.) read(lun, '(a)', IOSTAT = reason)new_string if(reason < 0)exit if( (index(new_string,'2018')==0) .and. (index(new_string,'2019')==0) .and. len(trim(new_string)) /= 7)cycle dir=trim(new_string) dir_file = trim(dir)//'/Time_Dep_Moments_at_END.dat' print '(a)',dir_file inquire (file=trim(dir_file), exist=itexists) if (.not.itexists) cycle n_files = n_files+1 file_string = trim(file_string)//' '//trim(trim(dir)//'/EndOfTracking_phase_space.dat') call get_data_from_EndOfTracking_phase_space(dir, time_off, qx, qy, betay) lines = 0 lun1=lunget() open(unit=lun1, file = trim(dir_file) ) do while(.true.) read(lun1,'(a)', end=99)string if(index(string,'<')/= 0)cycle lines = lines + 1 if(lines > nbins_max)then print *,' Moment array has size smaller than number of bins in data file' stop endif ! read(string,*,end=99)tbin, time, number, orb%vec(1:6) ! print '(a)',string read(string,'(i10)')i if(i > nbins_max)cycle read(string,'(i10,es18.8,i10,21(es12.4))',end=99)i, time(i),running_moment(i)%n, moments(i)%ave(:),& (moments(i)%sigma(1+j,1+j),moments(i)%sigma(1+j,2+j),moments(i)%sigma(2+j,2+j), j=0,4,2), & moments(i)%third(1:6) moments_sum(i)%ave(:) = (moments_sum(i)%ave(:) * running_moments_sum(i)%n + moments(i)%ave(:) * running_moment(i)%n) & /(running_moments_sum(i)%n + running_moment(i)%n) running_moments_sum(i)%n = running_moments_sum(i)%n+running_moment(i)%n ! if((lines/10000)*10000 == lines)print *, lines tbin_max = max(tbin_max, i) end do 99 continue close(lun1) end do !cycle through directories close(lun) if(n_files == 0)then lun_0=lunget() open(unit=lun_0, access='sequential', file = 'NoData.dat',position='append') write(lun_0,'(2a)')' No data in directory ',trim(cwd) close(unit=lun_0) goto 199 endif ! fourier transform of Number vs time n = log(float(tbin_max))/log(two) m = two**n allocate(data_fft_n(1:m), data_fft_x(1:m), data_fft_y(1:m), data_fft_xx(1:m), data_fft_yy(1:m), data_fft_xxx(1:m), data_fft_yyy(1:m)) do i=1,m data_fft_n(i) = cmplx (running_moment(i)%n, zero) data_fft_x(i) = cmplx(moments(i)%ave(1),zero) * 2*sin((twopi/2*i)/m) * sin((twopi/2.*i)/m) data_fft_y(i) = cmplx(moments(i)%ave(3),0.) * 2*sin((twopi/2*i)/m) * sin((twopi/2.*i)/m) end do call four1(data_fft_n(1:m),isign) call four1(data_fft_x(1:m),isign) call four1(data_fft_y(1:m),isign) print *,' tbin_max = ', tbin_max,' log(two) = ', log(two), ' n= ',n,' m = ',m,' size(data_fft_n) = ', size(data_fft_n) f_rev = 1./149.2e-3 !revolution frequency MHz m_rev = f_rev *1.e6*m*time_bin/2 peak_bin=maxloc(abs(data_fft_y(1:m_rev)),1) Q_vert = float(peak_bin)/m_rev/2. print '(a,es12.4,a,i10,a,i10)', 'f_rev = ',f_rev,' m_rev = ',m_rev,' peak_bin = ', peak_bin print *,' Q_vert = ', Q_vert lun = lunget() open(unit = lun, file = 'EndTime.dat') print *,' EndTime.dat' write(lun,'(a)')'#! Data from files ' write(lun,'(a)')'#!',dir_file write(lun, '(a9,1x,a12,3(3a14))')'time_bin','freq','Total Counts','Real fft','Imag fft',& '','Real fft','Imag fft',& '','Real fft','Imag fft' do i = 1, m freq=float(i)/time_bin/m/1.e6 write(lun,'(i10,es12.4,(i14,2es14.6),6(3es14.6))' ) i,freq,running_moments_sum(i)%n, data_fft_n(i),& moments_sum(i)%ave(1), data_fft_x(i), & moments_sum(i)%ave(3), data_fft_y(i) end do total_hits = sum(running_moments_sum(:)%n) print '(a,es12.4)','Total hits = ', total_hits end_flag = .true. efield_end_flag = .true. call vertical_offset_vs_s(y_last, s, end_flag) call pitch_vs_offset(y_last,yp_last,betay,Cp,end_flag) call efield_vs_radial_offset(x_average,pz,Ce,efield_end_flag) call get_positron_tunes(x_last,y_last,time,end_flag) print '(a)','cat '//trim(file_string)//' > all_EndOfTracking_phase_space.dat' call execute_command_line('cat '//trim(file_string)//' > all_EndOfTracking_phase_space.dat') 199 continue end program subroutine get_data_from_EndOfTracking_phase_space(dir, time_off, qx, qy, betay) use bmad implicit none character*100 dir character*50 dir_file character*600 long_string character*2 time_off_char logical itexists logical first/.true./ logical end_flag/.false./, efield_end_flag/.false./ real(rp) x,beta, Cp, Ce, Cp_2, Ce_2, Cp_rms, Ce_rms, CCe, CCp, CCe_correct real (rp), save::Cp_avg(1:300), Ce_avg(1:300), CCe_avg(1:300), CCp_avg(1:300), CCe_correct_avg(1:300), y_avg(1:300), x_avg(1:300) real(rp), save :: fe,fp,p0c,mass real(rp) pz, time, BetaXE, BetaXE_time, BetaDotB real(rp) Cp_avg_all, sig_p, Ce_avg_all, sig_e, CCe_avg_all, CCp_avg_all real(rp) time_off real(rp), save :: y2(1:300), mean_radial2(1:300) real(rp) y2_avg_all, y2_rms, mean_radial2_all, y_last, x_average, yp_last, y_amp2, y_average, x_last real(rp) qy, qx, betay real(rp), save :: quad_index,r0 real(rp) Ce_standard, Cp_standard real(rp) s integer ix,i,n integer, save:: lun,m integer, save :: lun1 ! beta2(x) = ((x+1)*p0c)**2/((x+1)**2*p0c**2+mass**2) ! pz=column(9) - pz ! time=column(10) - time ! BetaXE=column(17) - (beta X E)_y ! BetaXE-timecolumn(19) - (Beta X E) time ! BetaDotB=column(21) - (beta dot B)_y ! if(first)then if(int(time_off*1.e6) < 10)then write(time_off_char,'(i1)')int(time_off*1.e6) else write(time_off_char,'(i2)')int(time_off*1.e6) endif lun=lunget() open(lun,file='summary_efield_pitch_'//trim(time_off_char)//'us.dat') lun1=lunget() quad_index = 1.-qx**2 open(lun1,file='bad_runs.dat') write(lun,'(a,es12.4,a,es12.4,a,es12.4)')'time_off = ', time_off,' Qx = ',qx, ' quad_index = ', quad_index write( lun,'(a10,1x,a16,1x,a10,1x,20a16)')'file', 'dir', 'n muons', 'Cp_avg', 'Cp_rms', 'Ce_avg', 'Ce_rms', 'Cp_avg_all', 'sig_p', 'Ce_avg_all', & 'sig_e','y2_avg_all','y2_rms','mean_radial2_all','Ce_standard','Cp_standard','CCe_avg_all','CCp_avg_all','CCe_ced_avg_all','Qx','Qy','','' write( lun1,'(a10,1x,a16,1x,a10,1x,20a16)')'file', 'dir', 'n muons', 'Cp_avg', 'Cp_rms', 'Ce_avg', 'Ce_rms', 'Cp_avg_all', 'sig_p', 'Ce_avg_all', & 'sig_e','y2_avg_all','y2_rms','mean_radial2_all','Ce_standard','Cp_standard','CCe_avg_all','CCp_avg_all','CCe_ced_avg_all','Qx','Qy','','' first=.false. fe = 2/c_light/1.4513007*1.e6 fp = 1/1.4513007*1.e6/2. p0c = 3.094353005E9 mass=105.6583715E6 r0=7.112 m=0 endif dir_file = trim(dir)//'/EndOfTracking_phase_space.dat' inquire (file=dir_file, exist = itexists) if(.not. itexists)then print *,dir_file,' does not exist' return endif m=m+1 !files print '(2a)', 'open ',dir_file open(unit=5,file=dir_file) n=0 Cp_avg(m)=0 Cp_2=0 Ce_avg(m)=0 Ce_2 = 0 y2(m)=0 mean_radial2(m)=0 CCe_avg(m)=0 CCe_correct_avg(m)=0 CCp_avg(m)=0 x_avg(m)=0 y_avg(m)=0 do while(.true.) read(5,'(a)',end=99)long_string if(index(long_string,'muon')/=0)cycle i=0 ix=0 do while (i<26) call string_trim(long_string(ix+1:), long_string, ix) i=i+1 read(long_string(1:ix),*)x if(i == 4)then x_last = x elseif(i == 6)then y_last = x elseif(i == 7)then yp_last=x elseif(i == 9)then pz=x elseif (i == 10)then time=x elseif (i == 11)then s=x elseif(i == 17)then BetaXE=x elseif(i == 19)then BetaXE_time=x elseif(i == 21)then BetaDotB=x elseif(i == 24)then x_average = x elseif(i == 26)then y_average = x endif end do if(time > time_off .and. BetaXE_time > 149.e-9)then n=n+1 beta = ((pz+1)*p0c)**2/((pz+1)**2*p0c**2+mass**2) Cp = BetaDotB/beta/BetaXE_time * fp Ce = BetaXE*pz /BetaXE_time * fe Cp_avg(m) = (Cp_avg(m)*(n-1)+Cp)/float(n) Cp_2 = (Cp_2*(n-1) + Cp*Cp)/float(n) Ce_avg(m) = (Ce_avg(m)*(n-1)+Ce)/float(n) Ce_2 = (Ce_2*(n-1)+Ce*Ce)/float(n) y2(m) = (y2(m)*(n-1) + y_last**2)/float(n) mean_radial2(m) = (mean_radial2(m)*(n-1) + x_average**2)/float(n) y_amp2 = y_last**2+yp_last**2* betay**2 call convolute(x_average,y_last, y_amp2,CCe, CCp,quad_index, CCe_correct ) CCe_avg(m) = (CCe_avg(m)*(n-1)+CCe)/float(n) CCe_correct_avg(m) = (CCe_correct_avg(m)*(n-1)+CCe_correct)/float(n) CCp_avg(m) = (CCp_avg(m)*(n-1)+CCp)/float(n) x_avg(m)= (x_avg(m)*(n-1)+x_average)/float(n) y_avg(m) = (y_avg(m)*(n-1)+y_average)/float(n) call vertical_offset_vs_s(y_last,s,end_flag) call pitch_vs_offset(y_last,yp_last,betay,Cp,end_flag) call efield_vs_radial_offset(x_average,pz,Ce,efield_end_flag) call get_positron_tunes(x_last, y_last, time, end_flag) endif end do 99 continue close(unit=5) if(n == 0)return if(Ce_avg(m) <2. .and. Ce_avg(m) > 0.1)then Cp_rms = sqrt(Cp_2-Cp_avg(m)**2) Ce_rms = sqrt(Ce_2-Ce_avg(m)**2) Cp_avg_all = sum(Cp_avg(1:m))/float(m) x = sum(Cp_avg(1:m)**2)/float(m)-Cp_avg_all sig_p = sqrt(sum(Cp_avg(1:m)**2)/float(m)-Cp_avg_all**2) Ce_avg_all = sum(Ce_avg(1:m))/float(m) sig_e = sqrt(sum(Ce_avg(1:m)**2)/float(m)-Ce_avg_all**2) mean_radial2_all = sum(mean_radial2(1:m))/float(m) y2_avg_all = sum(y2(1:m))/float(m) y2_rms = sqrt(sum(y2(1:m))/float(m)-y2_avg_all**2) Ce_standard = 2*beta*quad_index*(1-quad_index)*mean_radial2_all/r0**2 * 1.e6 Cp_standard = y2_rms**2 * quad_index/2./r0**2 * 1.e6 CCe_avg_all = sum(CCe_avg(1:m)/float(m)) CCp_avg_all = sum(CCp_avg(1:m)/float(m)) print '(i10,1x,a,i10,1x,15es16.6)',m,dir, n, Cp_avg(m), Cp_rms, Ce_avg(m), Ce_rms, Cp_avg_all, sig_p, Ce_avg_all, sig_e, y2_avg_all, & y2_rms, mean_radial2_all, Ce_standard, Cp_standard, CCe_avg_all, CCp_avg_all write( lun,'(i10,1x,a16,1x,i10,1x,29es16.6)')m,dir, n, Cp_avg(m), Cp_rms, Ce_avg(m), Ce_rms , Cp_avg_all, sig_p, Ce_avg_all, sig_e, y2_avg_all, & y2_rms, mean_radial2_all, Ce_standard, Cp_standard, CCe_avg_all, CCp_avg_all,sum(CCe_correct_avg(1:m))/float(m), Qx, Qy, sum(y_avg(1:m))/float(m), sum(x_avg(1:m))/float(m) else write( lun1,'(i10,1x,a16,1x,i10,1x,20es16.6)')m,dir, n, Cp_avg(m), Cp_rms, Ce_avg(m), Ce_rms , Cp_avg_all, sig_p, Ce_avg_all, sig_e, y2_avg_all, & y2_rms, mean_radial2_all, Ce_standard, Cp_standard, CCe_avg_all, CCp_avg_all,sum(CCe_correct_avg(1:m))/float(m), Qx, Qy,sum(y_avg(1:m))/float(m), sum(x_avg(1:m))/float(m) m=m-1 endif return end subroutine get_data_from_EndOfTracking_phase_space subroutine convolute(x_average,y_last, y_amp2, Ce, Cp, quad_index, Ce_correct) use bmad implicit none real(rp) x_average, Ce,Cp, x, y_last, y_amp2, quad_index, Ce_correct logical first/.true./, itexists integer lun, i, ix, lines1, lines2 character*11 var_name(0:6)/'a0','a1','a2','a3','a4','wave_number','amp'/ ! character*11 pvar_name(0:4)/'b0','b1','b2','b3','b4'/ character*11 pvar_name(0:4)/'b0','b2','b4','b6','b8'/ character*100 string character*50 fit_parameters_file, pitch_fit_parameters_file real(rp), save :: a(0:6), b(0:4), quad_index_0, n1n if(first)then lun=lunget() fit_parameters_file = '../integral_fit_parameters.dat' fit_parameters_file = '../integral_fit_parameters_ref.dat' inquire (file=trim(fit_parameters_file), exist=itexists) if(.not. itexists) then print *, 'fit parameter file',trim(fit_parameters_file), ' does not exist' stop endif lines1 = 0 quad_index_0 =3.3701E-01**2 !from /ref n1n = quad_index_0*(1-quad_index_0) open(unit=lun,file=trim(fit_parameters_file)) do while(.true.) read(lun,'(a)',end=99)string lines1 = lines1 + 1 do i=0,6 ix =index(string,trim(var_name(i))//' =') if(ix /= 0)then call string_trim(string(ix:), string, ix) read(string(ix+4:),*)a(i) print '(a,es12.4)',var_name(i),a(i) endif end do end do 99 close(lun) lun=lunget() pitch_fit_parameters_file = '../spin_pitch_fit_parameters.dat' pitch_fit_parameters_file = '../spin_pitch_fit_parameters_ref.dat' inquire (file=trim(pitch_fit_parameters_file), exist=itexists) if(.not. itexists) then print *, 'fit parameter file', pitch_fit_parameters_file, ' does not exist' stop endif lines2=0 open(unit=lun,file=pitch_fit_parameters_file) do while(.true.) read(lun,'(a)',end=199)string lines2 = lines2 + 1 do i=0,4 ix =index(string,trim(pvar_name(i))//' =') if(ix /= 0)then call string_trim(string(ix:), string, ix) read(string(ix+4:),*)b(i) print '(a,es12.4)',pvar_name(i),b(i) endif end do end do 199 close(lun) if(lines1 == 0 .or. lines2 == 0)then print *,' fit parameter files are empty' stop endif first=.false. endif x = x_average*1000. ! Ce = a(0)+a(1)*x+a(2)*x**2 *(1+delta_a2*(quad_index-quad_index_0)) +a(3)*x**3 +a(4)*x**4+a(6)*(cos(a(5)*x)-1) Ce = a(0)+a(1)*x+a(2)*x**2 +a(3)*x**3 +a(4)*x**4+a(6)*(cos(a(5)*x)-1) Ce_correct=Ce*quad_index*(1-quad_index)/n1n ! x=sqrt(y_amp2) * 1000. ! Cp= -b(2)*x**2 +b(1)*x + b(3)*x**3+b(4)*x**4 x=y_last*1000 Cp = b(0)+b(1)*x**2+b(2)*x**4+b(3)*x**6+b(4)*x**8 return end subroutine convolute subroutine pitch_vs_offset(y, yp,betay,Cp, end_flag) use precision_def use sim_utils implicit none real(rp) y, yp, Cp real(rp) betay real(rp) amp, average_all real(rp), save :: yp2(-45:45), Cp_avg(-45:45), avg_cp(-45:45) logical first/.true./ logical end_flag integer ind integer, save:: lun integer, save :: m(-45:45) if(first)then lun=lunget() open(lun, file = 'pitch_vs_offset.dat') write(lun,'(a10,2a16,a10,a16)')'y_offset','','','number','' yp2(:)=0 Cp_avg(:) = 0 first = .false. m(:)=0 endif if(.not. end_flag)then ind = y*1000+sign(.5,y) if(-45 <= ind .and. ind <=45)then amp = (y**2/betay+yp**2*betay) m(ind)=m(ind)+1 yp2(ind) = (yp2(ind)*(m(ind)-1)+amp/betay)/float(m(ind)) Cp_avg(ind) = (Cp_avg(ind)*(m(ind)-1)+Cp)/float(m(ind)) endif ! yp2_all = yp2_all + amp/betay else do ind = 0,45 avg_cp(ind) = sum(Cp_avg(-ind:ind)*m(-ind:ind))/sum(m(-ind:ind)) end do average_all = sqrt(sum(yp2(-45:45)*m(-45:45))/sum(m(-45:45))) write(lun,'(a,es16.4,a,i10,a)')'# Average pitch^2 = ',average_all,' with ',sum(m(:)),' measurements' write(lun,'(a,es16.4,a,i10,a)')'# C_p = ', sum(Cp_avg(:)*m(:))/sum(m(:)),' with ',sum(m(:)),' measurements' write(lun,'(i10,2es16.8,i10, es16.8)')(ind,yp2(ind),Cp_avg(ind),m(ind), avg_cp(abs(ind)), ind=-45,45) endif return end subroutine pitch_vs_offset subroutine efield_vs_radial_offset(x_average,pz,Ce,efield_end_flag) use precision_def use sim_utils implicit none real(rp) x_average, pz, Ce real(rp) average_all real(rp), save :: Ce_avg(-45:45) logical first/.true./ logical efield_end_flag integer ind integer, save:: lun integer, save :: m(-45:45) if(first)then lun=lunget() open(lun, file = 'efield_vs_radial_offset.dat') write(lun,'(a10,a16,a10)')'x_e','','number' Ce_avg(:) = 0 first = .false. m(:)=0 endif if(.not. efield_end_flag)then ind = x_average*1000+sign(.5,x_average) if(-45 <= ind .and. ind <=45)then m(ind)=m(ind)+1 Ce_avg(ind) = (Ce_avg(ind)*(m(ind)-1)+Ce)/float(m(ind)) endif else write(lun,'(a,es12.4,1x,a,1x,i10)')'Average Ce = ', sum(Ce_avg(-45:45)*m(-45:45))/sum(m(-45:45)), ' Hits = ', sum(m(-45:45)) write(lun,'(i10,es16.8,i10)')(ind,Ce_avg(ind),m(ind), ind=-45,45) endif return end subroutine efield_vs_radial_offset subroutine get_tunes(dir, qx,qy,betax,betay,eta) use sim_utils implicit none logical itexists, twiss_found/.false./ character*290 string, new_string, out_string character*100 dir character*50 dir_file character*10 word(1:5)/'Qx =','Qy =','Betax =','Betay =','eta ='/ character*10 test_word real(rp) twiss(1:5), qx,qy,betax,betay,eta integer ix, reason, lun, lun1 integer iword, w_len if(twiss_found) return ix=-1 lun1=lunget() dir_file = trim(dir)//'/co_tunes_scraping_off.dat' inquire (file=dir_file, exist=itexists) if(.not.itexists) return open(lun1, file=trim(dir_file)) do while(ix <= 0) read(lun1,'(a)', IOSTAT=reason)new_string if(reason < 0) exit do iword = 1,5 w_len = len(trim(word(iword))) test_word = trim(word(iword)) ix = index(new_string,trim(word(iword))) if(ix /=0)then call string_trim(new_string(ix+len(trim(word(iword))):), out_string,ix) read(out_string(1:ix),*)twiss(iword) print *,word(iword),twiss(iword) endif qx=twiss(1) qy=twiss(2) betax=twiss(3) betay=twiss(4) eta=twiss(5) end do end do twiss_found = .true. close(lun1) return end subroutine get_tunes subroutine vertical_offset_vs_s(y_last,s, end_flag) use precision_def use sim_utils implicit none real(rp) y_last,s real(rp) s_turn real(rp) circumference/44.686/ real(rp), save :: y2_avg(1:45), y_avg(1:45), y2_rms(1:45) real(rp) angle logical first/.true./ logical end_flag integer ind integer i integer, save:: lun integer, save :: p(1:45) if(first)then lun=lunget() open(unit=lun, file = 'vertical_offset_vs_s.dat') write(lun,'(a1,a10,3a16)')'#','ring_index','','','rms_y' y2_avg(:)=0 y_avg(:)=0 p(:)=0 first=.false. endif if(.not. end_flag) then s_turn = mod(s, circumference) angle = s_turn/circumference * 360 ind = int(s_turn)+1 ! every 1 m p(ind) = p(ind) + 1 y2_avg(ind) = (y2_avg(ind)*(p(ind)-1) + y_last**2)/float(p(ind)) y_avg(ind) = (y_avg(ind)*(p(ind)-1)+y_last)/float(p(ind)) y2_rms(ind) = sqrt(y2_avg(ind)-y_avg(ind)**2) else do i = 1,45 write(lun,'(i10,1x,3es16.8)')i, y2_avg(i), y_avg(i), y2_rms(i) end do endif return end subroutine vertical_offset_vs_s subroutine get_positron_tunes(x,y,time, end_flag) use bmad use nr implicit none logical first/.true./, end_flag integer max_bins/10000/ integer, allocatable, save::n(:) integer ind, lun integer i integer isign/1/ integer nmax real(rp) x,y, time real(rp) time_bin/149.2e-9/ real(rp), allocatable,save :: y_avg(:), x_avg(:) complex(rp), allocatable, save :: data_fft_x(:), data_fft_y(:) real(rp) zero/0./ if(first)then allocate(n(1:max_bins), x_avg(1:max_bins), y_avg(1:max_bins)) allocate(data_fft_x(1:max_bins), data_fft_y(1:max_bins)) n(1:max_bins)=0 y_avg(1:max_bins)=0 x_avg(1:max_bins)=0 first=.false. endif if(end_flag)then nmax=2**int(log(float(max_bins))/log(2.d0)) do i=1,max_bins data_fft_x(i) = cmplx(x_avg(i), zero) data_fft_y(i) = cmplx(y_avg(i), zero) * 2*sin((twopi/2*i)/nmax) * sin((twopi/2.*i)/nmax) end do print *,' nmax = ', nmax call four1(data_fft_x(1:nmax),isign) call four1(data_fft_y(1:nmax),isign) lun=lunget() open(unit=lun, file = 'positron_pos_vs_time.dat') do i = 1, max_bins write(lun,'(i10,1x,i10,1x,6es16.8)')i,n(i), x_avg(i),y_avg(i), data_fft_x(i), data_fft_y(i) end do close(unit=lun) deallocate(data_fft_x, data_fft_y, n,x_avg, y_avg) else ind = time/time_bin+0.5 if(ind > max_bins)return n(ind) = n(ind)+1 y_avg(ind) = ((n(ind)-1)*y_avg(ind) + y)/n(ind) x_avg(ind) = ((n(ind)-1)*x_avg(ind) + x)/n(ind) endif return end subroutine get_positron_tunes