! Find all phase space files ! -command line input: first two arguments are initial and final timestamped subdirectory. ! To get all dated subdirectories '00000000_000000' '99999999_999999' ! The next narg-2 argmuments are the phase_space files to be combined (anything with 'phase_space" in the name) ! -namelist input: if there are zero arguments on the command line, read data from 'energy_spin_vs_time_input.dat' (an example is appended to this file ! Get phase space of distribution prior to injection into ring of those particles that are captured, (or survive at least until some number of turns) ! In aide of understanding how correlations of phase space coordinates and phase are converted into momentum phase correlations of stored beam ! ! Then do the same for distributions at fixed times - like the distribution that is collected with the calorimteres and used to fit omega_a ! all_Energy_vs_time_spin_0.dat - phase space, spin, phi of surviving partilces from 'DistStartInjLine' and 'MARK_INFLECTOR_DS' ! all_Energy_spin_time.dat - phase space, spin, phi of surviving partilces from all distributions collected at specific times ! Fit omega_phi for each particle that survives vs time ! all_omega_a.dat ! MARK_INFLECTOR_DS_phase_space.dat_trimmed - trimmed version of phase space from a single subdirectory. phi_omega is adjusted modulo twopi so entire distribution is in the same cycle ! DistStartInjLine_phase_space.dat_trimmed - trimmed version of phase space from a single subdirectory. phi_omega is adjusted modulo twopi so entire distribution is in the same cycle ! program energy_spin_vs_time use nr use energy_vs_time_mod use energy_vs_time_interface implicit none type muon_at_struct type(muon_alive_struct), allocatable :: muon_save(:) end type muon_at_struct type(muon_at_struct), allocatable :: muon_at(:) type(muon_alive_struct)muon character*300 string, new_string, word character*300, allocatable :: phase_space_file(:) character*300 dir_file, dir_file_t, psft character*300 dir(300) character*20 dir_min, dir_max character*200 cwd character*24000 file_string/' '/, file_string_3/' '/, file_string_4/' '/, file_string_5/' '/, file_string_6/' '/, file_string_7/' '/ character*24000, allocatable :: time_string(:) character*30 format_string character*100, allocatable :: heading(:) character*100 trimmed_file real(rp) slope, offset real(rp) a,b,chi2, siga, sigb, q real(rp), allocatable :: phase(:), t(:), err(:) real(rp) phi_ave real(rp) vec5/0./,location/0./, jx/0./, jy/0./ real(rp), allocatable :: phase_sum(:), time_sum(:) real(rp) delta_ix_mu logical itexists, first/.true./, first3/.true./ logical write_phi_fixed/.true./, first_write/.true./ logical op logical all_there integer reason, i, j integer, allocatable :: n_left(:), ix_mu(:) integer lun, lun2, lun3, lun4, lun5,lun6/0/ integer number_files integer n integer k integer nevents/1000000/ integer min_date,max_date, min_time, max_time, date, time integer nargs integer ix integer ix_numb integer, allocatable:: muon_number(:) integer ind integer total integer ios integer, allocatable :: number(:) namelist/energy_spin_vs_time_input/nargs, dir_min, dir_max, phase_space_file dir(1)=' ' ! phase_space_file(1)='END000' ! phase_space_file(2)='END028' time_string = ' ' nargs = command_argument_count() if (nargs == 1)then call cesr_getarg(1,dir(1)) ! the data is all in a single directory input directory name on command line ! print *, 'len(trim(dir(1)))= ',len(trim(dir(1))) number_files=1 print '(i10,1x,a19)', number_files,dir(1) else if( nargs == 2 .or. nargs>=4)then call cesr_getarg(1,dir_min) call cesr_getarg(2,dir_max) if(nargs >= 4)then ix=nargs-2 do i=3,nargs call cesr_getarg(i,phase_space_file(i-2)) end do endif endif if(nargs == 0) then OPEN (UNIT=5, FILE='energy_spin_vs_time_input.dat', STATUS='old', IOSTAT=ios) READ (5, NML=energy_spin_vs_time_input, IOSTAT=ios) print *, 'ios=', ios rewind(unit=5) CLOSE(5) WRITE(6,NML=energy_spin_vs_time_input) ix=nargs-2 allocate(phase_space_file(1:ix)) OPEN (UNIT=5, FILE='energy_spin_vs_time_input.dat', STATUS='old', IOSTAT=ios) READ (5, NML=energy_spin_vs_time_input, IOSTAT=ios) print *, 'ios=', ios rewind(unit=5) CLOSE(5) WRITE(6,NML=energy_spin_vs_time_input) ix=nargs-2 endif if(nargs == 2 .or. nargs >=4)then read(dir_min(1:8),*)min_date read(dir_min(10:15),*)min_time read(dir_max(1:8),*)max_date read(dir_max(10:15),*)max_time endif allocate(muon_at(1:ix), time_string(1:ix),heading(1:ix),phase(1:ix)) allocate(t(1:ix),err(1:ix), phase_sum(1:ix),time_sum(1:ix)) allocate(ix_mu(1:ix), number(1:ix)) allocate(n_left(1:ix)) n_left(:)=0 phase_sum(1:ix)=0 time_sum(1:ix)=0 number(1:ix)=0 print '(a,1x,a)',trim(phase_space_file(1)), trim(phase_space_file(2)) call getcwd(cwd) !if all directories are to be combined, input nothing on command line string='ls > out.dat' call execute_command_line(string) lun=lunget() open(unit=lun,file='out.dat', status='old') n=0 do while(.true.) read(lun, '(a)', IOSTAT = reason)new_string if(reason < 0)exit if( (index(new_string,'2025')==0) .and. (index(new_string,'2019')==0) .and. (index(new_string,'2021')==0) .and. (index(new_string,'2022')== 0).and. (index(new_string,'2023')== 0).and. (index(new_string,'2024')== 0) )cycle n=n+1 dir(n)=trim(new_string) print '(i10,1x,a19)',n, dir(n) number_files = n end do close(lun) endif !_______________________________________________________________________ k=0 allocate(muon_number(1:nevents)) do n=1, number_files if(nargs == 2 .or. nargs >= 4)then read(dir(n)(1:8),*)date read(dir(n)(10:15),*)time if(datemax_date)cycle if(time>max_time)cycle endif string='ls '//trim(dir(n))//' > '//trim(dir(n))//'/out.dat' call execute_command_line(string) lun=lunget() inquire (file = trim(dir(n))//'/'//'out.dat', exist = itexists) if(.not. itexists)then print *, trim(dir(n))//'/'//'out.dat', 'does not exist' cycle endif call find_all_files(trim(dir(n))//'/'//'out.dat' , phase_space_file, ix, all_there) if(.not. all_there)cycle open(unit=lun,file=trim(dir(n))//'/'//'out.dat', status='old') lun2=lunget() open (unit=lun2,file = trim(dir(n))//'/'//'Energy_vs_time_0.dat') if(first)write(lun2,'(a10,1x,4a16,1x,a10)')'muon index', 'momentum','time',' initial disp',' inital angle','ix' first=.false. lun3=lunget() open (unit=lun3,file = trim(dir(n))//'/'//'Energy_vs_time_spin_0.dat') if(nargs >4)then lun4=lunget() open(lun4,file=trim(dir(n))//'/'//'Energy_spin_time.dat') lun5=lunget() open (lun5,file=trim(dir(n))//'/'//'omega_a.dat') endif if(first3)then if(nargs >=4)then write(lun3,'(20x,a,36x,80x,a)') trim(phase_space_file(1)), trim(phase_space_file(2)) write(lun3,'(a10,1x,20a16)')'muon index', 'x','px/p','y','py/p','dp/p' ,'time','spin_x','spin_y', & 'spin_z','acos(s.p)', 'x','px/p','y','py/p','dp/p' ,'time','spin_x','spin_y','spin_z','acos(s.p)' write(lun4,'(a1,a)') ('#',trim(phase_space_file(i)), i=1,ix) if(ix<10)write(format_string,'("(a1,i9,",i1,"(6i12))")')ix if(ix>=10)write(format_string,'("(a1,i9,",i2,"(6i12))")')ix write(lun4, format_string)'#',(j, j=1,6*ix+1) heading(1:ix)= 'dp/p time spin_x spin_y spin_z acos(s.p)' if(ix<10)write(format_string,'("(a10,",i1,"(6a12))")' )ix if(ix>=10)write(format_string,'("(a10,a12,",i2,"(6a12))")' )ix write(lun4, format_string)'muon index','omega_a (b)', (heading(j)(1:4), heading(j)(6:9), & heading(j)(11:16), heading(j)(18:23), heading(j)(25:30), heading(j)(32:40), j=1,ix) write(lun5,'(a10,a12,a16,17a12)')'muon_numb','offset','slope','chi2','phi','spin_x','spin_y','spin_z','phi','spin_x', & 'spin_y','spin_z','phi','spin_x','spin_y','spin_z','phi','spin_x','spin_y','spin_z' endif first3=.false. endif !first3 do while(.true.) read(lun, '(a)', IOSTAT = reason)new_string if(reason < 0)exit do i = 1,ix if(index(new_string,phase_space_file(i))/=0)then dir_file = trim(dir(n))//'/'//trim(new_string) call get_data_from_END_phase_space(dir_file, muon_at(i)%muon_save, n_left(i)) if(i>2)call save_and_concatenate(dir_file, phase_space_file(i)) phi_ave = sum(muon_at(i)%muon_save(1:n_left(i))%phi,1,muon_at(i)%muon_save(1:n_left(i))%spin(3)/=0) total=0 do j=1,n_left(i) if(muon_at(i)%muon_save(j)%spin(3)/=0 .and. muon_at(i)%muon_save(j)%spin(2)/=0) total = total +1 end do phi_ave = phi_ave/total first_write = .true. do j=1,n_left(i) ! if(muon_at(i)%muon_save(j)%spin(3) /= 0)write(112,'(es12.4)')muon_at(i)%muon_save(j)%phi if(muon_at(i)%muon_save(j)%phi - phi_ave > twopi/2)muon_at(i)%muon_save(j)%phi = muon_at(i)%muon_save(j)%phi-twopi if(muon_at(i)%muon_save(j)%phi - phi_ave < -twopi/2)muon_at(i)%muon_save(j)%phi = muon_at(i)%muon_save(j)%phi+twopi if(muon_at(i)%muon_save(j)%phi - phi_ave > twopi/2)muon_at(i)%muon_save(j)%phi = muon_at(i)%muon_save(j)%phi-twopi if(muon_at(i)%muon_save(j)%phi - phi_ave < -twopi/2)muon_at(i)%muon_save(j)%phi = muon_at(i)%muon_save(j)%phi+twopi if(write_phi_fixed)then if(first_write)then inquire(lun6,opened=op) if(op)close(lun6) lun6=lunget() trimmed_file = trim(trim(dir_file)//'_trimmed') open(lun6,file = trimmed_file) write(lun6,'(a10,15a12)')'muon','jx','jy','x','xp','y','yp','vec5','dp/p','time','location','spin_x','spin_y','spin_z','phi' first_write = .false. endif muon=muon_at(i)%muon_save(j) write(lun6,'(i10,15es12.4)')muon%muon_numb, jx, jy, muon%x, muon%xp, muon%y,muon%yp, vec5,muon%pz,muon%time,location,muon%spin, muon%phi endif end do if(write_phi_fixed)close(lun6) dir_file_t = trim(trim(dir_file)//'_trimmed') psft=trim(trim(phase_space_file(i))//'_trimmed') if(i>2)call save_and_concatenate(dir_file_t, psft) exit endif end do ! if(index(new_string,'spin_vs_momentum_at_time.dat')/=0)then dir_file = trim(dir(n))//'/'//trim(new_string) call combine_spin_vs_mom_at_time_data(dir_file) endif ! end do close(unit=lun) if(any(n_left(1:ix) <=0))then j=minloc(n_left(1:ix) ,1) print '(a,a,a1,a)',' no events in file ',dir(n),'/',trim(phase_space_file(j)) print *,trim(phase_space_file(j)) stop endif do j=1,ix-1 muon_at(j)%muon_save(:)%mask=.true. end do do i = 1, n_left(ix) ! write(109,*) ! write(109,'(a1,i10)')'#',i ! write(109,*) ix_numb = muon_at(ix)%muon_save(i)%muon_numb ix_mu(ix) = i do j=1,ix-1 ix_mu(j)=minloc(abs(muon_at(j)%muon_save(:)%muon_numb - ix_numb),1,muon_at(j)%muon_save(:)%mask) muon_at(j)%muon_save(ix_mu(j))%mask = .false. delta_ix_mu = muon_at(j)%muon_save(ix_mu(j))%muon_numb - ix_numb if(delta_ix_mu /= 0)then print *, 'muon numbers out of synch' stop endif end do if(k+1>nevents) then print *,' more than ',nevents,' events' cycle endif k=k+1 muon_number(k) = muon_at(1)%muon_save(i)%muon_numb ind = ix_mu(1) write(lun2,'(i10, 1x, 4es16.8,1x,i16)')muon_at(1)%muon_save(ind)%muon_numb, muon_at(1)%muon_save(ind)%pz, & muon_at(1)%muon_save(ind)%time, muon_at(1)%muon_save(ind)%x, muon_at(1)%muon_save(ind)%xp, muon_at(1)%muon_save(ind)%ix write(lun3,'(i10, 1x, 20es16.8)')ix_numb, (muon_at(j)%muon_save(ix_mu(j))%x, & muon_at(j)%muon_save(ix_mu(j))%xp,muon_at(j)%muon_save(ix_mu(j))%y,muon_at(j)%muon_save(ix_mu(j))%yp, muon_at(j)%muon_save(ix_mu(j))%pz, & muon_at(j)%muon_save(ix_mu(j))%time, muon_at(j)%muon_save(ix_mu(j))%spin, muon_at(j)%muon_save(ix_mu(j))%phi, j=1,2) do j=1,ix t(j) = muon_at(j)%muon_save(ix_mu(j))%time phase(j) = muon_at(j)%muon_save(ix_mu(j))%phi ! write(109,'(2es12.4)')t(j),phase(j) end do call fit(t(3:ix),phase(3:ix),a,b,siga,sigb,chi2,q) do j=1,ix if(chi2 >0 .and. chi2 < 10)then phase_sum(j) = phase_sum(j)+ phase(j) time_sum(j) = time_sum(j) + t(j) number(j) = number(j)+1 endif end do write(format_string,'("(i10,es12.4,",i1,"(6es12.4))")') ix write(lun4, format_string) ix_numb, b, & (muon_at(j)%muon_save(ix_mu(j))%pz,muon_at(j)%muon_save(ix_mu(j))%time, muon_at(j)%muon_save(ix_mu(j))%spin, muon_at(j)%muon_save(ix_mu(j))%phi , j=1,ix) write(lun5,'(i10,es12.4,es16.8,17es12.4)')ix_numb, a,b,chi2,phase(1), muon_at(1)%muon_save(ix_mu(1))%spin, phase(2), muon_at(2)%muon_save(ix_mu(2))%spin, & phase(3), muon_at(3)%muon_save(ix_mu(3))%spin, phase(4), muon_at(4)%muon_save(ix_mu(4))%spin end do n_left = 0 close(unit=lun2) close(unit=lun3) close(unit=lun4) close(unit=lun5) ! if(ix == 2)then ! file_string = trim(file_string)//' '//trim(dir(n))//'/'//'Energy_vs_time_0.dat' ! endif if(ix >=2)then file_string = trim(file_string)//' '//trim(dir(n))//'/'//'Energy_vs_time_0.dat' file_string_3 = trim(file_string_3)//' '//trim(dir(n))//'/'//'Energy_vs_time_spin_0.dat' file_string_4 = trim(file_string_4)//' '//trim(dir(n))//'/'//'Energy_spin_time.dat' file_string_5 = trim(file_string_5)//' '//trim(dir(n))//'/'//'omega_a.dat' file_string_6 = trim(file_string_6)//' '//trim(dir(n))//'/'//trim(phase_space_file(1))//'_trimmed' file_string_7 = trim(file_string_7)//' '//trim(dir(n))//'/'//trim(phase_space_file(2))//'_trimmed' ! do i=3,ix ! word = phase_space_file(i) ! time_string(i) = trim(time_string(i))//' '//trim(dir(n))//'/'//trim(word) ! print '(i10,1x,2a)',i,trim(word),trim(time_string(i)) ! end do endif if(any(dir(1:n-1) == dir(n)))then print *, ' directory ', n, ' same as ', ' earlier' endif print '(i10, 1x,a)', n, trim(dir(n)) end do call fit_line(muon_at(1)%muon_save, k, slope, offset) open(lun6,file= 'phase_vs_time.dat') write(lun6,'(3a12)') 'time','phase','number' do i=1,ix if(number(i) /= 0)then write(lun6,'(3es12.4)')time_sum(i)/number(i), phase_sum(i)/number(i) else print '(a)',' Poor fit for phase vs time for '//trim(phase_space_file(i)) endif end do ! if(ix <=2)then ! print '(a)','cat '//trim(file_string)//' > all_Energy_vs_time_0.dat' ! call execute_command_line('cat '//trim(file_string)//' > all_Energy_vs_time_0.dat') ! endif if(ix >= 2)then print '(a)','cat '//trim(file_string)//' > all_Energy_vs_time_0.dat' call execute_command_line('cat '//trim(file_string)//' > all_Energy_vs_time_0.dat') print '(a)','cat '//trim(file_string_3)//' > all_Energy_vs_time_spin_.dat' call execute_command_line('cat '//trim(file_string_3)//' > all_Energy_vs_time_spin.dat') print '(a)','cat '//trim(file_string_4)//' > all_Energy_spin_time.dat' call execute_command_line('cat '//trim(file_string_4)//' > all_Energy_spin_time.dat') print '(a)','cat '//trim(file_string_5)//' > all_omega_a.dat' call execute_command_line('cat '//trim(file_string_5)//' > all_omega_a.dat') print '(a)','cat '//trim(file_string_6)//' > all_'//trim(phase_space_file(1))//'_trimmed' call execute_command_line('cat '//trim(file_string_6)//' > all_'//trim(phase_space_file(1))//'_trimmed') print '(a)','cat '//trim(file_string_7)//' > all_'//trim(phase_space_file(2))//'_trimmed' call execute_command_line('cat '//trim(file_string_7)//' > all_'//trim(phase_space_file(2))//'_trimmed') ! do i = 3,ix ! print '(a)','cat '//trim(time_string(i))//' > all_'//trim(phase_space_file(i)) ! call execute_command_line('cat '//trim(time_string(i))//' > all_'//trim(phase_space_file(i))) ! end do endif end program energy_spin_vs_time subroutine get_data_from_END_phase_space(dir_file,muon_alive, n_left) use energy_vs_time_mod use energy_vs_time_interface, dummy_except => get_data_from_END_phase_space implicit none type(muon_alive_struct), allocatable :: muon_alive(:) character*300 dir_file character*600 long_string logical itexists real(rp) x, pz, time,disp, ang real(rp) disp_y, ang_y real(rp) spin(3) real(rp) acos_sp real(rp) sxp(3), pvec(3) real(rp) sdotp real(rp) svec(3) integer ix,i,n,muon_number, n_left if(.not. allocated(muon_alive))allocate(muon_alive(1:1000000)) inquire (file=dir_file, exist = itexists) if(.not. itexists)return print '(2a)', 'open ',dir_file open(unit=5,file=dir_file) n=0 do while(.true.) read(5,'(a)',end=99)long_string if(index(long_string,'muon')/=0 .or. index(long_string,'#')/=0)cycle i=0 ix=0 do while (i<15) call string_trim(long_string(ix+1:), long_string, ix) i=i+1 read(long_string(1:ix),*)x if(i == 1) then muon_number = int(x) n=n+1 elseif(i==4)then disp = x elseif(i==5)then ang = x elseif(i==6)then disp_y = x elseif(i==7)then ang_y = x elseif( i == 10)then time = x elseif( i == 9)then pz = x elseif( i == 12)then svec(1) = x elseif( i == 13)then svec(2) = x elseif( i == 14)then svec(3) = x elseif( i == 15)then acos_sp = x if(index(dir_file,'INFLECTOR') /=0 .or. index(dir_file,'Start')/=0 .or. index(dir_file,'BACKLEG')/=0)then pvec(1:2)=(/ang, ang_y/) pvec(3) = sqrt((1+pz)**2-pvec(1)**2-pvec(2)**2) spin= svec if(dot_product(spin,spin)>0)then call RightOrLeft(pvec,spin,sxp) sdotp= dot_product(svec,pvec) acos_sp = atan2(sxp(2),sdotp) if(acos_sp < 0)acos_sp = acos_sp + twopi endif endif endif end do muon_alive(n)%muon_numb = muon_number muon_alive(n)%time = time muon_alive(n)%pz = pz muon_alive(n)%x= disp muon_alive(n)%yp = ang_y muon_alive(n)%y= disp_y muon_alive(n)%xp= ang muon_alive(n)%spin = svec muon_alive(n)%phi = acos_sp n_left = n end do 99 continue close(unit=5) return end subroutine get_data_from_END_phase_space subroutine save_and_concatenate(dir_file, phase_space_file) use sim_utils implicit none character*300, intent(in) :: dir_file character*600 long_string character*300 names(40), latest character*300, intent(in) :: phase_space_file logical itexists, new/.true./ integer lun, sofar/0/ integer i integer reason ! are we opening a new file or adding to existing latest = phase_space_file new=.true. do i=1,sofar !loop to see if this combination has already been started if(trim(latest) == trim(names(i)) )then new=.false. exit endif end do if(new)then !If it hasn't been started prepare to write a new , get rid of old version call execute_command_line('rm all_'//trim(latest)) !get rid of old version sofar=sofar+1 names(sofar) = latest endif lun=lunget() open(unit=lun,file='all_'//trim(latest),access='append') !open the file and prepare to write by adding on to whats already there inquire (file=dir_file, exist = itexists) !See if the file to be appended is actually in the subdirectory if(.not. itexists)return ! print '(2a)', 'open ',dir_file open(unit=5,file=dir_file) do while(.true.) read(5,'(a)',IOSTAT=reason)long_string if(reason<0)exit if(.not. new .and. (index(long_string,'muon')/=0 .or. index(long_string,'#')/=0 .or.index(long_string,'!')/=0 ))cycle write(lun,'(a)')trim(long_string) !start adding to whats already there end do close(unit=5) close(lun) return end subroutine save_and_concatenate subroutine fit_line(muon_save, k, slope, offset) use sim_utils use energy_vs_time_mod use energy_vs_time_interface, dummy_except => fit_line implicit none type(muon_alive_struct), allocatable :: muon_save(:) integer k,ix integer lun real(rp) slope, offset real(rp) sum_xp, sum_pp, sum_x, sum_p,sum_xsteps, det lun=lunget() open(unit=lun,file='line_fit.dat') ! then try x vs p sum_xp = sum( muon_save(:)%pz * muon_save(:)%x) sum_pp = sum(muon_save(:)%pz*muon_save(:)%pz) sum_x = sum(muon_save(:)%x) sum_p = sum(muon_save(:)%pz) det = k * sum_pp - sum_p**2 slope = (k*sum_xp - sum_x*sum_p)/det offset = (-sum_p * sum_xp + sum_pp * sum_x)/det print '(4(a,es12.4),a,i10)',' sum_xp =',sum_xp,' sum_pp =',sum_pp,' sum_x =',sum_x,' sum_p=',sum_p,' k=', k print '((a,es12.4))','det =', det print '(a,es12.4)', 'slope = ', slope, '; offset = ', offset write(lun,'(2(a,es12.4))') 'slope_x =',slope,'; offset_x= ', offset ! then try xp vs p sum_xp = sum( muon_save(:)%pz * muon_save(:)%xp) sum_pp = sum(muon_save(:)%pz*muon_save(:)%pz) sum_x = sum(muon_save(:)%xp) sum_p = sum(muon_save(:)%pz) det = k * sum_pp - sum_p**2 slope = (k*sum_xp - sum_x*sum_p)/det offset = (-sum_p * sum_xp + sum_pp * sum_x)/det print '(4(a,es12.4),a,i10)',' sum_xp =',sum_xp,' sum_pp =',sum_pp,' sum_x =',sum_x,' sum_p=',sum_p,' k=', k print '((a,es12.4))','det =', det print '(a,es12.4)', 'slope = ', slope, '; offset = ', offset write(lun,'(2(a,es12.4))') 'slope_xp =',slope,'; offset_xp= ', offset close(lun) return end subroutine fit_line subroutine find_all_files(dir_list , phase_space_file,ix, all_there) use sim_utils implicit none character*(*) dir_list character*300, allocatable ::phase_space_file(:) character*100 name(200) integer i, ix,j, tot_files,lun integer reason logical all_there all_there = .true. lun = lunget() open(unit=lun,file=trim(dir_list), status='old') j=1 do while(.true.) read(lun,'(a)', IOSTAT=reason)name(j) if(reason<0)exit j=j+1 end do tot_files = j close(lun) do 10 i=1,ix do j=1,tot_files if( trim(phase_space_file(i)) == trim(name(j)) )goto 10 end do all_there = .false. print '(a,a,a)', trim(phase_space_file(i)),' not found in directory ', dir_list(1:15) return 10 end do return end subroutine find_all_files ! combine spin_vs_momentum_at_time.dat files subroutine combine_spin_vs_mom_at_time_data(dir_file) use sim_utils implicit none character*300 dir_file character*300 string character*100 file_name_string, file_name(1:10) integer tot_files,lun integer, save :: lun_num(1:10) integer reason integer new_file integer ind integer unit_jj integer, save :: numb(20,-40:40) integer k real(rp), save :: time_new_file(1:10) real(rp) time, momentum, phase, rms, err_avg,fitted_phi real(rp), save :: avg_phase(20,-40:40),error(20,-40:40) logical first/.true./ if(first)then lun_num = (/501,502,503,504,505,506,507,508,509,510/) time_new_file = 0. first=.false. numb(:,:)=0 avg_phase(:,:)=0 error(:,:)=0 end if lun = lunget() new_file = 0 open(unit=lun,file=trim(dir_file)) do while(.true.) read(lun,'(a)', IOSTAT=reason)string if(reason<0)exit if(index(string,'slope')==0 )cycle do while(.true.) read(lun,'(a)', IOSTAT=reason)string if((index(string,"#")/=0) .and. index(string,'slope') == 0)exit if(string(1:10) == ' ')exit read(string, * )time, ind, momentum, phase, rms, err_avg, fitted_phi if(time > 26.e-6 .and. time < 26.1e-6)then !write a file if( all(time_new_file(:) /= time) )then new_file = new_file + 1 time_new_file(new_file) = time write(file_name_string,'(es12.4)')time file_name(new_file) = 'spin_vs_momentum_'//file_name_string(3:12)//'.dat' print '(a,es12.4,a,a)','time =',time,' file name = ', file_name(new_file) open (unit=lun_num(new_file), file= file_name(new_file)) write(lun_num(new_file),'(a1,a12,a10,5a12,a10,2a12)')'#','time','ind','momentum','phase','rms','err_avg','fitted_phi','number','avg_phase','sigma' else unit_jj = find_location_real(time_new_file,time) k = momentum*10000 numb(unit_jj,k) = numb(unit_jj,k)+1 avg_phase(unit_jj,k) = avg_phase(unit_jj,k)+phase error(unit_jj,k) = error(unit_jj,k) + (phase-avg_phase(unit_jj,k)/numb(unit_jj,k))**2 write(lun_num(unit_jj),'(e12.4,i10,5e12.4, i10, 2e12.4)')time, ind, momentum, phase, rms, err_avg, & fitted_phi , numb(unit_jj,k), avg_phase(unit_jj,k)/numb(unit_jj,k), sqrt(error(unit_jj,k)/numb(unit_jj,k)) endif endif end do end do close(unit=lun) return end subroutine combine_spin_vs_mom_at_time_data !example namelist input file ! &energy_spin_vs_time_input ! nargs = 22 ! dir_min = '00000000_000000' ! dir_max = '99999999_999999' ! phase_space_fille(1) = 'DistStartInjLine_phase_space.dat' ! phase_space_file(2) = 'MARK_INFLECTOR_DS_phase_space.dat ! phase_space_file(3) ='phase_space_at_01098_ns.dat' ! phase_space_file(4) ='phase_space_at_02216_ns.dat' ! phase_space_file(5) ='phase_space_at_03334_ns.dat' ! phase_space_file(6) ='phase_space_at_04452_ns.dat' ! phase_space_file(7) ='phase_space_at_05570_ns.dat' ! phase_space_file(8) ='phase_space_at_06688_ns.dat' ! phase_space_file(9) ='phase_space_at_07806_ns.dat' ! phase_space_file(10) ='phase_space_at_08924_ns.dat' ! phase_space_file(11) ='phase_space_at_10042_ns.dat' ! phase_space_file(12) ='phase_space_at_11160_ns.dat' ! phase_space_file(13) ='phase_space_at_12278_ns.dat' ! phase_space_file(14) ='phase_space_at_13396_ns.dat' ! phase_space_file(15) ='phase_space_at_14514_ns.dat' ! phase_space_file(16) ='phase_space_at_15632_ns.dat' ! phase_space_file(17) ='phase_space_at_16750_ns.dat' ! phase_space_file(18) ='phase_space_at_17868_ns.dat' ! phase_space_file(19) ='phase_space_at_18986_ns.dat' ! phase_space_file(20) ='phase_space_at_20104_ns.dat' ! phase_space_file(21) ='phase_space_at_21222_ns.dat' ! phase_space_file(22) ='phase_space_at_22340_ns.dat' !/