! Compile data from subdirectories labelled as ! Read ,, and collect data ! Summarize data in three files ! - selected data from and ! - selected data from and ! program compile_input use parameters_bmad use quad_scrape_parameters use muon_mod implicit none type (g2twiss_struct) twiss, twiss_opt, twiss_match, twiss_start_inj_line type (initial_offsets_struct) initial_offsets, initial_offsets_ref, inflector_end_target type (kicker_params_struct) kicker_params, kicker_params_0 type (field_file_struct) fringe_file, inflector_file integer lun,lun2,lun3,lun1,lun4,lun5 integer end, nturns, nmuons, turn_that_counts/10/ integer track_state integer seed/1/ integer vparam_id/0/, vparam_cbo_turns integer ios integer n/0/ integer m integer icolumn(100,1:38) integer ikick_max integer column1 integer ix integer p integer tot_turns/0/, tot_stored/0/ real(rp) xmean_m, xsigma_m, ymean_m, ysigma_m, tmean, tsigma, pxmean, pxsigma, pymean, pysigma, pzmean, pzsigma real(rp) epsx, epsy, tlength, pz real(rp) time_bin_width/0./ real(rp) ring_theta/0./ real(rp) vparam, vparam_max, delta_vparam, vparam_min real(rp) inflector_angle real(rp) column(100,1:38) real(rp) s,x,xp,y,yp,z real(rp) x_avg, sig_x, sig_p, p_ave real(rp) x_cbo, max_x/-99./, min_x/99./ real(rp) freq, betaCrossE(1:3),sum_time,BetaDotB(1:3), B_field(1:3), cyc_freq, BetaCrossB(1:3), h_freq, v_freq real(rp) x_turn_avg character*120 line, lat_file, lat_file_name/''/,new_file/' '/, muon_file/' '/ character*300 string, new_string character*300 long_string, another_long_string character*15 dir character*50 dir_file, moments_file character*16 element_name character *120 azimuthal_exp_z, azimuthal_exp_r character *120 azimuthal_exp_z_nom/'BzFourier20170628_LogID983.dat'/, azimuthal_exp_r_nom/'BrFourier2016.dat'/ logical err_flag, inf_end_us/.false./, inf_end_ds/.false./, quad_plate/.false./ logical first/.true./, twiss_file/.false./ logical create_new_distribution/.false./ logical loop/.false./,gnu_input_only/.false./ logical enerloss/.false./, quad_p/.false./ logical spin_tracking_on/.true./,muon_decay_on/.true./ logical local_ref/.false./,calcd/.false./ !why was this removed logical opt_incident/.false./,inj_matrix_tracking/.false./ logical start_tracking_at_inflector_exit/.false./ logical make_movie/.false./ logical use_lattice_twiss/.false./ logical write_phase_space_file/.false./ logical save_element_by_element_info/.true./ logical FiberScatter0/.false./ logical error_fields/.false./ logical itexists character*16 epsdistr, tdistr, pzdistr, inf_aperture, twiss_ref, ring_twiss !stuff for reading Beam_moments_tbt type (g2moment_struct) moment integer number_remaining(0:10000) integer turn, nmu integer k, max_turn integer reason real(rp) moment_ave_sqr, turn_x character*16 ele_name ! end of stuff for reading Beam_moments_tbt ! Input namelist structure ("/g-2/test/input.dat") namelist /input/ lat_file_name, nmuons, nturns, & ! GENERAL create_new_distribution, new_file, muon_file, & ! BEAM (write and/or read from file) tdistr, tlength, tsigma, & ! BEAM (longitudinal) pzdistr, pz, pzsigma, & ! BEAM (energy) epsdistr, epsx, epsy, twiss, twiss_ref, & ! BEAM (transverse) inf_aperture, inf_end_us, inf_end_ds, initial_offsets, initial_offsets_ref, enerloss, inflector_angle, inflector_field, & ! INFLECTOR inflector_end_target, & !target coordinates fro opt incident at end of inflector ring_theta, & ! ELEMENT POSITIONING kickerPlates, kickerCircuit, kickerPulseFile, kickerFieldType, & ! KICKER kicker_params, & ! KICKER kicker_tStart, & ! KICKER quad_plate, & ! QUAD SCATTERING quadPlates, quadCircuit, quadFieldType, & ! QUAD quad_params, & ! QUAD twiss_file, & loop, vparam_max, vparam_id, delta_vparam, vparam_min, vparam_cbo_turns, & gnu_input_only,spin_tracking_on, muon_decay_on, inflector_width, opt_incident, & inj_matrix_tracking, & fringe_file, inflector_file, & ! names of field maps for fringe and inflector and a few other details about the maps start_tracking_at_inflector_exit, & !track distribution from inflector exit rather than start of injection line ring_twiss, & ! ring_twiss = 'open' use input twiss as starting point for propagation around ring. ring_twiss='closed' compute closed ring seed, & make_movie, & use_lattice_twiss, & ! to compute twiss parameters through injection line, if true, start with values defined in lattice file write_phase_space_file, & ! to write a file with muon phase space at each element on first turn save_element_by_element_info, & !if true then allocate array muons_ele(:,:) which saves distribution at each element rf_quad, & !parameters of rf quads time_bin_width, & !bin for tbt moments vs time scraping_on, init_quad_focus, init_quad_steer, quad_ramp_start_time, quad_ramp_end_time, & B_radial, & !turn on quad scraping, scraping focus fraction, scraping steering fraction, start time and end time FiberScatter0, & !if true turn on scattering in harp fibers, default = false error_fields, & ! if true call field_errors to introduce errors in g-2 dipole field, default = false azimuthal_exp_z, azimuthal_exp_r ! fourier expansion of z and r fields scraping_on=.false. kicker_params%dtRise = 0. kicker_params%dtFall = 0. quad_params%short_quad_field_index(1) = -99. quad_params%long_quad_field_index(1) = -99. rf_quad(1:4)%amp_h=0 rf_quad(1:4)%amp_v=0 string='ls > out.dat' call execute_command_line(string) lun=lunget() open(unit=lun,file='out.dat', status='old') lun2=lunget() open(unit=lun2, file='compilation_vparam.dat') write(lun2,'(27a12)')'sub_dir','nmuons','nturns','muon_file','tsigma','tdistr','pzsigma','epsx','epsy','Q1_long','start_inf','rf_freq_h','rf_amp_h','phi_h','rf durat', & 'scraping','init_off%x','init_off%px','init_off%pz','inf_angle','inf_field','inf_exit_x','inf_exit_ang','vparam_id','x_inf','xp_inf','betax' lun3=lunget() open(unit=lun3, file='compilation_vparam_brief.dat') write(lun3,'(13a12)')'sub_dir','nmuons','nturns','init_off%x','init_off%px','init_off%pz','inf_angle','inf_field','inf_exit_x','inf_exit_ang','vparam_id','x_inf','xp_inf' lun1=lunget() open(unit=lun1, file='compilation.dat') write(lun1,'(26a12)')'sub_dir','nmuons','nturns','init_off%x','init_off%px','init_off%pz','inf_angle','inf_field','inf_exit_x','inf_exit_ang','kicker1','kicker2','kicker3','Decay','x_inf','xp_inf', & 'tot_turns','tot_stored','loop','muon_file','new_dist','x_avg','sig_x','sig_p','p_ave','x_cbo' lun5=lunget() open(unit=lun5,file='compile_spin_data.dat') write(lun5,'(a1,a16,9x,a10,6a12,a16,3a16,a16,2(3a16),a16,a16,a16,a16,a16,4a16)')'#','sub_dir','nturns','init_off%x','init_off%y','init_off%z','init_off%px','init_off%py','init_off%pz','omega_f',& 'betaCrossE-x', 'betaCrossE-y', 'betaCrossE-z','sum_time','BetaDotB-x','BetaDotB-y','BetaDotB-z',& 'B_field-x','B_field-y','B_field-z','cyc_freq','lattice', & 'BetaCrossB-x', 'BetaCrossB-y', 'BetaCrossB-z','h-betatron-freq','v-betatron-freq','spin tune','x_turn_avg' do while(.true.) read(lun, '(a)', IOSTAT = reason)new_string if(reason < 0)goto 99 if( (index(new_string,'2017')==0) .and. (index(new_string,'2018')==0) .and. (index(new_string,'2019')==0) .and. & (index(new_string,'2022')==0) .and. (index(new_string,'configuration')==0))cycle dir=trim(new_string) dir_file = dir//'/input.dat' print '(a)',dir_file inquire (file=dir_file, exist=itexists) if (.not.itexists) cycle OPEN (UNIT=5, FILE=dir_file, STATUS='old', IOSTAT=ios) READ (5, NML=input, IOSTAT=ios) print *, 'open: ', dir_file, ios WRITE(6,NML=input) ! print *, 'ios=', ios ! rewind(unit=5) CLOSE(5) x=-99 xp=-99 dir_file = dir//'/track_backleg_inflector.dat' inquire (file = dir_file, exist= itexists) if(itexists)then open(unit = 5, file=dir_file) print *,'open ', dir_file do while(.true.) read(5,'(a)',end=299)long_string print '(a50,a100)',dir_file,long_string(1:100) if(index(long_string,'END')/=0)then read(long_string,*)s,x,xp,y,yp,z,pz,element_name print '(7es12.4,a)',s,x,xp,y,yp,z,pz,element_name endif end do 299 continue close(unit=5) endif dir_file = dir//'/Beam_moments_tbt.dat' inquire (file = dir_file, exist= itexists) long_string = '' if(itexists)then open(unit = 5, file=dir_file) print *,'open ', dir_file max_x=-99. min_x = 99. do while(.true.) read(5,'(a)',end=399)another_long_string long_string = another_long_string read(long_string,*,IOSTAT = ios)column(1,1:16) if(column(1,2) > 10 .and. ios == 0)then max_x = max(max_x,column(1,3)) min_x = min(min_x, column(1,3)) endif x_cbo = max_x - min_x end do 399 continue print '(a50,a160)',dir_file,long_string if(index(long_string,'END')/=0)then read(long_string,*)column(1,1:16) ! print '(i10,1x,14es12.4,1x,i10)',column(1,1:16) tot_turns=column(1,2) tot_stored=column(1,16) x_avg = column(1,3) sig_x = sqrt(column(1,6)) sig_p = sqrt(column(1,14)) p_ave = column(1,15) endif close(unit=5) endif dir_file = dir//'/spin_tunes.dat' call get_data_from_spin_tunes(dir_file,freq, cyc_freq) dir_file = dir//'/BetaCrossE_turn.dat' call get_data_from_BetaCrossE_turn(dir_file,betaCrossE,sum_time,BetaDotB,B_field, BetaCrossB) dir_file = dir//'/betatron_tunes.dat' call get_data_from_betatron_tunes(dir_file,h_freq,v_freq, cyc_freq) dir_file = dir//'/EndOfTracking_phase_space.dat' call get_data_from_EndOfTracking(dir_file, x_turn_avg) m=0 p=0 dir_file = dir//'/VparamDependence.dat' inquire (file=dir_file, exist=itexists) if (itexists) then if( loop)then open(unit=5,file=dir_file) print *,' open ', dir_file, 'loop = ', loop do while(.true.) read(5,'(a)', end=199)long_string print *,long_string(1:92) ix = index(long_string,'column1') ! if(ix /= 0)print *,'column1 = ',long_string(ix+10:) if(ix /= 0)read(long_string(ix+10:),*)column1 m=m+1 ! if(m>=10)read(string,'(3es16.3,6es16.4,4es16.4,3i16,2es12.4,7es12.4, 2es12.4, 5es12.4, 6es12.4)')column(m,1:13),icolumn(m,14:16),column(m,17:38) if(m>=10)then p=p+1 read(long_string,*)column(m,1:13),icolumn(p,14:16), column(p,17:20) print '(2i10,13es12.4,3i10,17es12.4)',column1,p,column(p,1:13),icolumn(p,14:16), column(p,17:20) endif end do 199 continue close(unit=5) ikick_max = maxloc(icolumn(1:p,16),1) print *,' ikick_max = ', ikick_max, icolumn(1:p,16), 'column1= ',column1, column(1:p,column1) n=n+1 write(lun2,'(a15,1x,i10,i8,1x,a14,es12.4,1x,a11,4es12.4,l12,4es12.4,l12,3es12.4,4es12.4,1x,i11, es12.4,1x,i10,1x,3es12.4)') dir,nmuons, nturns, & muon_file, tsigma,tdistr,pzsigma,epsx,epsy,quad_params%long_quad_field_index(1),start_tracking_at_inflector_exit,& rf_quad(2)%freq_h,rf_quad(2)%amp_h,rf_quad(2)%phi_h,rf_quad(2)%t_max,scraping_on, initial_offsets%x_mean, initial_offsets%pxmean,& initial_offsets%pzmean,inflector_angle, inflector_field,inflector_end_target%x_mean,inflector_end_target%pxmean,vparam_id, column(ikick_max,column1), icolumn(ikick_max, 16) , x, xp, column(1,17) write(lun3,'(a15,1x,i10,i8,1x,3es12.4,4es12.4,1x,i11, 1x,es12.4,1x,i10,1x,2es12.4)') dir,nmuons, nturns, & initial_offsets%x_mean, initial_offsets%pxmean,& initial_offsets%pzmean,inflector_angle, inflector_field,inflector_end_target%x_mean,inflector_end_target%pxmean,vparam_id, column(ikick_max,column1), icolumn(ikick_max, 16) ,x, xp if((n/10)*10 == n)then write(lun2,*) write(lun3,*) endif endif !loop endif !VparamDependence exists write(lun1,'(a15,1x,i10,i8,1x,3es12.4,7es12.4,l12,1x,2es12.4,1x,2i12,1x,l6x,a20,1x,l9,1x,5es12.4)') dir,nmuons, nturns, & initial_offsets%x_mean, initial_offsets%pxmean,& initial_offsets%pzmean,inflector_angle, inflector_field,inflector_end_target%x_mean,inflector_end_target%pxmean,kicker_params%kicker_field(1),kicker_params%kicker_field(2),kicker_params%kicker_field(3), & muon_decay_on, x,xp, tot_turns, tot_stored, loop, muon_file, create_new_distribution, x_avg, sig_x, sig_p, p_ave, x_cbo write(lun5,'(a16,1x,i9,6es12.4,11es16.8, es16.8,1x,a, 7es16.8)')dir,nturns,initial_offsets%x_mean, initial_offsets%y_mean, initial_offsets%tmean, & initial_offsets%pxmean, initial_offsets%pymean, initial_offsets%pzmean, & freq,betaCrossE(1:3),sum_time,BetaDotB(1:3),B_field(1:3), cyc_freq, trim(lat_file_name), BetaCrossB, h_freq, v_freq,& freq/cyc_freq, x_turn_avg end do close(unit=lun2) close(unit=lun3) 99 close(unit=lun) end subroutine get_data_from_spin_tunes(dir_file,freq, cyc_freq) use bmad implicit none character*50 dir_file character*200 long_string logical itexists real(rp) freq,a_real,a_img, cyc_freq integer ix inquire (file=dir_file, exist = itexists) freq=0. if(.not. itexists)return print '(2a)', 'open ',dir_file open(unit=5,file=dir_file) do while(.true.) read(5,'(a)',end=99)long_string ix = index(long_string,'frequency') if(ix /=0)read(long_string(ix+13:),*)cyc_freq if(index(long_string,'z-x plane')/=0)then read(5,'(a)',end=99)long_string if(index(long_string,'spin_tune')/= 0)then read(5,*)freq,a_real,a_img exit else read(long_string,*)freq,a_real,a_img exit endif print '(3es12.4)',freq,a_real,a_img endif end do 99 close(unit=5) return end subroutine get_data_from_spin_tunes subroutine get_data_from_betatron_tunes(dir_file,h_freq,v_freq, cyc_freq) use bmad implicit none character*50 dir_file character*200 long_string character plane*1 logical itexists real(rp) h_freq,v_freq,a_real,a_img, cyc_freq integer ix inquire (file=dir_file, exist = itexists) h_freq=0. v_freq=0. if(.not. itexists)return print '(2a)', 'open ',dir_file open(unit=5,file=dir_file) do while(.true.) read(5,'(a)',end=99)long_string ix = index(long_string,'frequency') if(ix /=0)read(long_string(ix+13:),*)cyc_freq ix = index(long_string,'plane') if(ix/=0)then plane = long_string(ix-2:ix-1) read(5,'(a)',end=99)long_string if(index(long_string,'betatron_tune')/= 0)then if(plane == 'x')read(5,*)h_freq,a_real,a_img if(plane == 'y')read(5,*)v_freq,a_real,a_img cycle else if(plane == 'x')read(long_string,*)h_freq,a_real,a_img if(plane == 'y')read(long_string,*)v_freq,a_real,a_img exit endif if(plane == 'x')print '(a,3es12.4)','horizontal freq',h_freq,a_real,a_img if(plane == 'y')print '(a,3es12.4)','vertical freq',v_freq,a_real,a_img endif end do 99 close(unit=5) return end subroutine get_data_from_betatron_tunes subroutine get_data_from_BetaCrossE_turn(dir_file,betaCrossE,sum_time,BetaDotB,B_field, betaCrossB) use bmad implicit none character*50 dir_file character*600 long_string logical itexists real(rp) x,xp,y,yp,z,pz,t,betaCrossE(1:3),sum_time,BetaDotB(1:3), path, pathx, B_field(1:3), betacrossB(1:3) integer turn inquire (file=dir_file, exist = itexists) if(.not. itexists)return print '(2a)', 'open ',dir_file open(unit=5,file=dir_file) do while(.true.) read(5,'(a)',end=99)long_string if(index(long_string,'turn')/=0)cycle ! read(long_string,'(i10,1x,22es12.4)')turn,x,xp,y,yp,z,pz,t,betacrossE(1:3),sum_time,BetaDotB(1:3),path,pathx,B_field(1:3), betaCrossB(1:3) read(long_string,*)turn,x,xp,y,yp,z,pz,t,betacrossE(1:3),sum_time,BetaDotB(1:3),path,pathx,B_field(1:3), betaCrossB(1:3) ! print '(i10,1x,19es12.4)',turn,x,xp,y,yp,z,pz,t,betacrossE(1:3),sum_time,BetaDotB(1:3),path,pathx,B_field(1:3) end do 99 continue close(unit=5) return end subroutine get_data_from_BetaCrossE_turn subroutine get_data_from_EndOfTracking(dir_file,x_avg) use bmad implicit none character*50 dir_file character*600 long_string logical itexists real(rp) x_avg,x integer ix,i inquire (file=dir_file, exist = itexists) if(.not. itexists)return print '(2a)', 'open ',dir_file open(unit=5,file=dir_file) do while(.true.) read(5,'(a)',end=99)long_string if(index(long_string,'')/=0)then read(5,'(a)', end=99)long_string else cycle endif i=0 ix=0 do while (i<24) call string_trim(long_string(ix+1:), long_string, ix) i=i+1 read(long_string(1:ix),*)x if(i < 24)cycle x_avg=x goto 99 end do end do 99 continue close(unit=5) return end subroutine get_data_from_EndOfTracking