!+
! Subroutine read_phase (data_or_ref, num_in, u, graph, err_flag, reanal)
!
! Subroutine to read in phase/coupling data from an input file.
!-

subroutine read_phase (data_or_ref, num_in, u, graph, err_flag, reanal)

use cesrv_struct
use cesrv_interface
use phase_mod
use sim_utils
use cesr_basic_mod
use cesr_read_data_mod

implicit none

type my_beta_struct
  real(rp) beta_a, beta_b
  logical ok_beta_a, ok_beta_b
end type

type (my_beta_struct) beta_meas(0:n_det_maxx)
type (universe_struct), target :: u
type (phase_cbar_data_struct) pc_(0:n_det_maxx)
type (detector_struct) orbit_(0:n_det_maxx)
type (graph_struct) graph
type (cesr_phase_params_struct) p
type (cesr_det_plane_struct) :: horz(0:n_det_maxx), vert(0:n_det_maxx)
type (cesr_freq_struct) :: freq
type (var_struct), pointer :: vp
type (cesr_det_dc_position_struct) dc(0:n_det_maxx)
type (butns_struct) butns
type (raw_det_struct), pointer :: raw_det_x(:), raw_det_y(:)
type (cesr_data_params_struct), pointer :: data_params
!type but_value_struct
!  real(rp) value(4)
!end type

type (but_value_struct), save :: gain(0:120)

! for applying BPM tilts:
real(rp) :: tilt(0:n_det_maxx), cbar_mat_a(2,2), cbar_mat_b(2,2)
type (ele_struct) :: ele

namelist / phase_cbar_data / pc_
namelist / raworbit / orbit_
namelist / beta_data / beta_meas

real(rp) r_beta

integer  data_or_ref, i, j, ios, phase_num, num_in, det, iu
integer ix_99, ix_0, ix_ip, ie, raw(4, 120), det_type(120)

character(140) file_name
character(140) phase_number

logical :: reanal
logical err_flag, good, err

!-----------------------------------------------------------------------
! since we are reading in phase data we set the optimization to quad

logic%use_old_in_loading = .false.

if (logic%opt_vars == opt_steering$) then
  logic%opt_vars = opt_quad$
  call baseline_set (plot_design$, set_to$, graph%top1)
  call baseline_set (plot_design$, set_to$, graph%bottom1)
endif

! read a data file...
! first construct the file name

call fullfilename(&
     '$CESR_ONLINE/machine_data/mach_meas/phase/phase.number',phase_number)
call calc_file_number (phase_number, num_in, phase_num, err_flag)
if (err_flag) return
call form_file_name_with_number ('PHASE', phase_num, file_name, err_flag)
if (err_flag) return

! open the file and read the contents

err_flag = .true.
iu = lunget()
open (iu, file = file_name, status = 'old', action = 'read', iostat = ios)
if (ios /= 0) then       ! abort on open error
  print *, 'ERROR OPENING: ', trim(file_name)
  return
endif

call read_header (iu, data_or_ref, u%phase, u)

!

rewind (iu)

read (iu, nml = raworbit, iostat = ios)
if (ios /= 0) then
  print *, 'ERROR READING RAWORBIT'
  return
endif

pc_(:)%ok_x = .false.
pc_(:)%ok_y = .false.
read (iu, nml = phase_cbar_data, iostat = ios)

if (ios /= 0) then
  print *, 'ERROR READING PHASE_CBAR_DATA'
  return
endif

beta_meas%ok_beta_a = .false.
beta_meas%ok_beta_b = .false.
read (iu, nml = beta_data, iostat = ios)

close (iu)

! re-anal data?

if (reanal) then

  print *, 'RE-ANALIZING PHASE DATA.'

  p%lattice = logic%lattice
  p%species = u%ring%param%particle

  freq%rev = 390.14
  if (data_or_ref == data_file$) then
    freq%x%reflection = u%phase%data_params%horiz_reflection_shake
    freq%x%tune = u%tune%x%d(1)%meas * 390.1 / twopi
    freq%y%reflection = u%phase%data_params%vert_reflection_shake
    freq%y%tune = u%tune%y%d(1)%meas * 390.1 / twopi
    raw_det_x => u%raw_phase_x%det
    raw_det_y => u%raw_phase_y%det
  else
    freq%x%reflection = u%phase%ref_params%horiz_reflection_shake
    freq%x%tune = u%tune%x%d(1)%ref * 390.1 / twopi
    freq%y%reflection = u%phase%ref_params%vert_reflection_shake
    freq%y%tune = u%tune%y%d(1)%ref * 390.1 / twopi
    raw_det_x => u%raw_phase_x%det_ref
    raw_det_y => u%raw_phase_y%det_ref
  endif

  print *,' READ_PHASE: logic%gain_correction = ', logic%gain_correction
  call get_gain(gain)

  do det = 0, n_det_maxx
    horz(det)%but(:)%amp   = raw_det_x(det)%amp(:) * gain(det)%value(:)
    horz(det)%but(:)%phase = raw_det_x(det)%phase(:)
    horz(det)%ok           = pc_(det)%ok_x
    horz(det)%system_id    = raw_det_x(det)%system_id
    if (horz(det)%system_id == 0) horz(det)%but(:)%amp = 1000 * horz(det)%but(:)%amp
    call button_phase (horz(det), 0.0_rp, x_plane$)

    vert(det)%but(:)%amp   = raw_det_y(det)%amp(:) * gain(det)%value(:)
    vert(det)%but(:)%phase = raw_det_y(det)%phase(:)
    vert(det)%ok           = pc_(det)%ok_y
    vert(det)%system_id    = raw_det_y(det)%system_id
    if (vert(det)%system_id == 0) vert(det)%but(:)%amp = 1000 * vert(det)%but(:)%amp
    call button_phase (vert(det), 0.0_rp, y_plane$)
    dc(det)%signal = orbit_(det)%amp
    dc(det)%x = orbit_(det)%x_orb
    dc(det)%y = orbit_(det)%y_orb
  enddo

  call analyze_phase_data (horz, vert, freq, dc, p, u)

  do det = 0, n_det_maxx
    pc_(det)%x_phase  = horz(det)%phase_meas
    pc_(det)%x_cbar22 = horz(det)%cbar22
    pc_(det)%x_cbar12 = horz(det)%cbar12
    pc_(det)%y_phase  = vert(det)%phase_meas
    pc_(det)%y_cbar12 = vert(det)%cbar12
    pc_(det)%y_cbar11 = vert(det)%cbar11
    beta_meas(det)%beta_a    = horz(det)%beta_meas
    beta_meas(det)%beta_b    = vert(det)%beta_meas
    beta_meas(det)%ok_beta_a = horz(det)%ok
    beta_meas(det)%ok_beta_b = vert(det)%ok
  enddo

endif

! Convert orbit data

do i = 0, 120
  if (i == 100) cycle
  j = i
  if (j == 0) j = 100
  raw(:, j) = orbit_(i)%amp
  if (data_or_ref == data_file$) then
    det_type(j) = u%raw_phase_x%det(i)%system_id
  else
    det_type(j) = u%raw_phase_x%det_ref(i)%system_id
  endif
enddo

call raw_butns_to_position (raw, det_type, butns, err, &
    logic%nonlinear_calc, logic%offset_correction, logic%gain_correction)


!---------------------------------------------------------------------------
! Veto all quads that have a command less then 1000 from varying.

do i = lbound(u%quad_k1%v, 1), ubound(u%quad_k1%v, 1)
  vp => u%quad_k1%v(i)
  if (abs(vp%cu_saved) < 1000 .and. vp%good_var) then
    print *
    print *, 'NOTICE: QUAD COMMAND < 1000 CU FOR: ', vp%name
    print *, '        I AM AUTO VETOING THIS QUAD!'
    vp%good_var = .false.
  endif
enddo

!---------------------------------------------------------------------------
! apply tilts to Cbar, if necessary

if (logic%tilt_correction) then
   call get_tilt(tilt)
   write(*,*) "Applying tilts from ", trim(logic%tilt_correction_file)
    if (logic%tilt_correction) then

       do i = 0, n_det_maxx
          
          ie = u%orbit%x%d(i)%ix_ele
          r_beta = sqrt(u%ring%ele(ie)%b%beta / u%ring%ele(ie)%a%beta)

          ele = u%ring%ele(ie)

          ! rotate by -tilt(i); tilts are saved from fitted model in CESRv
          ! Note: have to repeat twice, to accomodate cbar12 from both 
          ! a-mode and b-mode measurements

          ele%c_mat(2,2) = pc_(i)%x_cbar22! * r_beta
          ele%c_mat(1,2) = pc_(i)%x_cbar12! * r_beta
          ele%c_mat(1,1) = pc_(i)%y_cbar11! / r_beta
          ele%c_mat(2,1) = 0.

          call rotate_cbar (ele, -tilt(i), cbar_mat_a, .true.)
          
          ele%c_mat(2,2) = pc_(i)%x_cbar22! * r_beta
          ele%c_mat(1,2) = pc_(i)%y_cbar12! / r_beta
          ele%c_mat(1,1) = pc_(i)%y_cbar11! / r_beta
          ele%c_mat(2,1) = 0.

          call rotate_cbar (ele, -tilt(i), cbar_mat_b, .true.)

          !write(*,'(a,i6,4es18.5)') "foo - ", i, pc_(i)%x_cbar12, cbar_mat_a(1,2), pc_(i)%y_cbar12, cbar_mat_b(1,2)
          pc_(i)%x_cbar22 = cbar_mat_a(2,2) 
          pc_(i)%x_cbar12 = cbar_mat_a(1,2) 
          pc_(i)%y_cbar11 = cbar_mat_b(1,1) 
          pc_(i)%y_cbar12 = cbar_mat_b(1,2) 
          
       enddo
           

    endif

endif


!---------------------------------------------------------------------------
! Read Phase data

select case (data_or_ref)
case (data_file$)

  u%phase%measured = .true.
  u%cbar%measured  = .true.
  u%cmat_a%measured = .true.
  u%cmat_b%measured = .true.
  u%beta%measured  = .true.
  u%orbit%measured = .true.

  u%phase%file_name = file_name
  u%cbar%file_name  = file_name
  u%cmat_a%file_name  = file_name
  u%cmat_b%file_name  = file_name
  u%beta%file_name  = file_name
  u%orbit%file_name = file_name

  u%phase%ix_meas = phase_num
  u%cbar%ix_meas  = phase_num
  u%cmat_a%ix_meas  = phase_num
  u%cmat_b%ix_meas  = phase_num
  u%beta%ix_meas  = phase_num

  u%cbar%m11%d(:)%good_dat = .false.
  u%cbar%m12%d(:)%good_dat = .false.
  u%cbar%m22%d(:)%good_dat = .false.

  u%cmat_a%m22%d(:)%good_dat = .false.
  u%cmat_a%m12%d(:)%good_dat = .false.
  u%cmat_b%m12%d(:)%good_dat = .false.
  u%cmat_b%m11%d(:)%good_dat = .false.

  u%orbit%x%d(:)%meas = butns%det%x_orb
  u%orbit%y%d(:)%meas = butns%det%y_orb
  u%orbit%x%d%good_dat = butns%det%ok
  u%orbit%y%d%good_dat = butns%det%ok

  u%phase%x%d(:)%good_dat = pc_%ok_x .and. pc_%ok_y
  u%phase%y%d(:)%good_dat = pc_%ok_x .and. pc_%ok_y
  u%phase%x%d(:)%meas = pc_%x_phase * twopi / 360
  u%phase%y%d(:)%meas = pc_%y_phase * twopi / 360

  u%cbar%m22%d(:)%meas = pc_%x_cbar22
  u%cbar%m12%d(:)%meas = pc_%x_cbar12
  u%cbar%m11%d(:)%meas = pc_%y_cbar11

  u%cbar%m11%d(:)%good_dat = pc_%ok_x .and. pc_%ok_y
  u%cbar%m12%d(:)%good_dat = pc_%ok_x .and. pc_%ok_y
  u%cbar%m22%d(:)%good_dat = pc_%ok_x .and. pc_%ok_y

  u%cmat_a%m22%d(:)%good_dat = pc_%ok_x .and. pc_%ok_y
  u%cmat_a%m12%d(:)%good_dat = pc_%ok_x .and. pc_%ok_y
  u%cmat_b%m11%d(:)%good_dat = pc_%ok_x .and. pc_%ok_y
  u%cmat_b%m12%d(:)%good_dat = pc_%ok_x .and. pc_%ok_y

  u%beta%x%d(:)%meas = beta_meas%beta_a
  u%beta%y%d(:)%meas = beta_meas%beta_b
  u%beta%x%d(:)%good_dat = beta_meas%ok_beta_a
  u%beta%y%d(:)%good_dat = beta_meas%ok_beta_b

  do i = 0, n_det_maxx

    ie = u%orbit%x%d(i)%ix_ele
    r_beta = sqrt(u%ring%ele(ie)%b%beta / u%ring%ele(ie)%a%beta)
    u%cmat_a%m22%d(i)%meas = pc_(i)%x_cbar22 * r_beta
    u%cmat_a%m12%d(i)%meas = pc_(i)%x_cbar12 * r_beta
    u%cmat_b%m11%d(i)%meas = pc_(i)%y_cbar11 / r_beta
    u%cmat_b%m12%d(i)%meas = pc_(i)%y_cbar12 / r_beta

  enddo

  !

  print *, 'Phase/Coupling/Orbit Data Read In.'

!!  print *, '    NOTE: I am automatically VETOING all Detector 5 Data.'
!!  u%phase%x%d(5)%good_dat = .false.
!!  u%phase%y%d(5)%good_dat = .false.
!!  u%cbar%m11%d(5)%good_dat = .false.
!!  u%cbar%m12%d(5)%good_dat = .false.
!!  u%cbar%m22%d(5)%good_dat = .false.

  ! Reverse phases for electrons

  if (u%ring%param%particle == electron$) then
    do i = 0, n_det_maxx
      u%phase%x%d(i)%meas = u%tune%x%d(1)%meas - u%phase%x%d(i)%meas
      u%phase%y%d(i)%meas = u%tune%y%d(1)%meas - u%phase%y%d(i)%meas
    enddo
  endif                    

  ! Set plots.

  select case (graph%top1%d2%type)
  case (phase_data$, cbar_data$, beta_data$, cmat_a_data$, cmat_b_data$) 
  case default
    call set_plot (graph%top1, u%phase)
    call plot_data_set (graph%top1, plot_meas$)
  end select

  select case (graph%bottom1%d2%type)
  case (phase_data$, cbar_data$, beta_data$, cmat_a_data$, cmat_b_data$) 
  case default
    call set_plot (graph%bottom1, u%cbar)
    call plot_data_set (graph%bottom1, plot_meas$)
  end select

  u%beta%p2%plot_data = plot_meas$
  u%orbit%p2%plot_data = plot_meas$
  u%cmat_a%p2%plot_data = plot_meas$
  u%cmat_b%p2%plot_data = plot_meas$

!--------------------------------------------------------------------
! Read phase ref

case (ref_file$) 

  u%phase%ix_ref = phase_num
  u%cbar%ix_ref  = phase_num
  u%cmat_a%ix_ref  = phase_num
  u%cmat_b%ix_ref  = phase_num
  u%beta%ix_ref  = phase_num

  u%cbar%m11%d(:)%good_ref = .false.
  u%cbar%m12%d(:)%good_ref = .false.
  u%cbar%m22%d(:)%good_ref = .false.

  u%cmat_a%m22%d(:)%good_ref = .false.
  u%cmat_a%m12%d(:)%good_ref = .false.
  u%cmat_b%m12%d(:)%good_ref = .false.
  u%cmat_b%m11%d(:)%good_ref = .false.

  u%orbit%x%d%good_ref = butns%det%ok
  u%orbit%y%d%good_ref = butns%det%ok
  u%orbit%x%d(:)%ref = butns%det%x_orb
  u%orbit%y%d(:)%ref = butns%det%y_orb

  u%phase%x%d(:)%good_ref = pc_%ok_x .and. pc_%ok_y
  u%phase%y%d(:)%good_ref = pc_%ok_x .and. pc_%ok_y
  u%phase%x%d(:)%ref = pc_%x_phase * twopi / 360
  u%phase%y%d(:)%ref = pc_%y_phase * twopi / 360

  u%cbar%m22%d(:)%ref = pc_%x_cbar22
  u%cbar%m12%d(:)%ref = pc_%x_cbar12
  u%cbar%m11%d(:)%ref = pc_%y_cbar11

  u%cbar%m11%d(:)%good_ref = pc_%ok_x .and. pc_%ok_y
  u%cbar%m12%d(:)%good_ref = pc_%ok_x .and. pc_%ok_y
  u%cbar%m22%d(:)%good_ref = pc_%ok_x .and. pc_%ok_y

  u%cmat_a%m22%d(:)%good_ref = pc_%ok_x .and. pc_%ok_y
  u%cmat_a%m12%d(:)%good_ref = pc_%ok_x .and. pc_%ok_y
  u%cmat_b%m11%d(:)%good_ref = pc_%ok_x .and. pc_%ok_y
  u%cmat_b%m12%d(:)%good_ref = pc_%ok_x .and. pc_%ok_y

  do i = 0, n_det_maxx

    ie = u%orbit%x%d(i)%ix_ele
    r_beta = sqrt(u%ring%ele(ie)%b%beta / u%ring%ele(ie)%a%beta)
    u%cmat_a%m22%d(i)%ref = pc_(i)%x_cbar22 * r_beta
    u%cmat_a%m12%d(i)%ref = pc_(i)%x_cbar12 * r_beta
    u%cmat_b%m11%d(i)%ref = pc_(i)%y_cbar11 / r_beta
    u%cmat_b%m12%d(i)%ref = pc_(i)%y_cbar12 / r_beta

    if (all(orbit_(i)%amp(:) > 100)) then
      !! u%orbit%x%d(i)%good_ref = .true.
      !! u%orbit%y%d(i)%good_ref = .true.
    endif
  enddo

  u%phase%x%d(5)%good_ref = .false.
  u%phase%y%d(5)%good_ref = .false.
  u%cbar%m11%d(5)%good_ref = .false.
  u%cbar%m12%d(5)%good_ref = .false.
  u%cbar%m22%d(5)%good_ref = .false.

  u%beta%x%d(:)%ref = beta_meas%beta_a
  u%beta%y%d(:)%ref = beta_meas%beta_b
  u%beta%x%d(:)%good_ref = beta_meas%ok_beta_a
  u%beta%y%d(:)%good_ref = beta_meas%ok_beta_b

  print *, 'Reference Phase/Coupling/Orbit Data Read In'

  if (logic%ref_particle == electron$) then  ! reverse phases for electrons
    u%phase%x%d(:)%ref = u%tune%x%d(1)%meas - u%phase%x%d(:)%ref
    u%phase%y%d(:)%ref = u%tune%y%d(1)%meas - u%phase%y%d(:)%ref
  endif

  u%phase%ref_file_name = file_name
  u%cbar%ref_file_name = file_name
  u%cmat_a%ref_file_name = file_name
  u%cmat_b%ref_file_name = file_name
  u%orbit%ref_file_name = file_name

  u%phase%ref_measured = .true.
  u%cbar%ref_measured = .true.
  u%beta%ref_measured = .true.
  u%orbit%ref_measured = .true.

  call baseline_set (plot_ref$, add$, graph%bottom1, graph%top1)
  u%beta%p2%base = plot_ref$

end select

!---------------------------------------------------------------------
! misc bookkeeping

if (logic%opt_vars /= opt_custom$ .and. logic%opt_vars /= opt_sex$) &
                                    call opt_vars_set(opt_quad$, .false.)

err_flag = .false.

end subroutine