! WARNING: THIS CODE IS DEPRECATED AS OF 2016.01.08 --JSh ! 2009.10.15 - j. shanks ! program to simulate all measurement data ! using the same set of misalignments / errors program sim_all use bmad use random_mod use dr_misalign_mod use radiation_mod use cesr_basic_mod use sim_utils use nonlin_bpm_mod use sim_bpm_mod use sim_decay_mod use mode3_mod use tune_tracker_mod use sim_ac_and_cbar_mod implicit none type (lat_struct), target :: ring type (ele_struct), pointer :: ele type (db_struct) cesr type (coord_struct), allocatable :: orb(:) logical err integer i, j, ix, jx, kx, dummy real(rp) harvest integer :: ran_seed = 0. type(meas_beam_struct) :: beam real(rp) :: current = 1000000, skq02w_kick = 0.! amount of horizontal shearing on bpm's, in meters integer :: n_bpms = 0 character(40) :: bpm_mask = "^DET\_[0-9]{2}[ewEW]$" ! default to using CESR naming convention ! For ORM data simulation: type(orm_data_struct), allocatable :: orm_data(:) integer, allocatable :: hkick_ix(:), vkick_ix(:) ! array of element numbers of the steering magnets integer :: n_steering_tot = 0, n_steering_h = 0, n_steering_v = 0 !break into horizontal and vertical, for tallying real(rp) :: orm_mag_hkick = 1.e-7, orm_mag_vkick = 1.e-7 !magnitude of the kick; identical at each kicker ! for manual dispersion measurement simulation real(rp), allocatable :: eta(:,:) real(rp) :: delta, etaysq(1:99) real(rp) :: etay_rms = 0. ! For gain map data simulation: type(gain_map_data_struct), allocatable :: gain_map_data(:) real(rp) :: dx = 1.e-3, dy = 1.e-3 ! gain map grid spacing integer :: nx = 3, ny = 3 !nx, ny = number of x,y data points; dx,dy = grid spacing ! For TBT tracking: real(rp) :: init_vec(6)=0. ! for starting oscillations for phase 'measurement' logical :: track6x6 = .false., damping = .false. integer :: n_turns = 1024, n_damping = 0 !if 6x6 tracking, n_damping = # of damping times to wait before taking measurement ! Related to phase/coupling simulation real(rp), allocatable :: phi(:,:,:), delta_phi_ax(:), delta_phi_by(:) real(rp) :: dphi_ax_avg = 0, dphi_by_avg = 0 real(rp), allocatable :: cbar_meas(:,:,:), cbar_bmad(:,:,:) ! (:,2,2) = across all BPMs real(rp), allocatable :: A(:,:,:) ! = (n_bpms,2,2) = (i_bpm, (Aa,Ab),(Ba,Bb)) real(rp) :: rms_cbar(2,2)=0, rms_phia = 0, rms_phib = 0 ! parameters for misalignment; to be read in and loaded into bpm_error_sigmas and ma_params structures type(det_error_struct) :: bpm_error_sigmas type(det_struct), allocatable :: bpm(:) type(ma_struct) ma_params(20) logical :: misalign_magnets = .false. real(rp) :: shear_x = 0., gain_sigma = 0., timing_sigma = 0. real(rp) :: bpm_offset = 0., bpm_rotation = 0., bpm_noise = 0. ! parameters for tune tracker logical tunetracker_on TYPE(tt_param_struct) :: tt_params(max_tt) !max_tt=3 set on tune_tracker module INTEGER :: n_TTs = 0 ! number of tune trackers INTEGER id(max_tt) ! File handling: integer ios, version, arg_num, iargc, readstatus character(20) :: cbar_output='cesrv', orm_output='tao_cesr' character(100) lat_file, file_name, lat_file_name, path, output_file_name, rms_file_name character(50) orbit_file namelist /parameters/ lat_file, bpm_offset, bpm_rotation, bpm_noise, gain_sigma, timing_sigma, shear_x, misalign_magnets, & current, dx, dy, nx, ny, orm_output, orm_mag_hkick, orm_mag_vkick, cbar_output, init_vec, & track6x6, damping, n_damping, n_turns, delta, ma_params, skq02w_kick, & n_TTs, tt_params, bpm_mask write(*,'(a)') 'WARNING: THIS CODE IS DEPRECATED AS OF 2016.01.08 --JSh' cbar_output = trim(cbar_output) orm_output = trim(orm_output) !initialize randomizer's seed (harvest) ran_seed = 0. ! ran_seed = 0 <--> seed randomizer w/ current CPU time call ran_seed_put(ran_seed) arg_num=iargc() if(arg_num==0) then file_name='sim_all.in' else call getarg(1, file_name) end if call string_trim (file_name, file_name, ix) open (unit= 1, file = file_name, status = 'old', iostat = ios) if(ios.ne.0)then write(*,*) "" write(*,*) 'ERROR: CANNOT OPEN FILE: ', trim(file_name) endif ! read in the parameters rewind(1) read(1, nml = parameters,iostat=readstatus) if(readstatus > 0) then print *,"CAN NOT READ FROM ",file_name stop end if ! load parameters into bpm_error_sigmas structure: bpm_error_sigmas%bpm_offset = bpm_offset bpm_error_sigmas%bpm_rotation = bpm_rotation bpm_error_sigmas%shear_x = shear_x bpm_error_sigmas%gain_sigma = gain_sigma bpm_error_sigmas%timing_sigma = timing_sigma if ((init_vec(1) .eq. 0.) .and. (init_vec(2) .eq. 0.) .and. (init_vec(3) .eq. 0.) .and. & (init_vec(4) .eq. 0.) .and. (init_vec(5) .eq. 0.) .and. (init_vec(6) .eq. 0.)) then init_vec(1) = 1.e-4 init_vec(3) = 1.e-4 init_vec(5) = 1.e-6 endif ! init lattice if(lat_file(1:8) == 'digested')then call read_digested_bmad_file (lat_file, ring, version) else call fullfilename(lat_file,lat_file) call bmad_parser (lat_file, ring) endif dummy = splitfilename(lat_file,path,lat_file_name) lat_file_name = lat_file_name(1:len_trim(lat_file_name)-4) if (track6x6 .eqv. .true.) then if (n_TTs .gt. 0) then ! use damping and excitation bmad_com%radiation_damping_on = .true. bmad_com%radiation_fluctuations_on = .false. ! tune trackers will dominate well over these effects else bmad_com%radiation_damping_on = .false. bmad_com%radiation_fluctuations_on = .false. endif call set_on_off(rfcavity$, ring, on$) call twiss3_at_start(ring,err) ! twiss3_at_start, twiss3_propagate_all return zero beta at L0? call twiss3_propagate_all(ring) call closed_orbit_calc(ring, orb, 6) call twiss_at_start(ring) call twiss_propagate_all(ring) else call set_on_off(rfcavity$, ring, off$) bmad_com%radiation_damping_on = .false. bmad_com%radiation_fluctuations_on = .false. call twiss_at_start(ring) call twiss_propagate_all(ring) call closed_orbit_calc(ring, orb, 4) endif ! Locate BPMs we wish to study: call find_bpms(ring, bpm_mask, bpm) n_bpms = ubound(bpm,1) beam%current = current beam%I0 = current ! find all horiz, vert steerings: call find_steerings(cesr,hkick_ix,vkick_ix) n_steering_h = ubound(hkick_ix,1) n_steering_v = ubound(vkick_ix,1) n_steering_tot = n_steering_h + n_steering_v call resolution_to_button_error(bpm(1), beam%current, bpm_noise) bpm(:)%butn_res_sigma = bpm(1)%butn_res_sigma if ((n_TTs .gt. 0) .and. (track6x6 .eqv. .true.)) then ! Check if the element for each tune tracker's kicker has kick attributes. call check_tt_kickers(ring,tt_params,n_TTs) ! initialize each tune tracker call tt_init(ring, orb, tt_params, n_TTs, id) endif !=== Generate BPM errors ============================================= ! Note that putting values here have no effect on tracking. They're ! only for storage open(unit=47, file='bpms.ma', status='replace') write(47, '(a7,7a14)') "!index", "x-offset", "y-offset", "tilt", "g1", "g2", "g3", "g4" write(47, '(a7, 7a14)') "" ! call dr_misalign_...(blah) ! directly from ring_ma: write(*,*) "Introducing misalignments..." dr_misalign_params%sigma_cutoff = 3. dr_misalign_params%alignment_multiplier = 1 if (misalign_magnets .eqv. .true.) then call dr_misalign(ring, ma_params) endif do i=1,n_bpms call bpm_errors(bpm(i),bpm_error_sigmas, ring%ele(bpm(i)%ix_quad)) write(47, '(i7,7e14.4)') bpm(i)%ix_db, bpm(i)%x_offset, & bpm(i)%y_offset, bpm(i)%tilt, bpm(i)%gain(:) enddo close(unit=47) open(unit=31, file='quad.ma', status='replace') open(unit=32, file='sext.ma', status='replace') open(unit=33, file='bend.ma', status='replace') write (31, '(a7,a7,3a14)') "name", "index", "x-offset", "y-offset", "tilt" write (32, '(a7,a7,3a14)') "name", "index", "x-offset", "y-offset", "tilt" write (31, '(a7,a7,3a14)') "" write (32, '(a7,a7,3a14)') "" write (33, '(a7,a7,2a14)') "name", "index", "x-offset", "y-offset" write (33, '(a7,a7,3a14)') "" do i=0, 99 ! write quad misalignments to output if (cesr%quad(i)%ix_lat == 0) cycle ele => ring%ele(cesr%quad(i)%ix_lat) write(31, '(a7,3e14.5)') ele%name, value_of_attribute(ele, 'X_OFFSET', err), & value_of_attribute(ele, 'Y_OFFSET', err),value_of_attribute(ele, 'TILT', err) enddo do i=0, 99 ! write sextupole misalignments to output if (cesr%sex(i)%ix_lat == 0) cycle ele => ring%ele(cesr%sex(i)%ix_lat) write(32, '(a7,3e14.5)') ele%name, value_of_attribute(ele, 'X_OFFSET', err), & value_of_attribute(ele, 'Y_OFFSET', err) enddo do i=1, ring%n_ele_max if (match_reg(ring%ele(i)%name, 'RAW_BET.*') .or. match_reg(ring%ele(i)%name, 'BUMP.*')) cycle if (match_reg(ring%ele(i)%name, '^B[0-9]*')) then ele => ring%ele(i) write(33, '(a7,2e14.5)') ele%name, value_of_attribute(ele, 'X_OFFSET', err), & value_of_attribute(ele, 'Y_OFFSET', err) endif enddo close(unit=31) close(unit=32) close(unit=33) !===================================================================== !==============BEGIN GAIN MAP SIMULATION============================== !===================================================================== open(unit=23, file="bpm_gain_anal.input", status="replace") write(23,'(a)') "&bpm_gain_anal_params" write(23,'(a)') " read_current_from_butns_file = .false." write(23,'(a)') " read_fake_data = .true." write(23,'(a)') " sigma_cutoff = 100" write(23,'(a)') " chisq_cutoff = 100" ! write(23,'(a,i3)') " fit_bpm_idx(",ix, " = ", bpm%ix_db !neglecting this input = fit all BPMs. do ix=1,n_bpms if (.not. allocated(gain_map_data)) allocate(gain_map_data(nx*ny)) call sim_gain_map_data(bpm(ix), beam, nx, ny, dx, dy, gain_map_data) do j=1, nx*ny write(23, '(a,i2,a,i2)') " orbit(",j,")%butns_index =", j write(23, '(a,i2,a,es12.3)') " orbit(",j,")%current = ", beam%current write(23,'(a,i2,a,i2,a,4f12.3)') " orbit(",j,")%butns%det(", bpm(ix)%ix_db,")%amp = ", & gain_map_data(j)%amp(1:4) write(23,'(a,i2,a,i2,a,es12.3)') " orbit(",j,")%butns%det(",bpm(ix)%ix_db,")%vec(1) = ", & gain_map_data(j)%vec(1) write(23,'(a,i2,a,i2,a,es12.3)') " orbit(",j,")%butns%det(",bpm(ix)%ix_db,")%vec(3) = ", & gain_map_data(j)%vec(3) enddo enddo write(23,'(a)') "/" close(23) !===================================================================== !============== END GAIN MAP SIMULATION ============================ !============== BEGIN ORM SIMULATION ============================ !===================================================================== ! write mdsa_fakedata.in for tao_cesr if(orm_output == 'tao_cesr') then open(unit=47,file='mdsa_fakedata.in', status='replace') write(47,'(a,a,a)') "&lattice file_name = '$BMAD_LAT/", trim(lat_file_name),".lat' /" do i=1, n_steering_tot write(47,'(a,i3.3,a,i3.3,a)') "&measurement who = orbit 'orbit_files/orbit.1",2*i,"' 'orbit_files/orbit.1",2*i-1, "' /" enddo close(unit=47) endif ! prep orm_data_struct for passing information to sim_orm_data: if (.not. allocated(orm_data)) allocate(orm_data(n_steering_tot*2)) do jx = 1, ubound(hkick_ix,1) orm_data(2*jx-1:2*jx)%mag_hkick = orm_mag_hkick orm_data(2*jx-1:2*jx)%mag_vkick = 0. orm_data(2*jx-1:2:jx)%hkick_ix = hkick_ix(jx) orm_data(2*jx-1:2:jx)%ele_name = ring%ele(hkick_ix(jx))%name enddo do jx = 1, ubound(vkick_ix,1) kx = ubound(hkick_ix,1) + jx orm_data(2*kx-1:2*kx)%mag_hkick = 0. orm_data(2*kx-1:2*kx)%mag_vkick = orm_mag_vkick orm_data(2*kx-1:2:kx)%vkick_ix = vkick_ix(jx) orm_data(2*kx-1:2:kx)%ele_name = ring%ele(vkick_ix(jx))%name enddo ! simulate all ORM data-- save in orm_data_struct call sim_orm_data(ring, bpm, beam, orm_data) do i=1, 2*n_steering_tot write(orbit_file, fmt='(A,i3.3)') 'orbit_files/orbit.1', i open(unit=47, file=orbit_file, status='replace') if(orm_output=='tao_cesr') then !write headers for orbit files write(47,'(a)') "&data_parameters ! Fortran90 namelist." write(47,'(a)') "param%data_type = 'ORBIT' ! Identifies file as an orbit file." write(47,'(a,a,a)') "param%lattice = '", trim(lat_file_name), "'" !write(47,'(a)') "param%data_date = '1998-JAN-06 18:05:36'" !write(47,'(a)') "param%comment = 'This optional comment is used for a title when plotting'" if(i .lt. n_steering_h*2+1) then write(47,'(a,a,a)') "param%var_attrib_name = '", "HKICK", "'" write(47,'(a,a,a)') "param%var_ele_name = '", trim(orm_data(i)%ele_name),"' !name of element being varied" write(47,'(a,es14.4,a)') "param%dvar = ", orm_mag_hkick , " ! radians" else write(47,'(a,a,a)') "param%var_attrib_name = '", "VKICK", "'" write(47,'(a,a,a)') "param%var_ele_name = '", trim(orm_data(i)%ele_name),"' !name of element being varied" write(47,'(a,es14.4,a)') "param%dvar = ", orm_mag_vkick , " ! radians" endif write(47,'(a)') "/" write(47,'(a)') "" write(47,'(a)') "&cesr_data" endif do j=1, n_bpms if(orm_output=='loco') then ! output for LOCO is in meters write(47,'(i6, 2f12.7)') bpm(j)%ix_db, orm_data(i)%bpm(j)%vec(1), orm_data(i)%bpm(j)%vec(3) elseif(orm_output=='tao_cesr') then ! output for tao_cesr must be in meters write(47,'(a,i2,a,2es18.5)') " orbit(", bpm(j)%ix_db, ") =", orm_data(i)%bpm(j)%vec(1), orm_data(i)%bpm(j)%vec(3) endif enddo if(orm_output=='tao_cesr') then write(47,'(a)') "/" write(47,'(a)') "" endif close(unit=47) enddo ! writing ORM data to files. !===================================================================== !============== END ORM SIMULATION ============================ !==============BEGIN PHASE/CBAR SIMULATION============================ !===================================================================== call string_trim (file_name, file_name, ix) file_name = file_name(:ix-3) write(rms_file_name,'(a4,a,a4)') 'rms_', trim(file_name), '.out' open(unit=32,file=trim(rms_file_name)) write(32,'(a)') '!kx cbar22 cbar12 cbar11 phi_a phi_b' allocate(cbar_bmad(n_bpms,2,2), delta_phi_ax(n_bpms), delta_phi_by(n_bpms)) call sim_cbar_data(ring, bpm, n_turns, track6x6, n_TTs, tt_params, id, n_damping, init_vec, phi, cbar_meas, A) open(unit=34, file='amplitude.out', status='replace') write(34,'(a7,8a18)') 'bpm_ix', 'Aa', 'Ab', 'Ba', 'Bb', 'Aa/sqrt(beta.a)', 'Ab/sqrt(beta.b)', & 'Ba/sqrt(beta.a)', 'Bb/sqrt(beta.b)' do ix=1,n_bpms write(34,'(i7,8es18.4)') bpm(ix)%ix_db, A(ix,1,1), A(ix,1,2), A(ix,2,1), A(ix,2,2), & A(ix,1,1) / sqrt(ring%ele(bpm(ix)%ix_lat)%a%beta), & A(ix,1,2) / sqrt(ring%ele(bpm(ix)%ix_lat)%b%beta), & A(ix,2,1) / sqrt(ring%ele(bpm(ix)%ix_lat)%a%beta), & A(ix,2,2) / sqrt(ring%ele(bpm(ix)%ix_lat)%b%beta) enddo close(34) do ix = 1, n_bpms call c_to_cbar(ring%ele(bpm(ix)%ix_lat),cbar_bmad(ix,:,:)) enddo rms_cbar(:,:) = 0 do ix=1,n_bpms if (bpm(ix)%ix_db == 50) cycle ! veto for now rms_cbar(1,1) = rms_cbar(1,1) + (cbar_bmad(ix,1,1)-cbar_meas(ix,1,1))**2 rms_cbar(1,2) = rms_cbar(1,2) + (cbar_bmad(ix,1,2)-cbar_meas(ix,1,2))**2 rms_cbar(2,2) = rms_cbar(2,2) + (cbar_bmad(ix,2,2)-cbar_meas(ix,2,2))**2 enddo rms_cbar(1,1) = sqrt(rms_cbar(1,1)/(n_bpms-1)) ! n_bpms - 1 because we veto BPM 50 for now. rms_cbar(1,2) = sqrt(rms_cbar(1,2)/(n_bpms-1)) rms_cbar(2,2) = sqrt(rms_cbar(2,2)/(n_bpms-1)) ! Correct for the offset in phase from the 'design' phase by adding a constant shift dphi_avg: delta_phi_ax(:) = ring%ele(bpm(:)%ix_lat)%a%phi - phi(:,1,1) delta_phi_by(:) = ring%ele(bpm(:)%ix_lat)%b%phi - phi(:,2,2) dphi_ax_avg = (sum(delta_phi_ax(1:49))+sum(delta_phi_ax(52:))) / (size(delta_phi_ax)-2) ! VETOING BPM 50. indices are CORRECT dphi_by_avg = (sum(delta_phi_by(1:49))+sum(delta_phi_by(52:))) / (size(delta_phi_by)-2) ! VETOING BPM 50. rms_phia = 0 rms_phib = 0 do i=1, n_bpms if (bpm(i)%ix_db == 50) cycle ! veto for now rms_phia = rms_phia + ((phi(i,1,1)+dphi_ax_avg) - ring%ele(bpm(i)%ix_lat)%a%phi)**2 rms_phib = rms_phib + ((phi(i,2,2)+dphi_by_avg) - ring%ele(bpm(i)%ix_lat)%b%phi)**2 enddo rms_phia = sqrt(rms_phia/(n_bpms-1)) rms_phib = sqrt(rms_phib/(n_bpms-1)) write(32,'(i10,5es18.4)') kx, rms_cbar(2,2), rms_cbar(1,2), rms_cbar(1,1), rms_phia, rms_phib write(*,'(a,es10.3,a,es10.3,a,es10.3)') " RMS cbar22 = ", rms_cbar(2,2), & " cbar12 = ", rms_cbar(1,2), " cbar11 = ", rms_cbar(1,1) write(*,'(a19,es10.3,a23,es10.3)') "RMS delta(phi_a) = ", rms_phia, " RMS delta(phi_b) = ", rms_phib ! enddo ! matches [do kx=1, nIter] close(32) !=============================================================== open(unit=74,file='cbar.out',status='replace') write(74,'(a7,6a18)') '!bpm_ix', 'cbar12 (meas)', 'cbar12 (bmad)', 'cbar22 (meas)', & 'cbar22 (bmad)', 'cbar11 (meas)', 'cbar11 (bmad)' do ix=1,n_bpms if (bpm(ix)%ix_db == 50) cycle ! debug-- veto BPM 50 for now. nonlin_bpm model is bad there? write(74,'(i7,6es18.4)') bpm(ix)%ix_db, cbar_meas(ix,1,2), cbar_bmad(ix,1,2), & cbar_meas(ix,2,2), cbar_bmad(ix,2,2), cbar_meas(ix,1,1), cbar_bmad(ix,1,1) enddo close(74) ! open(unit=72,status='replace') open(unit=73, file='phase.out',status='replace') write(73,'(a7,6a18)')"bpm_ix","phi_ax","phi_by","phi_ax (des)","phi_by (des)", "phi_ax (meas-des)", "phi_by (meas-des)" do i=1,n_bpms if (bpm(i)%ix_db==50) cycle ! skip bpm 50 for now.... write(73,'(i7,6es18.4)') bpm(i)%ix_db, phi(i,1,1)+dphi_ax_avg, phi(i,2,2)+dphi_by_avg, & ring%ele(bpm(i)%ix_lat)%a%phi, ring%ele(bpm(i)%ix_lat)%b%phi, & phi(i,1,1)+dphi_ax_avg-ring%ele(bpm(i)%ix_lat)%a%phi, & phi(i,2,2)+dphi_by_avg-ring%ele(bpm(i)%ix_lat)%b%phi enddo close(73) close(72) write(output_file_name,'(a,a4)') trim(file_name),'.dat' !write(*,*) output_file_name open(unit=41,file=output_file_name, status='replace') if (cbar_output == 'cesrv') then write(41,'(a)') " &DATA_PARAMETERS ! Fortran90 namelist." write(41,'(a)') " file_type = 'PHASE DATA' ! Identifies file as an phase/cbar file." write(41,'(a,es18.4)') " horiz_shake_freq = ", ring%a%tune write(41,'(a,es18.4)') " horiz_shake_freq = ", ring%b%tune ! write(41,'(a)') " lattice = 'N9A18A000.FD93S_5269_15KG' ! This is optional." ! write(41,'(a)') " data_date = '1998-JAN-06 18:05:36' ! This is optional." ! write(41,'(a)') " save_set = 51220 ! This is optional." ! write(41,'(a)') " comment = 'This optional comment is used for a title when plotting'" write(41,'(a)') " /END" write(41,'(a)') "" write(41,'(a)') "! det phase_a phase_b cbar_11 cbar_12, cbar_22" write(41,'(a)') "! (radians) (radians)" do i=1,n_bpms if(bpm(i)%ix_db==50) cycle ! veto 50 for now... write(41,'(i7,5es18.4)') bpm(i)%ix_db, phi(i,1,1)+dphi_ax_avg, phi(i,2,2)+dphi_by_avg, & cbar_meas(i,1,1), cbar_meas(i,1,2), cbar_meas(i,2,2) enddo elseif (cbar_output == 'tao') then write(41,'(a)') "set global lattice_calc_on = F" write(41,'(a)') "set global plot_on = F" do i=1, n_bpms if (bpm(i)%ix_db==50) cycle ! veto 50 for now... bad nonlin_bpm map? write(41, '(a17,i2,a9,es14.4)') "set data phase.a[",bpm(i)%ix_db,"]|meas = ", phi(i,1,1) +dphi_ax_avg write(41, '(a17,i2,a9,es14.4)') "set data phase.b[",bpm(i)%ix_db,"]|meas = ", phi(i,2,2) +dphi_by_avg write(41, '(a17,i2,a9,es14.4)') "set data cbar.11[",bpm(i)%ix_db,"]|meas = ", cbar_meas(i,1,1) write(41, '(a17,i2,a9,es14.4)') "set data cbar.12[",bpm(i)%ix_db,"]|meas = ", cbar_meas(i,1,2) write(41, '(a17,i2,a9,es14.4)') "set data cbar.22[",bpm(i)%ix_db,"]|meas = ", cbar_meas(i,2,2) enddo write(41,'(a)') "set global lattice_calc_on = T" write(41,'(a)') "set global plot_on = T" write(41,'(a)') "place top phase" write(41,'(a)') "place bottom cbar" write(41,'(a)') "plot top meas-design" write(41,'(a)') "plot bottom meas-design" endif close(41) !================================================================================= !==================END PHASE/CBAR SIMULATION====================================== !==================BEGIN ETA SIMULATION=========================================== !================================================================================= open(unit=31, file='rms_etay.out', status='replace') open(unit=32, file='eta.out', status='replace') write(31,'(a6,4a15)') 'kx', 'bpm_noise', 'bpm_rotation', 'gain_sigma', 'etay_rms(kx)' write(32,'(a6,5a15)') 'bpm_ix', 'bpm_noise', 'bpm_rotation', 'gain_sigma', 'etax(ix)', 'etay(ix)' ! if (.not. allocated(etay_rms)) allocate(etay_rms(n_iter)) call sim_eta_data(ring, bpm, delta, eta) do ix = 1, n_bpms write(32,'(i6,5es15.3)') ix, bpm_noise, bpm_rotation, gain_sigma, eta(ix,1), eta(ix,2) etaysq(ix) = (eta(ix,2))**2 enddo etay_rms = sqrt(sum(etaysq)/n_bpms) ! RMS write(31,'(i6,4es15.3)') kx, bpm_noise, bpm_rotation, gain_sigma, etay_rms close(31) close(32) !================================================================================= !===================END ETA SIMULATION============================================ !================================================================================= end program sim_all