subroutine calibrate_bpm (ix_bpm_in, ix_var_in, how, cu_in, u, graph, err_flag, ix_loop, finished) use cesrv_struct use cesrv_interface use cesrv_utils_mod use quick_plot use cesr_read_data_mod implicit none type (var_struct), pointer :: k1_var, x_kick_var, y_kick_var, v(:) type (universe_struct), target :: u type (ele_struct), pointer :: q_ele, bpm_ele type (graph_struct) graph type (data_struct), pointer :: d(:), orb_x, orb_y type (cesr_xy_data_struct), pointer :: bpm_offset integer cu_delta, cu_orig, iv, i, cu1, cu2, cu_in(2), ix_loop integer n, iu, n_dat, ix_var, type, ix_quad, ix_var_in, ix_bpm_in integer ix, ix_ele_quad, ix_db, iv_x, iv_y, ie, ix_bpm integer, save :: iu_raw = 0 real(rp) dk1, cal_meas, cal_old, slope, s_quad real(rp) dmeas, dk1_err, merit0, psum, ds, ds_min real(rp) sum2_meas_ref, sum2_all, rms_x, rms_y, dx, dy, w_factor, l_quad, orb_x_orig, orb_y_orig character(*) how character(12) node_name, bpm_name, ans character(20) date character(140) cal_file, file_name, raw_file_name character(72) comment character(7) ratio_str logical make_meas, err_flag, is_there, does_exist, hit_limit, loaded_steerings logical, optional :: finished logical, save :: init_needed = .true. ! Init if (present(finished)) finished = .true. ! Locate quad / bpm if needed ix_var = ix_var_in ix_bpm = ix_bpm_in if (how == 'RECAL') then ! Find nearest bpm q_ele => u%ring%ele(u%var(ix_var)%ix_ele) if (ix_bpm < 0) then ds_min = 1e10 d => u%orbit%x%d do i = 0, 99 if (.not. d(i)%exists) cycle ds = abs(d(i)%s - (q_ele%s - q_ele%value(l$)/2)) if (ds > u%ring%param%total_length/2) ds = u%ring%param%total_length - ds if (ds < ds_min) then ds_min = ds ix_bpm = i endif enddo else ds = abs(u%orbit%x%d(ix_bpm)%s - (q_ele%s - q_ele%value(l$)/2)) if (ds > u%ring%param%total_length/2) ds = u%ring%param%total_length - ds endif elseif (how == 'NEW_CAL') then ! Find nearest quad if (.not. u%orbit%x%d(ix_bpm)%exists) then print *, 'BPM does not exist!' return endif ds_min = 1e10 v => u%quad_k1%v do i = lbound(v, 1), ubound(v, 1) if (.not. v(i)%exists) cycle q_ele => u%ring%ele(v(i)%ix_ele) ds = abs(u%orbit%x%d(ix_bpm)%s - (q_ele%s - q_ele%value(l$)/2)) if (ds > u%ring%param%total_length/2) ds = u%ring%param%total_length - ds if (ds < ds_min) then ds_min = ds ix_var = v(i)%ix_var endif enddo endif bpm_name = u%orbit%x%d(ix_bpm)%ele_name(5:) if (.not. u%orbit%x%d(ix_bpm)%good_user) then print *, 'BPM is vetoed! Stopping calibration for BPM ', bpm_name return endif if (.not. u%orbit%x%d(ix_bpm)%good_dat) then print *, 'The data is not good for BPM! Stopping calibration for BPM ', bpm_name return endif ! Init if (init_needed) then call read_bpm_quad_offsets (u%bpm_offset) init_needed = .false. endif k1_var => u%var(ix_var) node_name = k1_var%db_node_name ix_db = k1_var%ix_db print '(1x, 5a, i4, a)', 'Quadrupole: ', k1_var%name, '[', node_name, ',', ix_db, ']' print '(1x, 2a)', 'Calibrating BPM: ', bpm_name print '(1x, a, f10.1)', 'Distance between Quad center and BPM:', ds_min if (.not. any(node_name == ['CSR QUAD CUR', 'CSR QADD CUR'])) then print *, 'ERROR IN CALIBRATE_BPM: THIS IS NOT A QUADRUPOLE: ', node_name return endif if (k1_var%dvar_dcu == 0) then print *, 'ERROR IN CALIBRATE_BPM: QUAD CALIBRATION FACTOR IS ZERO!', node_name return endif ix_quad = k1_var%ix_index ix_ele_quad = k1_var%ix_ele q_ele => u%ring%ele(ix_ele_quad) !--------------------------------------------------------------- ! Recalibration using old data. ! In this case ix_in is the varying quad index if (how == 'RECAL') then cu1 = k1_var%cu_saved_ref cu2 = k1_var%cu_saved cu_delta = cu2 - cu1 print *, 'Range used CU(ref)/CU(data): ', cu1, cu2 print *, 'Delta CU: ', cu_delta cu_orig = (cu1 + cu2) / 2 ! Just to have a number !--------------------------------------------------------------- ! This section if we need to take a measurement elseif (how == 'NEW_CAL') then call vxgetn (node_name, ix_db, ix_db, cu_orig) ! Find delta. ! If cu_in(1) /= 0 but cu_in(2) == 0 then cu_in(1) is the delta. ! If cu_in(1) == 0 but cu_in(2) /= 0 then cu_in(2) is percent of nominal delta. if (cu_in(1) /= 0 .and. cu_in(2) /= 0) then cu1 = cu_in(1) cu2 = cu_in(2) else if (cu_in(1) == 0) then dk1 = logic%dQ_max_quad_calib * fourpi / (max(q_ele%a%beta, q_ele%b%beta) * q_ele%value(l$)) cu_delta = abs(dk1 / k1_var%dvar_dcu) if (cu_in(2) /= 0) cu_delta = cu_delta * cu_in(2) / 100.0 else cu_delta = cu_in(1) endif cu_delta = min (cu_delta, k1_var%cu_high_lim - k1_var%cu_low_lim) ! calc what max and min should be. ! try to make max and min not near 0 if (cu_orig + cu_delta/2 > k1_var%cu_high_lim) then cu2 = k1_var%cu_high_lim cu1 = cu2 - cu_delta elseif (cu_orig - cu_delta/2 < k1_var%cu_low_lim) then cu1 = k1_var%cu_low_lim cu2 = cu1 + cu_delta else cu1 = cu_orig - cu_delta / 2 cu2 = cu1 + cu_delta endif endif ! Type out endpoint info cu_delta = cu2 - cu1 call inverse_cross_corr (node_name, cu1, cu2, cu_delta, ix_db, ix_db) print *, 'CU now: ', cu_orig print *, 'Range used CU(ref)/CU(data): ', cu1, cu2 print *, 'Delta CU: ', cu_delta print '(a, 2f10.1)', 'H Tune variation +/- kHz:', & 390.1 * q_ele%value(l$) * q_ele%a%beta * (cu1-cu_orig) * k1_var%dvar_dcu / fourpi, & 390.1 * q_ele%value(l$) * q_ele%a%beta * (cu2-cu_orig) * k1_var%dvar_dcu / fourpi print '(a, 2f10.1)', 'V Tune variation +/- kHz:', & -390.1 * q_ele%value(l$) * q_ele%b%beta * (cu1-cu_orig) * k1_var%dvar_dcu / fourpi, & -390.1 * q_ele%value(l$) * q_ele%b%beta * (cu2-cu_orig) * k1_var%dvar_dcu / fourpi if (max(cu1, cu2) > k1_var%cu_high_lim .or. min(cu1, cu2) < k1_var%cu_low_lim) then print *, 'ERROR IN CALIBRATE_BPM: A SET POINT IS OUTSIDE A LIMIT' print *, ' LIMITS ARE:', k1_var%cu_high_lim, k1_var%cu_low_lim err_flag = .true. return endif ! Measure phase at low point and high point make_meas = .false. if (logic%auto_measurement) then make_meas = .true. else call cesrv_logic_get ('Y', 'N', 'Make an Phase measurement?', make_meas) if (.not. make_meas) then print *, 'No measurement, No Analysis.' print *, 'If you want to analyze old data use the RECALIBRATE command.' return endif endif if (make_meas) then write (comment, '(2a, 4x, 2a, 4x, a, i6)') & 'BPM: ', trim(bpm_name), 'Quad: ', trim(k1_var%name), 'CU:', cu1 call set_and_meas (node_name, ix_db, cu1, phase_data$, (/ix_var/), 1, ref_file$, & u, graph, err_flag, comment) if (err_flag) then print *, 'RESETTING QUADRUPOLE AND QUITTING...' if (.not. logic%debug) call vxputn (node_name, ix_db, ix_db, cu_orig) return endif write (comment, '(2a, 4x, 2a, 4x, a, i6)') & 'BPM: ', trim(bpm_name), 'Quad: ', trim(k1_var%name), 'CU:', cu2 call set_and_meas (node_name, ix_db, cu2, phase_data$, (/ix_var/), 1, data_file$, & u, graph, err_flag, comment) if (err_flag) then print *, 'RESETTING QUADRUPOLE AND QUITTING...' if (.not. logic%debug) call vxputn (node_name, ix_db, ix_db, cu_orig) return endif if (.not. logic%debug) call vxputn (node_name, ix_db, ix_db, cu_orig) endif endif !--------------------------------------------------------------- ! Make sure we have enough data points n = count(u%phase%x%d%exists .and. u%phase%x%d%good_dat .and. u%phase%x%d%good_ref) if (n < 40) then print * print *, 'NOT ENOUGH PHASE DATUMS (< 40) FOR ANALYSIS. STOPPING HERE.' return endif !--------------------------------------------------------------- ! Save current values logic%opt_base_locked = .true. ! So baseline_set will not set this. u%data%val_temp = u%data%meas u%data%good_temp = u%data%good_user u%var%val_temp = u%var%saved u%var%good_temp = u%var%good_user !!! call dmerit_calc ('normal') ! to initialize other vars if needed ! Zero steerings do i = lbound(u%hsteer_kick%v, 1), ubound(u%hsteer_kick%v, 1) if (u%hsteer_kick%v(i)%exists) then u%hsteer_kick%v(i)%model = 0 endif if (u%vsteer_kick%v(i)%exists) then u%vsteer_kick%v(i)%model = 0 endif enddo !--------------------------------------------------------------- ! Fit the phase reference data... print *, 'Fitting the base_model to the reference phase data...' u%data%meas = u%data%ref u%cbar%m12%d%good_user = .false. u%skew_quad_k1%v%good_user = .false. u%quad_k1%v%good_user = .true. u%var%saved = u%var%saved_ref call opt_vars_set (opt_quad$, .true.) logic%opt_base = plot_model$ ! For merit calc, etc. call set_data_useit_opt (u%data) ! Run the optimizer !!! call dmerit_calc ('reinit') call run_optimizer ! plot the results call set_plot (graph%top1, u%phase, .true.) call set_plot (graph%bottom1, u%phase, .true.) call plot_data_set (graph%top1, plot_meas$) call plot_data_set (graph%bottom1, plot_meas$) call baseline_set (plot_model$, set_to$, graph%top1) call baseline_set (plot_design$, set_to$, graph%bottom1) call plotdo ('X', graph, .false., u) ! Emergency stop? if (.not. logic%auto_measurement) then call get_input_string ('"exit" for emergency exit. Anything else to continue...', ans) print * if (len_trim(ans) > 0 .and. index('exit', trim(ans)) == 1) return endif ! transfer the model to the base. call transfer ('MODEL', 'BASE_MODEL', 1.0_rp, u) ! reset u%data%meas = u%data%val_temp u%var%saved = u%var%val_temp !--------------------------------------------------------------- ! Fit the phase data - ref... u%data%good_user = .false. u%phase%x%d%good_user = u%phase%x%d%good_temp u%phase%y%d%good_user = u%phase%y%d%good_temp u%tune%x%d%good_user = u%tune%x%d%good_temp u%tune%y%d%good_user = u%tune%y%d%good_temp u%var%good_user = .false. call opt_vars_set (opt_all_vars$, .true.) logic%opt_base = plot_ref_and_model_and_base$ ! For merit calc, etc. call set_data_useit_opt (u%data) if (.not. k1_var%good_opt) then print *, 'ERROR: I CANNOT COMPUTE THE DERIVITIVE FOR: ', k1_var%name return endif ! Recompute dmerit vector for the varying quad. ! an initial call to dmerit_calc speeds things up if the other vars ! have not had their dmerit vectors calculated. k1_var%good_user = .true. k1_var%good_var = .true. k1_var%useit = .true. ! Run the optimizer. ! Turn off the limits since this can mess up the analysis. logic%limit_on = .false. !!! call dmerit_calc ('reinit') call run_optimizer !!! call dmerit_calc ('reinit') call run_optimizer ! Print results dk1 = k1_var%model - k1_var%base_model cal_meas = (u%ring%ele(0)%value(E_TOT$) / 1e9) * dk1 / cu_delta cal_old = (u%ring%ele(0)%value(E_TOT$) / 1e9) * k1_var%dvar_dcu d => u%phase%x%d sum2_meas_ref = sum2_calc (d%meas - d%ref, d%useit_opt) sum2_all = sum2_calc ((d%meas - d%ref) - (d%model - d%base_model), d%useit_opt) d => u%phase%y%d sum2_meas_ref = sum2_meas_ref + sum2_calc (d%meas - d%ref, d%useit_opt) sum2_all = sum2_all + sum2_calc ((d%meas - d%ref) - (d%model - d%base_model), d%useit_opt) dk1_err = sqrt(sum2_all/sum2_meas_ref) print * print '(1x, 5a, i4, a)', 'Calibrating: ', k1_var%name, '[', node_name, ',', k1_var%ix_db, ']' print *, 'Computed change in K1 needed to give Data:', dk1 print *, 'Computed change in CU needed to give Data:', nint(dk1 / k1_var%dvar_dcu) print '(a, f10.1)', ' Error in Fit (%):', 100 * dk1_err print * print *, ' Calibration in K1-GeV / CU:' print *, ' From the Measurement: ', cal_meas print *, ' From the Calibration File: ', cal_old print *, ' Ratio (Meas / Cal_File): ', cal_meas / cal_old !write (iu_raw, '(a, f10.1)') ' Error in Fit (%):', 100 * dk1_err !write (iu_raw, *) !write (iu_raw, *) ' Calibration in K1-GeV / CU:' !write (iu_raw, *) ' From the Measurement: ', cal_meas !write (iu_raw, *) ' From the Calibration File: ', cal_old !write (iu_raw, *) ' Ratio (Meas / Cal_File): ', cal_meas / cal_old ! plot the results call set_plot (graph%top1, u%phase, .true.) call set_plot (graph%bottom1, u%phase, .true.) call baseline_set (plot_ref_and_model_and_base$, set_to$, graph%top1) call baseline_set (plot_ref$, set_to$, graph%bottom1) call plotdo ('X', graph, .false., u) if (.not. logic%auto_measurement) then call get_input_string ('"exit" for emergency exit. Anything else to continue...', ans) print * if (len_trim(ans) > 0 .and. index('exit', trim(ans)) == 1) return endif ! Reset u%var%good_user = u%var%good_temp !--------------------------------------------------------------- ! Fit the quad kicks... x_kick_var => u%x_kick_quad%v(ix_quad) y_kick_var => u%y_kick_quad%v(ix_quad) u%data%good_user = .false. u%orbit%x%d%good_user = u%orbit%x%d%good_temp u%orbit%y%d%good_user = u%orbit%y%d%good_temp call transfer ('MODEL', 'BASE_MODEL', 1.0_rp, u) logic%opt_base = plot_ref_and_model_and_base$ ! For merit calc. call set_data_useit_opt (u%data) ! Recompute dmerit vector for the varying quad kicks. u%var%good_user = .false. x_kick_var%good_user = .true. y_kick_var%good_user = .true. do i = lbound(u%hsteer_kick%v, 1), ubound(u%hsteer_kick%v, 1) u%hsteer_kick%v(i)%model = 0 u%vsteer_kick%v(i)%model = 0 enddo call set_var_useit (u) !!! call dmerit_calc ('reinit') call run_optimizer ! u%data%fit = u%data%model - u%data%base_model ! Plot call set_plot (graph%top1, u%orbit, .true.) call set_plot (graph%bottom1, u%orbit, .true.) call plot_data_set (graph%top1, plot_meas$) call plot_data_set (graph%bottom1, plot_meas$) call baseline_set (plot_ref_and_model_and_base$, set_to$, graph%top1) call baseline_set (plot_ref$, set_to$, graph%bottom1) call setup_data_plot (graph%bottom1, u) call scale_data (0.0_rp, graph%bottom1%plot1, .false., graph%bottom1%plot2) call qp_calc_axis_places (graph%bottom1%plot1%y_axis) call qp_calc_axis_places (graph%bottom1%plot2%y_axis) graph%top1%plot1%y_axis = graph%bottom1%plot1%y_axis graph%top1%plot2%y_axis = graph%bottom1%plot2%y_axis call plotdo ('X', graph, .false., u) ! Calculate offset l_quad = u%ring%ele(ix_ele_quad)%value(l$) orb_x => u%orbit%x%d(ix_bpm) orb_y => u%orbit%y%d(ix_bpm) dx = x_kick_var%model / (l_quad * dk1) + orb_x%ref dy = -y_kick_var%model / (l_quad * dk1) + orb_y%ref orb_x_orig = (orb_x%meas * (cu_orig - cu1) + orb_x%ref * (cu2 - cu_orig)) / (cu2 - cu1) orb_y_orig = (orb_y%meas * (cu_orig - cu1) + orb_y%ref * (cu2 - cu_orig)) / (cu2 - cu1) bpm_offset => u%bpm_offset(ix_bpm) d => u%orbit%x%d rms_x = rms_calc (d%meas - d%ref - d%model + d%base_model, d%useit_opt) d => u%orbit%y%d rms_y = rms_calc (d%meas - d%ref - d%model + d%base_model, d%useit_opt) print * print '(a, 2f10.3)', ' X Y' print '(a, 2f10.3)', 'kick (urad): ', 1e6 * x_kick_var%model, 1e6 * y_kick_var%model print '(a, 2f10.3)', 'orbit_dat (mm) ', 1e3 * orb_x%meas, 1e3 * orb_y%meas print '(a, 2f10.3)', 'orbit_ref (mm) ', 1e3 * orb_x%ref, 1e3 * orb_y%ref print '(a, 2f10.3)', 'q_center-q_center_old (mm) ', 1e3 * dx, 1e3 * dy print '(a, 2f10.3)', 'q_center-orbit_orig (mm) ', 1e3 * (dx - orb_x_orig), 1e3 * (dy - orb_y_orig) print '(a, 2f10.3)', 'q_center-bpm_center (old mm):', 1e3 * bpm_offset%x, 1e3 * bpm_offset%y print '(a, 2f10.3)', 'q_center-bpm_center (new mm):', 1e3 * (bpm_offset%x + dx), 1e3 * (bpm_offset%y + dy) print '(a, 2f10.3)', 'RMS: ', 1e3 * rms_x, 1e3 * rms_y !write (iu_raw, *) !write (iu_raw, '(a, 2f10.3)') ' X Y' !write (iu_raw, '(a, 2f10.3)') 'kick (urad): ', 1e6 * x_kick_var%model, 1e6 * y_kick_var%model !write (iu_raw, '(a, 2f10.3)') 'orbit_dat (mm) ', 1e3 * orb_x%meas, 1e3 * orb_y%meas !write (iu_raw, '(a, 2f10.3)') 'orbit_ref (mm) ', 1e3 * orb_x%ref, 1e3 * orb_y%ref !write (iu_raw, '(a, 2f10.3)') 'q_center-q_center_old (mm) ', 1e3 * dx, 1e3 * dy !write (iu_raw, '(a, 2f10.3)') 'q_center-orbit_orig (mm) ', 1e3 * (dx - orb_x_orig), 1e3 * (dy - orb_y_orig) !write (iu_raw, '(a, 2f10.3)') 'q_center-bpm_center (old mm):', 1e3 * bpm_offset%x, 1e3 * bpm_offset%y !write (iu_raw, '(a, 2f10.3)') 'q_center-bpm_center (new mm):', 1e3 * (bpm_offset%x + dx), 1e3 * (bpm_offset%y + dy) !write (iu_raw, '(a, 2f10.3)') 'RMS: ', 1e3 * rms_x, 1e3 * rms_y ! Pause if (.not. logic%auto_measurement) then call get_input_string ('"exit" for emergency exit. Anything else to continue...', ans) print * if (len_trim(ans) > 0 .and. index('exit', trim(ans)) == 1) return endif ! Reset u%var%good_user = u%var%good_temp u%data%good_user = u%data%good_temp call set_data_useit_opt (u%data) x_kick_var%model = 0 y_kick_var%model = 0 k1_var%model = k1_var%base_model logic%limit_on = .true. call var_bookkeeper (x_kick_var, u%ring, u%orb) call var_bookkeeper (y_kick_var, u%ring, u%orb) call var_bookkeeper (k1_var, u%ring, u%orb) !------------------------------------------------------ ! Calculate bump to displace orbit to the quad center. ! First find steerings to use in the bump... call opt_vars_set (opt_steering$, .true.) u%hsteer_kick%v%weight = u%hsteer_kick%v%weight / 1000 u%vsteer_kick%v%weight = u%vsteer_kick%v%weight / 1000 logic%opt_base = plot_model_and_base$ ! For merit calc. ! Set the reference orbit. w_factor = 100 u%orbit%x%d%meas = 0 u%orbit%y%d%meas = 0 orb_x%meas = orb_x_orig - dx orb_y%meas = orb_y_orig - dy orb_x%weight = w_factor * orb_x%weight orb_y%weight = w_factor * orb_y%weight orb_x%good_user = .true. orb_y%good_user = .true. u%energy_data%d1%d(1)%good_user = .false. call transfer ('MODEL', 'BASE_MODEL', 1.0_rp, u) call set_data_useit_opt (u%data) ! Find the steering values for the bump. ! If the bump runs into a limit then need to use different steerings u%var%good_user = .false. u%var%hit_limit = .false. s_quad = u%ring%ele(ix_ele_quad)%s do call choose_close_steerings (u%hsteer_kick, s_quad, u) call choose_close_steerings (u%vsteer_kick, s_quad, u) call set_var_useit (u) call run_optimizer(hit_limit) if (.not. hit_limit) exit enddo !! call showit ('TOP10', u) ! Reset u%var%good_user = u%var%good_temp call set_var_useit(u) u%data%meas = u%data%val_temp orb_x%weight = orb_x%weight / w_factor orb_y%weight = orb_y%weight / w_factor u%hsteer_kick%v%weight = u%hsteer_kick%v%weight * 1000 u%vsteer_kick%v%weight = u%vsteer_kick%v%weight * 1000 ! Plot logic%opt_base_locked = .false. ! Unlock for normal running. call plot_data_set (graph%top1, plot_model$) call baseline_set (plot_base$, set_to$, graph%top1) call setup_data_plot (graph%top1, u) call scale_data (0.0_rp, graph%top1%plot1, .false., graph%top1%plot2) call plotdo ('X', graph, .false., u) !---------------------------------------------------------------------- ! Load bump if needed loaded_steerings = .false. if (ix_loop > 0 .and. ix_loop < logic%bpm_calib_max_cycles) then if (abs(dx-orb_x_orig) > logic%bpm_calib_max_dorbit .or. & abs(dy-orb_y_orig) > logic%bpm_calib_max_dorbit) then call load_steerings (u, 1.0_rp) loaded_steerings = .true. if (present(finished)) finished = .false. endif endif !---------------------------------------------------------------------------------- ! Save u%new_bpm_offset%ix_data = ix_bpm u%new_bpm_offset%good = .true. u%new_bpm_offset%x = bpm_offset%x + dx u%new_bpm_offset%y = bpm_offset%y + dy u%new_bpm_offset%sig_x = rms_x u%new_bpm_offset%sig_y = rms_y call date_and_time_stamp (date, .true.) u%new_bpm_offset%comment = date write (u%new_bpm_offset%comment, '(2a, i5.5, a, i5.5)') trim(u%new_bpm_offset%comment), & ' p:', u%phase%ix_ref, ' p:', u%phase%ix_meas print *, 'Note: Use SAVE BPM_CALIBRATION to save this calibration' ! If steerings have been loaded then do not write to detcal.dat iu = lunget() if (.not. loaded_steerings) then open (iu, file = 'detcal.dat', access = 'append') write (iu, '(i4, 3p, 4(a, f8.2), 2a)') ix_bpm, ',', u%new_bpm_offset%y, ',', u%new_bpm_offset%x, ',', & u%new_bpm_offset%sig_y, ',', u%new_bpm_offset%sig_x, ', ', trim(u%new_bpm_offset%comment) close (iu) endif ! Write to calibrate_bpm.dat.yyyy-mm-dd file_name = 'calibrate_bpm.dat' inquire (file = file_name, exist = does_exist) open (iu, file = file_name, access = 'append', recl = 160) if (.not. does_exist) then write (iu, '(a)') ' Phase_File Offset Old_Off Delta_Off dOrb ' write (iu, '(a)') 'Ix_bpm #1 #2 x y x y x y x y' endif write (iu, '(i6, 2i7, 3p, 4(2x, 2f6.2), 2x, a)') ix_bpm, u%phase%ix_ref, u%phase%ix_meas, & (bpm_offset%x + dx), (bpm_offset%y + dy), & bpm_offset%x, bpm_offset%y, dx, dy, (dx - orb_x_orig), (dy - orb_y_orig), trim(date) close (iu) !---------------------------------------------------------------------- contains function sum2_calc (data, useit_opt) result (sum2) real(rp) data(:), sum2 logical useit_opt(:) ! sum2 = sum(data**2, mask = useit_opt) - sum(data, mask = useit_opt)**2 / count(useit_opt) end function !---------------------------------------------------------------------- ! contains function rms_calc (data, useit_opt) result (rms) real(rp) data(:), rms integer n logical useit_opt(:) ! n = count(useit_opt) rms = sqrt(sum(data**2, mask = useit_opt) / n - (sum(data, mask = useit_opt) / n)**2) end function end subroutine