subroutine calibrate_steering (ix_var, how, cu_in, u, graph, err_flag) use cesrv_struct use cesrv_interface implicit none type (var_struct), pointer :: this_var type (universe_struct), target :: u type (ele_struct) ele type (graph_struct) graph integer cu_delta, cu_orig, iv, i, cu1, cu2, cu_in(2) integer iu, n_dat, ix_var, cu_delta_actual real(rp) sum_sd, sum_s2, dvar, cal_meas, cal_old, slope, beta real(rp) dvar_err, dmeas, factor, sum_resid, sig_resid character(*) how character date*20, err_string*40 character(140) cal_file logical make_meas, err_flag, is_there, old_ok_to_read_dmerit ! Init this_var => u%var(ix_var) if (.not. any(this_var%db_node_name == ['CSR HORZ CUR', 'CSR VERT CUR', & 'CSR HBND CUR', 'CSR HSP VOLT', 'UND VERT CUR'])) then print *, 'ERROR IN CALIBRATE_STEERING: THIS IS NOT A STEERING: ', this_var%db_node_name return endif !--------------------------------------------------------------- ! Recalibration using old data. if (how == 'RECAL') then cu1 = this_var%cu_saved_ref cu2 = this_var%cu_saved cu_delta_actual = cu2 - cu1 call inverse_cross_corr (this_var%db_node_name, cu1, cu2, cu_delta, & this_var%ix_db, this_var%ix_db) print *, 'Steering Range used CU(ref)/CU(data): ', cu1, cu2 print *, 'Actual Delta CU: ', cu_delta_actual print *, 'Effective Delta CU: (zero crossing corrected)', cu_delta !--------------------------------------------------------------- ! This section if we need to take a measurement elseif (how == 'NEW_CAL') then call vxgetn (this_var%db_node_name, this_var%ix_db, this_var%ix_db, cu_orig) ! Find delta. Default is to have 4 mm max amplitude at beta = 20 m ignoring ! the resonance denominator, rounded off to the nearest 100 CU if (cu_in(1) == 0 .and. cu_in(2) == 0) then call twiss_at_element (u%ring, this_var%ix_ele, average = ele) if (this_var%ix_attrib == hkick$) then beta = ele%a%beta factor = abs(sin(u%ring%a%tune/2)) else beta = ele%b%beta factor = abs(sin(u%ring%b%tune/2)) endif cu_delta_actual = 2 * 0.004 * factor / (sqrt(beta * 20) * & abs(this_var%dvar_dcu)) cu_delta_actual = 100 * nint (cu_delta_actual / 100.0_rp) cu_delta_actual = max(100, cu_delta_actual) cu_delta_actual = min (cu_delta_actual, & int(0.9 * this_var%cu_high_lim - this_var%cu_low_lim)) ! calc what max and min should be. ! try to make max and min not near 0 cu1 = (this_var%cu_high_lim + this_var%cu_low_lim - cu_delta_actual) / 2 cu2 = cu1 + cu_delta_actual if (cu_orig < cu1) then cu1 = cu_orig cu2 = cu1 + cu_delta_actual elseif (cu_orig > cu2) then cu2 = cu_orig cu1 = cu2 - cu_delta_actual endif if (abs(cu1) < 200) then if (this_var%cu_low_lim < -225) then cu1 = -200 cu2 = cu1 + cu_delta_actual endif elseif (abs(cu2) < 200) then if (this_var%cu_high_lim > 225) then cu2 = 200 cu1 = cu2 - cu_delta_actual endif endif if (this_var%db_node_name == 'CSR HBND CUR' .and. cu1 < 1000) then cu1 = 1000 cu2 = min(cu1 + cu_delta_actual, this_var%cu_high_lim) endif else cu1 = cu_in(1) cu2 = cu_in(2) endif ! Type out endpoint info cu_delta_actual = cu2 - cu1 call inverse_cross_corr (this_var%db_node_name, cu1, cu2, cu_delta, & this_var%ix_db, this_var%ix_db) print *, 'Range used CU(ref)/CU(data): ', cu1, cu2 print *, 'Actual Delta CU: ', cu_delta_actual print *, 'Effective Delta CU (zero crossing corrected):', cu_delta if (max(cu1,cu2) > this_var%cu_high_lim .or. & min(cu1,cu2) < this_var%cu_low_lim) then print *, 'ERROR IN CALIBRATE_STEERING: A SET POINT IS OUTSIDE A LIMIT' print *, ' LIMITS ARE:', this_var%cu_high_lim, this_var%cu_low_lim err_flag = .true. err_string = '*** Error: Set point outside limit' goto 8000 ! write error message to cal file endif ! measure orbit 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 orbit 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 call set_and_meas (this_var%db_node_name, this_var%ix_db, cu1, & orbit_data$, (/ix_var/), 1, ref_file$, u, graph, err_flag) if (err_flag) then err_string = '*** Error: Set to first Set Point failed' goto 8000 endif call set_and_meas (this_var%db_node_name, this_var%ix_db, cu2, & orbit_data$, (/ix_var/), 1, data_file$, u, graph, err_flag) if (err_flag) then err_string = '*** Error: Set to second Set Point failed' goto 8000 endif if (.not. logic%debug) call vxputn (this_var%db_node_name, this_var%ix_db, this_var%ix_db, cu_orig) endif !--------------------------------------------------------------- ! error else print *, 'INTERNAL ERROR IN CALIBRATE_STEERING: BAD "HOW" ARGUMENT.' print *, ' GET HELP.' return endif !--------------------------------------------------------------- ! Analyze data. if (.not. this_var%good_opt) then print *, 'ERROR: I CANNOT COMPUTE THE DERIVITIVE FOR: ', this_var%name print *, ' SINCE THIS VARIABLE IS NOT IN THE OPIMIZATION LIST' print *, ' YOU PROBABLY NEED TO USE THE COMMAND: "OPT STEERING"' return endif ! recompute dmerit vector for the varying steering. ! an initial call to dmerit_calc speeds things up if the other vars ! have not had their dmerit vectors calculated. call dmerit_calc ('normal') ! to initialize other vars if needed old_ok_to_read_dmerit = logic%ok_to_read_dmerit_file ! save this_var%good_user = .true. this_var%good_var = .true. this_var%useit = .true. this_var%dm_dv_computed = .false. logic%ok_to_read_dmerit_file = .false. call dmerit_calc ('normal') ! sum_sd = 0 sum_s2 = 0 n_dat = 0 iv = this_var%ix_dvar do i = 1, size(u%data) if (.not. associated(u%data(i)%d1)) cycle if (u%data(i)%d1%d2%type /= orbit_data$) cycle if (.not. u%data(i)%useit_opt) cycle slope = u%dm_dv(u%data(i)%ix_dmeas, iv) dmeas = (u%data(i)%meas - u%data(i)%ref) sum_sd = sum_sd + slope * dmeas sum_s2 = sum_s2 + slope**2 n_dat = n_dat + 1 enddo if (n_dat == 0) then print *, 'EGADS! No GOOD DATA! Aborting analysis!' return endif factor = 1.0e6 * (u%ring%ele(this_var%ix_ele)%value(E_TOT$) / 1e9) if (this_var%ix_attrib == hkick$) factor = -factor dvar = (sum_sd / sum_s2) cal_meas = factor * dvar / cu_delta cal_old = factor * this_var%dvar_dcu sum_resid = 0 do i = 1, size(u%data) if (.not. associated(u%data(i)%d1)) cycle if (u%data(i)%d1%d2%type /= orbit_data$) cycle if (.not. u%data(i)%useit_opt) cycle slope = u%dm_dv(u%data(i)%ix_dmeas, iv) dmeas = (u%data(i)%meas - u%data(i)%ref) sum_resid = sum_resid + (slope * dvar - dmeas)**2 enddo sig_resid = sqrt(sum_resid * sum_s2 / (n_dat * sum_sd**2)) dvar_err = 100 * sig_resid print * print '(1x, 5a, i4, a)', 'Calibrating: ', this_var%name, '[', & this_var%db_node_name, ',', this_var%ix_db, ']' print *, 'Computed Change in Kick needed to give Data (mrad):', 1000*dvar print '(a, f10.1)', ' Error in Fit (%):', dvar_err print * print *, ' Calibration in mrad-GeV / 1000 CU:' print *, ' From the Measurement: ', cal_meas print *, ' From the Calibration File: ', cal_old print *, ' Ratio (Meas / Cal_File):', cal_meas / cal_old err_flag = .false. ! set for plotting u%data%old = u%data%model this_var%model = this_var%model + dvar call var_bookkeeper (this_var, u%ring, u%orb) call ring_calc (u) u%data%fit = u%data%model - u%data%old call set_plot (graph%top1, u%orbit) call set_plot (graph%bottom1, u%orbit) call plot_data_set (graph%top1, plot_meas$) call plot_data_set (graph%bottom1, plot_meas$) call baseline_set (plot_ref_and_fit$, set_to$, graph%top1) call baseline_set (plot_ref$, set_to$, graph%bottom1) u%main_title1 = 'Top Plot: Orbit_Difference - Model_Fit' write (u%main_title2, '(3a, i8)') 'Bottom Plot: Difference Orbit for ', & trim(u%ring%ele(this_var%ix_ele)%name), ' With delta CU =', cu_delta call plotdo ('X', graph, .false., u) ! write results to the data file 8000 continue iu = lunget() call date_and_time_stamp (date) cal_file = 'steering_' // date(1:11) // '.cal_meas' inquire (file = cal_file, exist = is_there) open (iu, file = cal_file, status = 'unknown', access = 'append') if (.not. is_there) write (iu, '(3(/, 2x, a))') & ' Fit CORSTR %Fit', & 'Steer CU1 CU2 dCU Orb1 Orb2 Cal Cal Ratio Error Ix', & '----- ---- ---- ---- ----- ----- ----- ----- ------ ----- --' if (err_flag) then write (iu, '(a7, 3i7, 2x, a)') this_var%name, & cu1, cu2, cu_delta, trim(err_string) print *, 'Error recorded in: ', trim(cal_file) else write (iu, '(a7, 5i7, 3f8.3, f7.1, i4)') this_var%name, & cu1, cu2, cu_delta, u%orbit%ix_ref, u%orbit%ix_meas, & cal_meas, cal_old, cal_meas/cal_old, dvar_err print *, 'Calibration recorded in: ', trim(cal_file), ix_var endif close (iu) ! reset call set_data_useit_opt (u%data) this_var%model = this_var%model - dvar logic%ok_to_read_dmerit_file = old_ok_to_read_dmerit end subroutine