subroutine read_tbt (data_or_ref, num_in, u, graph, err_flag)

use cesrv_struct
use cesrv_interface
use tbt_mod
use file_number_mod

implicit none

type (universe_struct), target :: u
type (graph_struct) graph
type (tbt_all_data_struct), pointer :: tbt_all
type (tbt_bunch_struct), pointer :: tbt_bunch
type (sex_res_all_data_struct), target :: res_all
type (sex_res_bunch_struct), pointer :: res_bunch
type (sex_res_bpm_struct), pointer :: res_bpm

character(140) file0, file1, path, basename

real(rp) sqrt_ba, sqrt_bb, sum_data, sum_model, b2_ratio, a2_x, a2_y

integer  data_or_ref, num_in, num, idum, i, ix, ie

logical err_flag, err

! Read

if (data_or_ref == data_file$) then
  tbt_all => u%tbt_meas
else
  tbt_all => u%tbt_ref
endif

call number_to_true_number (num_in, 'TBT', num, err_flag) 
if (err_flag) return
call tbt_read_data_file(num-1, gain_and_pedestal_correct$, .true., tbt_all, err_flag)
if (err_flag) return
call tbt_read_data_file(num, gain_and_pedestal_correct$, .false., tbt_all, err_flag)
if (err_flag) return

call tbt_resonance_analysis (tbt_all, res_all, logic%filter_beta_res)

res_bunch => res_all%bunch(logic%ix_bunch_fft)
tbt_bunch => tbt_all%bunch(logic%ix_bunch_fft)

call number_to_file_name (num-1, 'tbt', file0, idum, err)
call number_to_file_name (num,   'tbt', file1, idum, err)

tbt_all%data_file = file1
tbt_all%file_number = num

ix = splitfilename (file1, path, basename)

call read_save_set_cu ('CONDX', tbt_all%condx_number, data_or_ref, u, err_flag) 

!---------------------------------------------------------------------
! Data

if (data_or_ref == data_file$) then

  logic%opt_vars = opt_sex$

  print *, ' TBT DATA READ IN: ' 
  print *, '    ', trim(file0)
  print *, '    ', trim(file1)

  u%beta%file_name        = basename
  u%x_fft%file_name       = basename
  u%y_fft%file_name       = basename
  u%q2x%file_name         = basename
  u%q2y%file_name         = basename
  u%qx_plus_qy%file_name  = basename
  u%qx_minus_qy%file_name = basename

  u%beta%measured        = .true.
  u%x_fft%measured       = .true.
  u%y_fft%measured       = .true.
  u%q2x%measured         = .true.
  u%q2y%measured         = .true.
  u%qx_plus_qy%measured  = .true.
  u%qx_minus_qy%measured = .true.

  u%beta%ix_meas        = num
  u%x_fft%ix_meas       = num
  u%y_fft%ix_meas       = num
  u%q2x%ix_meas         = num
  u%q2y%ix_meas         = num
  u%qx_plus_qy%ix_meas  = num
  u%qx_minus_qy%ix_meas = num

  do i = 0, 120
    ie = u%q2x%a_in%d(i)%ix_ele

    A2_x = res_bunch%bpm(i)%qa_x%A_in**2 + res_bunch%bpm(i)%qa_x%A_out**2
    A2_y = res_bunch%bpm(i)%qb_y%A_in**2 + res_bunch%bpm(i)%qb_y%A_out**2

    if (A2_x == 0 .or. A2_y == 0) cycle

    sqrt_ba = sqrt(u%ring%ele(ie)%a%beta / A2_x)
    sqrt_bb = sqrt(u%ring%ele(ie)%b%beta / A2_y)

    u%beta%x%d(i)%meas          = res_bunch%bpm(i)%qa_x%A_in**2 + res_bunch%bpm(i)%qa_x%A_out**2
    u%beta%y%d(i)%meas          = res_bunch%bpm(i)%qb_y%A_in**2 + res_bunch%bpm(i)%qb_y%A_out**2

    u%q2x%a_in%d(i)%meas         = res_bunch%bpm(i)%q2a%A_in         * sqrt_ba * sqrt_ba
    u%q2y%a_in%d(i)%meas         = res_bunch%bpm(i)%q2b%A_in         * sqrt_bb * sqrt_bb
    u%qx_plus_qy%a_in%d(i)%meas  = res_bunch%bpm(i)%qa_plus_qb%A_in  * sqrt_ba * sqrt_bb
    u%qx_minus_qy%a_in%d(i)%meas = res_bunch%bpm(i)%qa_minus_qb%A_in * sqrt_ba * sqrt_bb

    u%q2x%a_out%d(i)%meas         = res_bunch%bpm(i)%q2a%A_out         * sqrt_ba * sqrt_ba
    u%q2y%a_out%d(i)%meas         = res_bunch%bpm(i)%q2b%A_out         * sqrt_bb * sqrt_bb
    u%qx_plus_qy%a_out%d(i)%meas  = res_bunch%bpm(i)%qa_plus_qb%A_out  * sqrt_ba * sqrt_bb
    u%qx_minus_qy%a_out%d(i)%meas = res_bunch%bpm(i)%qa_minus_qb%A_out * sqrt_ba * sqrt_bb
  enddo

  u%beta%x%d%good_dat           = tbt_bunch%bpm%data_ok
  u%beta%y%d%good_dat           = tbt_bunch%bpm%data_ok

  u%q2x%a_in%d%good_dat         = tbt_bunch%bpm%data_ok
  u%q2y%a_in%d%good_dat         = tbt_bunch%bpm%data_ok
  u%qx_plus_qy%a_in%d%good_dat  = tbt_bunch%bpm%data_ok
  u%qx_minus_qy%a_in%d%good_dat = tbt_bunch%bpm%data_ok

  u%q2x%a_out%d%good_dat         = tbt_bunch%bpm%data_ok
  u%q2y%a_out%d%good_dat         = tbt_bunch%bpm%data_ok
  u%qx_plus_qy%a_out%d%good_dat  = tbt_bunch%bpm%data_ok
  u%qx_minus_qy%a_out%d%good_dat = tbt_bunch%bpm%data_ok

  u%q2x%data_params%x_amp         = sum(sqrt(res_bunch%bpm%q2a%A_in**2 + res_bunch%bpm%q2a%A_out**2), &
                                                  mask = tbt_bunch%bpm%data_ok) / count(tbt_bunch%bpm%data_ok)
  u%q2y%data_params%x_amp         = sum(sqrt(res_bunch%bpm%q2b%A_in**2 + res_bunch%bpm%q2b%A_out**2), &
                                                  mask = tbt_bunch%bpm%data_ok) / count(tbt_bunch%bpm%data_ok)
  u%qx_plus_qy%data_params%x_amp  = sum(sqrt(res_bunch%bpm%qa_plus_qb%A_in**2 + res_bunch%bpm%qa_plus_qb%A_out**2), &
                                                  mask = tbt_bunch%bpm%data_ok) / count(tbt_bunch%bpm%data_ok)
  u%qx_minus_qy%data_params%x_amp = sum(sqrt(res_bunch%bpm%qa_minus_qb%A_in**2 + res_bunch%bpm%qa_minus_qb%A_out**2), &
                                                  mask = tbt_bunch%bpm%data_ok) / count(tbt_bunch%bpm%data_ok)

  ! Renormalize beta

  if (count(u%beta%x%d%good_user) == 0) then
    u%beta%x%d%good_temp = u%beta%x%d%good_dat
  else
    u%beta%x%d%good_temp = u%beta%x%d%good_dat .and. u%beta%x%d%good_user 
  endif

  sum_data  = sum(sqrt(u%beta%x%d%meas), mask = u%beta%x%d%good_temp)
  sum_model = sum(sqrt(u%beta%x%d%model), mask = u%beta%x%d%good_temp)
  u%beta%x%d%meas = u%beta%x%d%meas * (sum_model / sum_data)**2

  sum_data  = sum(sqrt(u%beta%y%d%meas), mask = u%beta%x%d%good_temp)
  sum_model = sum(sqrt(u%beta%y%d%model), mask = u%beta%x%d%good_temp)
  u%beta%y%d%meas = u%beta%y%d%meas * (sum_model / sum_data)**2

  !

  u%tune%x%d(1)%meas = res_bunch%tune_a 
  u%tune%x%d(1)%meas = u%tune%x%d(1)%meas + int(u%tune%x%d(1)%model/twopi) * twopi
  u%tune%x%d(1)%good_dat = .true.
  u%tune%data_params%x_amp = sum(sqrt(res_bunch%bpm%qa_x%A_in**2 + res_bunch%bpm%qa_x%A_out**2), &
                                                  mask = tbt_bunch%bpm%data_ok) / count(tbt_bunch%bpm%data_ok)


  u%tune%y%d(1)%meas = res_bunch%tune_b 
  u%tune%y%d(1)%meas = u%tune%y%d(1)%meas + int(u%tune%y%d(1)%model/twopi) * twopi
  u%tune%y%d(1)%good_dat = .true.
  u%tune%data_params%y_amp = sum(sqrt(res_bunch%bpm%qb_y%A_in**2 + res_bunch%bpm%qb_y%A_out**2), &
                                                  mask = tbt_bunch%bpm%data_ok) / count(tbt_bunch%bpm%data_ok)

  !

  if (logic%read_twiss_with_tbt) then

    u%phase%file_name = basename
    u%phase%measured  = .true.
    u%phase%ix_meas   = num
    u%phase%x%d%good_dat = tbt_bunch%bpm%data_ok
    u%phase%y%d%good_dat = tbt_bunch%bpm%data_ok

    u%cbar%file_name = basename
    u%cbar%measured  = .true.
    u%cbar%ix_meas   = num
    u%cbar%m11%d%good_dat = tbt_bunch%bpm%data_ok
    u%cbar%m12%d%good_dat = tbt_bunch%bpm%data_ok
    u%cbar%m22%d%good_dat = tbt_bunch%bpm%data_ok

    u%orbit%file_name = basename
    u%orbit%measured  = .true.
    u%orbit%ix_meas   = num
    u%orbit%x%d%good_dat = tbt_bunch%bpm%data_ok
    u%orbit%y%d%good_dat = tbt_bunch%bpm%data_ok

    do i = 0, 120

      if (.not. u%orbit%x%d(i)%exists) cycle
      if (.not. tbt_bunch%bpm(i)%data_ok) cycle
      res_bpm => res_bunch%bpm(i)

      u%phase%x%d(i)%meas = atan2(res_bpm%qa_x%A_out, res_bpm%qa_x%A_in)
      call this_phase_shift (u%phase%x%d%meas, u%phase%x%d%model, u%phase%x%d%good_dat)
 
      u%phase%y%d(i)%meas = atan2(res_bpm%qb_y%A_out, res_bpm%qb_y%A_in)
      call this_phase_shift (u%phase%y%d%meas, u%phase%y%d%model, u%phase%y%d%good_dat)

      b2_ratio = sqrt(u%beta%y%d(i)%model / u%beta%x%d(i)%model)

      u%cbar%m11%d(i)%meas =  (res_bpm%qb_x%A_in * res_bpm%qb_y%A_in + res_bpm%qb_x%A_out * res_bpm%qb_y%A_out) * b2_ratio / &
                                             (res_bpm%qb_y%A_in**2 + res_bpm%qb_y%A_out**2) 
      u%cbar%m12%d(i)%meas =  (res_bpm%qa_x%A_in * res_bpm%qa_y%A_out - res_bpm%qa_x%A_out * res_bpm%qa_y%A_in) / b2_ratio / &
                                             (res_bpm%qa_x%A_in**2 + res_bpm%qa_x%A_out**2)      
      u%cbar%m22%d(i)%meas = -(res_bpm%qa_x%A_in * res_bpm%qa_y%A_in + res_bpm%qa_x%A_out * res_bpm%qa_y%A_out) / b2_ratio / &
                                             (res_bpm%qa_x%A_in**2 + res_bpm%qa_x%A_out**2) 
      
      u%orbit%x%d(i)%meas = res_bpm%x_ave
      u%orbit%y%d(i)%meas = res_bpm%y_ave

    enddo

  endif

  ! Set plots
  ! set_plot will override %plot_data values so have to do things in the correct order.

  if (.not. this_graph_type (graph%top1%d2%type, .false.)) then
    call set_plot (graph%top1, u%q2x)
    call plot_data_set (graph%top1, plot_meas$)
  endif

  if (.not. this_graph_type (graph%bottom1%d2%type, .false.)) then
    call set_plot (graph%bottom1, u%q2y)
    call plot_data_set (graph%bottom1, plot_meas$)
  endif

  if (logic%wide_plot_window) then
    if (.not. this_graph_type (graph%top2%d2%type, .false.)) then
      call set_plot (graph%top2, u%qx_plus_qy)
      call plot_data_set (graph%top2, plot_meas$)
    endif
    if (.not. this_graph_type (graph%bottom2%d2%type, .false.)) then
      call set_plot (graph%bottom2, u%qx_minus_qy)
      call plot_data_set (graph%bottom2, plot_meas$)
    endif
  endif

  u%beta%p2%plot_data = plot_meas$
  if (logic%read_twiss_with_tbt) then
    u%phase%p2%plot_data = plot_meas$
    u%cbar%p2%plot_data = plot_meas$
    u%orbit%p2%plot_data = plot_meas$
  endif

  print '(a, f10.6)', 'Qx Amplitude:    ', u%tune%data_params%x_amp
  print '(a, f10.6)', 'Qy Amplitude:    ', u%tune%data_params%y_amp
  print '(a, f10.6)', '2Qx Amplitude:   ', u%q2x%data_params%x_amp
  print '(a, f10.6)', '2Qy Amplitude:   ', u%q2y%data_params%x_amp
  print '(a, f10.6)', 'Qx+Qy Amplitude: ', u%qx_plus_qy%data_params%x_amp
  print '(a, f10.6)', 'Qx-Qy Amplitude: ', u%qx_minus_qy%data_params%x_amp

!---------------------------------------------------------------------
! Ref

else

  print *, ' TBT REFERENCE READ IN:' 
  print *, '    ', trim(file0)
  print *, '    ', trim(file1)

  if (this_graph_type (graph%top1%d2%type, .true.))    call baseline_set (plot_ref$, add$, graph%top1)
  if (this_graph_type (graph%bottom1%d2%type, .true.)) call baseline_set (plot_ref$, add$, graph%bottom1)
  if (this_graph_type (graph%top2%d2%type, .true.))    call baseline_set (plot_ref$, add$, graph%top2)
  if (this_graph_type (graph%bottom2%d2%type, .true.)) call baseline_set (plot_ref$, add$, graph%bottom2)

  u%beta%ref_file_name        = basename
  u%x_fft%ref_file_name       = basename
  u%y_fft%ref_file_name       = basename
  u%q2x%ref_file_name         = basename
  u%q2y%ref_file_name         = basename
  u%qx_plus_qy%ref_file_name  = basename
  u%qx_minus_qy%ref_file_name = basename

  u%beta%ref_measured        = .true.
  u%x_fft%ref_measured       = .true.
  u%y_fft%ref_measured       = .true.
  u%q2x%ref_measured         = .true.
  u%q2y%ref_measured         = .true.
  u%qx_plus_qy%ref_measured  = .true.
  u%qx_minus_qy%ref_measured = .true.

  u%beta%ix_ref        = num
  u%x_fft%ix_ref       = num
  u%y_fft%ix_ref       = num
  u%q2x%ix_ref         = num
  u%q2y%ix_ref         = num
  u%qx_plus_qy%ix_ref  = num
  u%qx_minus_qy%ix_ref = num

  do i = 0, 120
    ie = u%q2x%a_in%d(i)%ix_ele

    A2_x = res_bunch%bpm(i)%qa_x%A_in**2 + res_bunch%bpm(i)%qa_x%A_out**2
    A2_y = res_bunch%bpm(i)%qb_y%A_in**2 + res_bunch%bpm(i)%qb_y%A_out**2

    if (A2_x == 0 .or. A2_y == 0) cycle

    sqrt_ba = sqrt(u%ring%ele(ie)%a%beta / A2_x)
    sqrt_bb = sqrt(u%ring%ele(ie)%b%beta / A2_y)

    u%beta%x%d(i)%ref          = res_bunch%bpm(i)%qa_x%A_in**2 + res_bunch%bpm(i)%qa_x%A_out**2
    u%beta%y%d(i)%ref          = res_bunch%bpm(i)%qb_y%A_in**2 + res_bunch%bpm(i)%qb_y%A_out**2

    u%q2x%a_in%d(i)%ref         = res_bunch%bpm(i)%q2a%A_in         * sqrt_ba * sqrt_ba
    u%q2y%a_in%d(i)%ref         = res_bunch%bpm(i)%q2b%A_in         * sqrt_bb * sqrt_bb
    u%qx_plus_qy%a_in%d(i)%ref  = res_bunch%bpm(i)%qa_plus_qb%A_in  * sqrt_ba * sqrt_bb
    u%qx_minus_qy%a_in%d(i)%ref = res_bunch%bpm(i)%qa_minus_qb%A_in * sqrt_ba * sqrt_bb

    u%q2x%a_out%d(i)%ref         = res_bunch%bpm(i)%q2a%A_out         * sqrt_ba * sqrt_ba
    u%q2y%a_out%d(i)%ref         = res_bunch%bpm(i)%q2b%A_out         * sqrt_bb * sqrt_bb
    u%qx_plus_qy%a_out%d(i)%ref  = res_bunch%bpm(i)%qa_plus_qb%A_out  * sqrt_ba * sqrt_bb
    u%qx_minus_qy%a_out%d(i)%ref = res_bunch%bpm(i)%qa_minus_qb%A_out * sqrt_ba * sqrt_bb

  enddo

  u%beta%x%d%good_ref           = tbt_bunch%bpm%data_ok
  u%beta%y%d%good_ref           = tbt_bunch%bpm%data_ok

  u%q2x%a_in%d%good_ref         = tbt_bunch%bpm%data_ok
  u%q2y%a_in%d%good_ref         = tbt_bunch%bpm%data_ok
  u%qx_plus_qy%a_in%d%good_ref  = tbt_bunch%bpm%data_ok
  u%qx_minus_qy%a_in%d%good_ref = tbt_bunch%bpm%data_ok

  u%q2x%a_out%d%good_ref         = tbt_bunch%bpm%data_ok
  u%q2y%a_out%d%good_ref         = tbt_bunch%bpm%data_ok
  u%qx_plus_qy%a_out%d%good_ref  = tbt_bunch%bpm%data_ok
  u%qx_minus_qy%a_out%d%good_ref = tbt_bunch%bpm%data_ok

  u%q2x%ref_params%x_amp         = sum(sqrt(res_bunch%bpm%q2a%A_in**2 + res_bunch%bpm%q2a%A_out**2), &
                                                  mask = tbt_bunch%bpm%data_ok) / count(tbt_bunch%bpm%data_ok)
  u%q2y%ref_params%x_amp         = sum(sqrt(res_bunch%bpm%q2b%A_in**2 + res_bunch%bpm%q2b%A_out**2), &
                                                  mask = tbt_bunch%bpm%data_ok) / count(tbt_bunch%bpm%data_ok)
  u%qx_plus_qy%ref_params%x_amp  = sum(sqrt(res_bunch%bpm%qa_plus_qb%A_in**2 + res_bunch%bpm%qa_plus_qb%A_out**2), &
                                                  mask = tbt_bunch%bpm%data_ok) / count(tbt_bunch%bpm%data_ok)
  u%qx_minus_qy%ref_params%x_amp = sum(sqrt(res_bunch%bpm%qa_minus_qb%A_in**2 + res_bunch%bpm%qa_minus_qb%A_out**2), &
                                                  mask = tbt_bunch%bpm%data_ok) / count(tbt_bunch%bpm%data_ok)

  ! Renormalize beta

  if (count(u%beta%x%d%good_user) == 0) then
    u%beta%x%d%good_temp = u%beta%x%d%good_dat
  else
    u%beta%x%d%good_temp = u%beta%x%d%good_dat .and. u%beta%x%d%good_user 
  endif

  sum_data  = sum(sqrt(u%beta%x%d%ref), mask = u%beta%x%d%good_temp)
  sum_model = sum(sqrt(u%beta%x%d%model), mask = u%beta%x%d%good_temp)
  u%beta%x%d%ref = u%beta%x%d%ref * (sum_model / sum_data)**2

  sum_data  = sum(sqrt(u%beta%y%d%ref), mask = u%beta%x%d%good_temp)
  sum_model = sum(sqrt(u%beta%y%d%model), mask = u%beta%x%d%good_temp)
  u%beta%y%d%ref = u%beta%y%d%ref * (sum_model / sum_data)**2

  !

  u%tune%x%d(1)%ref = res_bunch%tune_a 
  u%tune%x%d(1)%ref = u%tune%x%d(1)%ref + int(u%tune%x%d(1)%model/twopi) * twopi
  u%tune%x%d(1)%good_ref = .true.
  u%tune%ref_params%x_amp = sum(sqrt(res_bunch%bpm%qa_x%A_in**2 + res_bunch%bpm%qa_x%A_out**2), &
                                                  mask = tbt_bunch%bpm%data_ok) / count(tbt_bunch%bpm%data_ok)

  u%tune%y%d(1)%ref = res_bunch%tune_b 
  u%tune%y%d(1)%ref = u%tune%y%d(1)%ref + int(u%tune%y%d(1)%model/twopi) * twopi
  u%tune%y%d(1)%good_ref = .true.
  u%tune%ref_params%y_amp = sum(sqrt(res_bunch%bpm%qb_y%A_in**2 + res_bunch%bpm%qb_y%A_out**2), &
                                                  mask = tbt_bunch%bpm%data_ok) / count(tbt_bunch%bpm%data_ok)

  !

  if (logic%read_twiss_with_tbt) then

    u%phase%ref_file_name = basename
    u%phase%ref_measured  = .true.
    u%phase%ix_ref   = num
    u%phase%x%d%good_ref = tbt_bunch%bpm%data_ok
    u%phase%y%d%good_ref = tbt_bunch%bpm%data_ok

    u%cbar%ref_file_name = basename
    u%cbar%ref_measured  = .true.
    u%cbar%ix_ref   = num
    u%cbar%m11%d%good_ref = tbt_bunch%bpm%data_ok
    u%cbar%m12%d%good_ref = tbt_bunch%bpm%data_ok
    u%cbar%m22%d%good_ref = tbt_bunch%bpm%data_ok

    u%orbit%ref_file_name = basename
    u%orbit%ref_measured  = .true.
    u%orbit%ix_ref   = num
    u%orbit%x%d%good_ref = tbt_bunch%bpm%data_ok
    u%orbit%y%d%good_ref = tbt_bunch%bpm%data_ok

    do i = 0, 120

      if (.not. u%orbit%x%d(i)%exists) cycle
      if (.not. tbt_bunch%bpm(i)%data_ok) cycle
      res_bpm => res_bunch%bpm(i)

      u%phase%x%d(i)%ref = atan2(res_bpm%qa_x%A_out, res_bpm%qa_x%A_in)
      call this_phase_shift (u%phase%x%d%ref, u%phase%x%d%model, u%phase%x%d%good_ref)
 
      u%phase%y%d(i)%ref = atan2(res_bpm%qb_y%A_out, res_bpm%qb_y%A_in)
      call this_phase_shift (u%phase%y%d%ref, u%phase%y%d%model, u%phase%y%d%good_ref)

      b2_ratio = sqrt(u%beta%y%d(i)%model / u%beta%x%d(i)%model)

      u%cbar%m11%d(i)%ref = (res_bpm%qb_x%A_in * res_bpm%qb_y%A_in + res_bpm%qb_x%A_out * res_bpm%qb_y%A_out) / &
                                             (res_bpm%qb_y%A_in**2 + res_bpm%qb_y%A_out**2) / b2_ratio
      u%cbar%m12%d(i)%ref = (res_bpm%qa_x%A_in * res_bpm%qa_y%A_out - res_bpm%qa_x%A_out * res_bpm%qa_y%A_in) * b2_ratio / &
                                             (res_bpm%qa_x%A_in**2 + res_bpm%qa_x%A_out**2)      
      u%cbar%m22%d(i)%ref = (res_bpm%qa_x%A_in * res_bpm%qa_y%A_in + res_bpm%qa_x%A_out * res_bpm%qa_y%A_out) * b2_ratio / &
                                             (res_bpm%qa_x%A_in**2 + res_bpm%qa_x%A_out**2) 
      
      u%orbit%x%d(i)%ref = res_bpm%x_ave
      u%orbit%y%d(i)%ref = res_bpm%y_ave

    enddo

  endif

  print '(a, f10.6)', 'Qx Amplitude:    ', u%tune%ref_params%x_amp
  print '(a, f10.6)', 'Qy Amplitude:    ', u%tune%ref_params%y_amp
  print '(a, f10.6)', '2Qx Amplitude:   ', u%q2x%ref_params%x_amp
  print '(a, f10.6)', '2Qy Amplitude:   ', u%q2y%ref_params%x_amp
  print '(a, f10.6)', 'Qx+Qy Amplitude: ', u%qx_plus_qy%ref_params%x_amp
  print '(a, f10.6)', 'Qx-Qy Amplitude: ', u%qx_minus_qy%ref_params%x_amp

endif

!-----------------------------------------------------------------------------
contains
subroutine this_phase_shift (meas, design, good_dat)

implicit none

real(rp) meas(:), design(:)
real(rp) offset, rms, rms_min

integer n, i, i_min

logical good_dat(:)

!

rms_min = 1e10
n = count(good_dat)
do i = 0, 7
  offset = i * twopi / 8
  meas = (design + offset) + modulo2(meas - (design + offset), pi)
  rms = sum(meas - (design + offset))**2
  if (rms < rms_min) then
    rms_min = rms
    i_min = i
  endif
enddo

offset = i_min * twopi / 8
meas = (design + offset) + modulo2(meas - (design + offset), pi)

end subroutine

!-----------------------------------------------------------------------------
! contains

function this_graph_type (graph_type, include_twiss) result (is_res_type)

implicit none

integer graph_type
logical include_twiss, is_res_type

!

select case (graph_type)
case (q2x_data$, q2y_data$, qx_minus_qy_data$, qx_plus_qy_data$)
  is_res_type = .true.
case default
  is_res_type = .false.
end select

if (include_twiss) is_res_type = .true.

end function

end subroutine