subroutine calibrate_skewquad (ix_var, how, cu_in, u, graph, err_flag)

use cesrv_struct
use cesrv_interface

implicit none

type (var_struct), pointer :: q_var
type (universe_struct), target :: u
type (ele_struct), pointer :: q_ele
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, dk1, dcbar_max
real(rp) dQ_data, dQ_ref, dQ_model

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

q_var => u%var(ix_var)
q_ele => u%ring%ele(u%var(ix_var)%ix_ele)
if (q_ele%key == overlay$) then
  q_ele => pointer_to_slave(q_ele, 1)
endif

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

if (.not. any(q_var%db_node_name == ['CSR SQEWQUAD'])) then
  print *, 'ERROR IN CALIBRATE_SKEWQUAD: THIS IS NOT A SKEW QUADRUPOLE: ', q_var%db_node_name
  return
endif

if (q_var%dvar_dcu == 0) then
  print *, 'ERROR IN CALIBRATE_SKEWQUAD: QUAD CALIBRATION FACTOR IS ZERO!', q_var%db_node_name
  return
endif


!---------------------------------------------------------------
! Recalibration using old data.

if (how == 'RECAL') then

  cu1 = q_var%cu_saved_ref
  cu2 = q_var%cu_saved
  cu_delta_actual = cu2 - cu1
  call inverse_cross_corr (q_var%db_node_name, cu1, cu2, cu_delta, q_var%ix_db, q_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 (q_var%db_node_name, q_var%ix_db, q_var%ix_db, cu_orig)

  ! Find delta. Default is 0.1 peak cbar11 wave.

  dcbar_max = 0.1

  if (cu_in(1) == 0 .and. cu_in(2) == 0) then
    dk1 = dcbar_max * abs(sin(u%ring%a%tune - u%ring%b%tune)) / &
                                  (sqrt(q_ele%a%beta * q_ele%b%beta) * q_ele%value(l$))
    cu_delta_actual = abs(dk1 / q_var%dvar_dcu) 
    cu_delta_actual = min (cu_delta_actual, q_var%cu_high_lim - q_var%cu_low_lim - 50)

    ! calc what max and min should be.
    ! try to make max and min not near 0

    cu1 = (q_var%cu_high_lim + q_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) < 100) then
      if (q_var%cu_low_lim < -115) then
        cu1 = -100
        cu2 = cu1 + cu_delta_actual
      endif
    elseif (abs(cu2) < 100) then
      if (q_var%cu_high_lim > 115) then
        cu2 = 100
        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 (q_var%db_node_name, cu1, cu2, cu_delta, &
                                              q_var%ix_db, q_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) > q_var%cu_high_lim .or. min(cu1,cu2) < q_var%cu_low_lim) then
    print *, 'ERROR IN CALIBRATE_SKEWQUAD: A SET POINT IS OUTSIDE A LIMIT'
    print *, '      LIMITS ARE:', q_var%cu_high_lim, q_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 (q_var%db_node_name, q_var%ix_db, cu1, &
                          phase_data$, (/ix_var/), 1, ref_file$, u, graph, err_flag)
    if (err_flag) then
      if (.not. logic%debug) call vxputn (q_var%db_node_name, q_var%ix_db, q_var%ix_db, cu_orig)
      err_string = '*** Error: Set to first Set Point failed'
      goto 8000
    endif


    call set_and_meas (q_var%db_node_name, q_var%ix_db, cu2, &
                         phase_data$, (/ix_var/), 1, data_file$, u, graph, err_flag)
    if (err_flag) then
      if (.not. logic%debug) call vxputn (q_var%db_node_name, q_var%ix_db, q_var%ix_db, cu_orig)
      err_string = '*** Error: Set to second Set Point failed'
      goto 8000
    endif

    if (.not. logic%debug) call vxputn (q_var%db_node_name, q_var%ix_db, q_var%ix_db, cu_orig)

  endif

!---------------------------------------------------------------
! error

else
  print *, 'INTERNAL ERROR IN CALIBRATE_SKEWQUAD: BAD "HOW" ARGUMENT.'
  print *, '         GET HELP.'
  return

endif

! Check tune spread

dQ_data = u%tune%x%d(1)%meas - u%tune%y%d(1)%meas
dQ_ref = u%tune%x%d(1)%ref - u%tune%y%d(1)%ref
dQ_model = u%tune%x%d(1)%model - u%tune%y%d(1)%model
if (abs(dQ_data-dQ_model) > 0.01 * dQ_model .and. abs(dQ_ref-dQ_model) > 0.01 * dQ_model) then
  print *, 'WARNING! TUNE SPLIT IN MODEL OVER 1% DIFFERENT THAN ACTUAL TUNE SPLIT IN EITHER DATA SETS!'
  print *, 'REMEMBER! CALCULATED CALIBRATION IS DIRECTLY PROPORTIONAL TO dQ!'
endif

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

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

name = q_var%db_node_name
u%phase%x%d%good_user = .false.
u%phase%y%d%good_user = .false.
call set_plot (graph%top1, u%cbar)
call set_plot (graph%bottom1, u%cbar)

call set_data_useit_opt (u%data)

if (.not. q_var%good_opt) then
  print *, 'ERROR: I CANNOT COMPUTE THE DERIVITIVE FOR: ', q_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
q_var%good_user = .true.
q_var%good_var = .true.
q_var%useit = .true.
q_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 = q_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(q_var%ix_ele)%value(E_TOT$) / 1e9) * dvar / cu_delta
cal_old  = (u%ring%ele(q_var%ix_ele)%value(E_TOT$) / 1e9) * q_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: ', q_var%name, '[', &
                              q_var%db_node_name, ',', q_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
q_var%model = q_var%model + dvar
call var_bookkeeper (q_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(q_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 = 'skew_quad_' // 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)') q_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)') q_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)

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

logic%ok_to_read_dmerit_file = old_ok_to_read_dmerit

end subroutine