subroutine calibrate_quad (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 (graph_struct) graph

integer cu_delta, cu_orig, iv, i, cu1, cu2, cu_in(2)
integer iu, n_dat, ix_var, type, cu_delta_actual

real(rp) sum_d2, sum_sd, sum_s2, dvar, cal_meas, cal_old, slope
real(rp) dvar_err, dmeas, sum_resid, sig_resid

character(*) how
character(140) cal_file
character date*20, name*12, err_string*40
character(7) ratio_str

logical make_meas, err_flag, is_there, old_ok_to_read_dmerit

! Init

this_var => u%var(ix_var)

print '(1x, 5a, i4, a)', 'Calibrating: ', this_var%name, '[', &
                                this_var%db_node_name, ',', this_var%ix_db, ']'

if (.not. any(this_var%db_node_name == ['CSR QUAD CUR''CSR QUAD CUR'])) then
  print *, 'ERROR IN CALIBRATE_QUAD: THIS IS NOT A QUADRUPOLE: ', 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 *, '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 2000 CU.

  if (cu_in(1) == 0 .and. cu_in(2) == 0) then
    cu_delta_actual = 2000
    cu_delta_actual = min (cu_delta_actual, this_var%cu_high_lim - this_var%cu_low_lim - 50)

    ! 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

  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_QUAD: 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 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 a 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
    call set_and_meas (this_var%db_node_name, this_var%ix_db, cu1, &
            phase_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, &
             phase_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_QUAD: BAD "HOW" ARGUMENT.'
  print *, '         GET HELP.'
  return

endif

!---------------------------------------------------------------
! Analyze data.

u%data%good_temp = u%data%good_user  ! save current

name = this_var%db_node_name
u%cbar%m11%d%good_user = .false.
u%cbar%m12%d%good_user = .false.
u%cbar%m21%d%good_user = .false.
u%cbar%m22%d%good_user = .false.
call set_plot (graph%top1, u%phase)
call set_plot (graph%bottom1, u%phase)

call set_data_useit_opt (u%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 COMMOND: "OPT QUADRUPOLE"'
  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.

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
sum_d2 = 0
n_dat = 0

iv = this_var%ix_dvar
do i = 1, size(u%data)
  if (.not. associated(u%data(i)%d1)) cycle
  type = u%data(i)%d1%d2%type
  if (type /= phase_data$ .and. type /= cbar_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
  sum_d2 = sum_d2 + dmeas**2
  n_dat = n_dat + 1
enddo

if (n_dat == 0) then
  print *, 'EGADS! No GOOD DATA! Aborting analysis!'
  return
endif

dvar = (sum_sd / sum_s2)
cal_meas = (u%ring%ele(this_var%ix_ele)%value(E_TOT$) / 1e9) * dvar / cu_delta
cal_old  = (u%ring%ele(this_var%ix_ele)%value(E_TOT$) / 1e9) * this_var%dvar_dcu

sum_resid = 0
do i = 1, size(u%data)
  if (.not. associated(u%data(i)%d1)) cycle
  type = u%data(i)%d1%d2%type
  if (type /= phase_data$ .and. type /= cbar_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) * slope)**2
enddo

sig_resid = sqrt(sum_resid) / sum_s2
dvar_err = 100 * sig_resid / abs(dvar)

if (cal_old == 0) then
  ratio_str = '  -----'
else
  write (ratio_str, '(f7.2)') cal_meas / cal_old
endif

print *
print '(1x, 5a, i4, a)', 'Calibrating: ', this_var%name, '[', &
                              this_var%db_node_name, ',', this_var%ix_db, ']'
print *, 'Computed Change in K1 needed to give Data:', dvar
print '(a, f10.1)', ' Error in Fit (%):', dvar_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): ', ratio_str

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 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: Data - Ref - Model_Fit'
write (u%main_title2, '(3a, i8)') 'Bottom Plot: Data - Ref 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 = 'quadrupole_' // 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(/, a))') &
'                                Phase  Phase                        Meas/ %Fit', &
'     Quad    CU1    CU2    dCU    #1     #2     Fit Cal    Old Cal   Old  Error', &
'---------- -----  -----  -----  -----  -----   --------   --------  ----- -----'

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, '(a10, i6, 2i7, 2(2x, i5), 1p, 2e11.2, 0p, a, f6.1)') this_var%name, &
                cu1, cu2, cu_delta, u%phase%ix_ref, u%phase%ix_meas, &
                cal_meas, cal_old, ratio_str, dvar_err
  print *, 'Calibration recorded in: ', trim(cal_file)
endif

close (iu)

! reset

u%data%good_user = u%data%good_temp
call set_data_useit_opt (u%data)

this_var%model = this_var%model - dvar
call var_bookkeeper (this_var, u%ring, u%orb)
call ring_calc (u)

logic%ok_to_read_dmerit_file = old_ok_to_read_dmerit

end subroutine