module phase_mod

use cesr_basic_mod
use sim_utils
use cesrv_struct
use cesrv_interface
use cesr_tune_tracker_io_mod
use cesr_read_data_mod
use nr

include 'cbpm_cesr_interface.inc'

real(rp), save :: amp_min_default = 1000

contains

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Input:
!   who        -- Character(*): 'AC_ETA' or 'PHASE'.
!   auto_meas  -- Logical: If true then don't ask any questions.


subroutine take_phase_measurement (p, skew_quad_on, u, who, &
                                     measurement_taken, auto_meas, comment)

implicit none

type (cesr_freq_struct),save :: freq
type (cesr_phase_params_struct) :: p
type (cesr_det_plane_struct) horz(0:n_det_maxx)
type (cesr_det_plane_struct) vert(0:n_det_maxx)
type (cesr_det_dc_position_struct) dc(0:n_det_maxx)

type (universe_struct) u
type (cesr_all_data_struct) ac_eta

character(*) who
character(80) line, prompt
character(9) cur_name(-1:1) / 'Electrons', ' ', 'Positrons' /
character(5) string
character(140) data_file_number, number_file

real(rp) a_frac_tune, b_frac_tune
real(rp) tune_x, tune_y, cur

integer i, csrbpm_comd(10), atten, stat, detec, bpm_type
integer det, ix, ios, id, it, it1, it2
integer species_list(CBPM_MX_BUNCH_LIST)
integer train_list(CBPM_MX_BUNCH_LIST)
integer bunch_list(CBPM_MX_BUNCH_LIST)
integer take_cbpm_phase, take_cbpm_orbit

integer raw(4, 120), current(2)
integer :: comd(16) = (/ (0, i = 1, 16) /)             ! command array for BEARUN
integer :: termout = -1 

logical :: skew_quad_on(:), auto_meas, abort_wanted
logical not_standard
logical measurement_taken, no_current, err_flag, ok

character(*), optional :: comment
character(1) :: bell = char(7)

! init
                               
freq%rev = 390.14
measurement_taken = .false.

! electrons or positrons?

call set_species (p%species, err_flag)
if (err_flag) return

call check_current(ok, cur)
if (.not. ok) return

no_current = .false.
if (cur < 20e-6) then
  print *, 'WARNING: YOU DON''T SEEM TO HAVE ANY CURRENT!'
  no_current = .true.
endif

print *, 'Beam current setting is: ', cur * 1e3

! If the analog tune tracker is being used then assume we 
! are working at the reflected frequencies.

freq%x%reflection = .false.
freq%y%reflection = .false.

if (who == 'PHASE' .and. .not. logic%digital_shaker_on) then
  freq%x%reflection = .true.
  freq%y%reflection = .true.
endif

call getlat (p%lattice)

if (who == 'AC_ETA') then
  call set_phase_meas_freq_range (logic%ac_eta_freq_range) 
  number_file = '$CESR_ONLINE/machine_data/mach_meas/ac_eta/ac_eta.number'
elseif (who == 'PHASE') then
  call set_phase_meas_freq_range (0) 
  number_file = '$CESR_ONLINE/machine_data/mach_meas/phase/phase.number'
else
  print *, 'INTERNAL ERROR!'
  call err_exit
endif

! Standard settings?

10 continue

if (.not. auto_meas) then
  print *, 'Remember!'
  print *, '  * Take a CSR save set? (save set number is recorded in data file)'
  print *, '  * You need to use Train 1, Car 1.'
  print *, '  * Set EXTEC to look at Train 1, Car 1.'
  !! print *, '  * Set EXTEC to the correct bunch current.'
  print *, '  * Phase measurement: Nickolet mode 8, AC_Eta measurement: mode 3.'
  print *, '  * AC_Eta measurement: Lock on Synchrotron tune.'

  call get_input_string (' Use Standard Settings? <CR = Yes>: ', line)
  call string_trim (line, line, ix)
  call str_upcase (line, line)
  if (ix == 0 .or. line(1:1) == 'Y') then
    not_standard = .false.
  elseif (line(1:1) == 'N') then  
    not_standard = .true.
  else
    print *, 'I Don''t understand this.'
    goto 10
  endif
else
  not_standard = .false.
endif

! check the setup

call check_setup (skew_quad_on, p%species, who, abort_wanted)
if (abort_wanted) return

! In the case where there was current before: Check current again just in case.

if (.not. no_current) then
  call check_current(ok)
  if (.not. ok) return
endif


! Find if reflection shake

if (not_standard) then
  print *
  print *, 'Are the tune trackers locked to the Tune or the Reflection?'

  do
    if (.not. freq%x%reflection) then
       prompt = ' HORIZONTAL tracker: Tune or Reflection? <CR = T>: '
    else
       prompt = ' HORIZONTAL tracker: Tune or Reflection? <CR = R>: '
    endif
    call get_input_string (prompt, line)
    call str_upcase(line, line)
    call string_trim (line, line, ix)

    if (line(1:1) == 'R') then  
      freq%x%reflection = .true.
      exit
    elseif (line(1:1) == 'T') then
      freq%x%reflection = .false.
      exit
    elseif (ix == 0) then ! ix = 0 ==> nothing typed, use default
      exit
    else
      print *
      print *, 'Please enter "T" or "R"  (or nothing), Try again...'
    endif
  enddo

  do
    if (.not. freq%y%reflection) then
       prompt = ' VERTICAL tracker: Tune or Reflection? <CR = T>: '
    else
       prompt = ' VERTICAL tracker: Tune or Reflection? <CR = R>: '
    endif
    call get_input_string (prompt, line)
    call str_upcase(line, line)
    call string_trim (line, line, ix)

    if (line(1:1) == 'R') then  
      freq%y%reflection = .true.
      exit
    elseif (line(1:1) == 'T') then
      freq%y%reflection = .false.
      exit
    elseif (ix == 0) then ! ix = 0 ==> nothing typed, use default
      exit
    else
      print *
      print *, 'Please enter "T" or "R"  (or nothing), Try again...'
    endif
  enddo

endif

call set_digital_tracker_tune_status (freq%x%reflection, freq%y%reflection)

! Get tunes and tracking

call get_freq (freq, 5)
print *
print '(a, f8.3)', ' freq horz:', freq%x%tune
if (who == 'PHASE') print '(a, f8.3)', ' freq vert:', freq%y%tune

a_frac_tune = u%design%ele(u%design%n_ele_track)%a%phi / twopi
a_frac_tune = a_frac_tune - int(a_frac_tune)
b_frac_tune = u%design%ele(u%design%n_ele_track)%b%phi / twopi
b_frac_tune = b_frac_tune - int(b_frac_tune)
tune_x = freq%x%tune / 390.
tune_y = freq%y%tune / 390.


if (who == 'PHASE') then
  if (abs(tune_x - a_frac_tune) + abs(tune_y - b_frac_tune) .gt. &
          abs(tune_x - b_frac_tune) + abs(tune_y - a_frac_tune)) then
    print *
    print *, '!-------------------------------------------------------------'
    print *, '! Warning: Based on the tracker frequencies you seem to have -'
    print *, '!          Horizontal and Vertical modes switched!           -'
    print *, '!-------------------------------------------------------------'
    call get_input_string ('Continue? (Y/N)', line)
    call string_trim(line, line, ix)
    if (line(1:1) /= 'Y' .and. line(1:1) /= 'y') goto 10
  endif
endif

call vxgetn ('CSR SAVRECRD', 1, 1, p%save_set)  ! get save set number

!

1000 continue

species_list(1) = p%species
train_list(1) = 1
bunch_list(1) = 1

! complete measurement

print *
if (present(comment)) then
  print *, 'Comment is: ', trim(comment)
  p%comment = comment
else
  print *, 'Type "QUIT" to NOT take a measurement.'
  call get_input_string (' Optional comment:', p%comment) 
endif

if (p%comment == 'quit' .or. p%comment == 'QUIT') goto 8000  ! And end

call vxgetn ('CSR SAVRECRD', 1, 1, p%save_set)  ! get save set number

call vxputn ('CSRBPM PSTAT', 1, 120, CBPM_NO_DATA) ! Clear status
stat = take_cbpm_phase (species_list, train_list, bunch_list, CBPM_PROCESS_WAIT)

call all_phase_measurement (p, horz, vert, dc)

call set_pstat_stale()

measurement_taken = .true.
p%unit = lunget()
call fullfilename(number_file, data_file_number)
call increment_file_number (data_file_number, 5, ix, string)
call form_file_name_with_number (who, ix, p%file_name, err_flag)

if (who == 'AC_ETA') then
  ! Analyze without Wolski normal mode calc so that the data file is not affected.
  ! The Wolski calc is done when the ac_eta data is read in from the file.
  call analyze_ac_eta_data (horz, freq, dc, p, ac_eta, u, .false.)
else
  call analyze_phase_data (horz, vert, freq, dc, p, u)
endif

call write_shaking_data (horz, vert, freq, dc, p, ac_eta, who)

! type out info on non OK detectors, etc to the terminal
! and signal end of measurement

print *
print *, '      |   OK  | N_good | '
print *, '  det |  x  y |  x  y  | int(x_pos/5mm)'

do det = 0, 120
  if (det > 99 .and. all(horz(det)%but(:)%amp == 0) .and. &
                              all(vert(det)%but(:)%amp == 0)) cycle
  if (horz(det)%n_buts /= 4 .or. .not. horz(det)%ok .or.  &
      vert(det)%n_buts /= 4 .or. .not. vert(det)%ok) then
    print '(i5, 2x, 2l3, 2x, 2i3, i8)', det, horz(det)%ok, vert(det)%ok, &
             horz(det)%n_buts, vert(det)%n_buts, int(dc(det)%x/5)
  endif
enddo

print *
print *, 'Written: ', trim(p%file_name)
print *, bell, bell, bell, bell
print *
print *, '************* Measurement Complete ****************'
print *

! End

8000  continue
if (who == 'AC_ETA') then
  call set_phase_meas_freq_range (0) 
endif

end subroutine take_phase_measurement

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

subroutine all_phase_measurement (p, horz, vert, dc)

implicit none

type (cesr_phase_params_struct) p
type (cesr_det_plane_struct), target :: horz(0:n_det_maxx), vert(0:n_det_maxx)
type (cesr_det_dc_position_struct) dc(0:n_det_maxx)
type (butns_struct) butns
                              
integer i, bpm_type(120), bpm_flag(120), data(1920), stat(120), raw(4,120)
integer det, id, ix, but

real(rp) ch1, cv1, ch2, cv2, dummy

logical err_flag

! init 

call vxgetn ('CSRBPM TYPE ', 1, 120, bpm_type)
call vxgetn ('CSRBPM FLAGS', 1, 120, bpm_flag)

horz(1:120)%system_id = bpm_type; horz(0)%system_id = bpm_type(100)
vert%system_id = horz%system_id

if (.not. logic%auto_measurement) then
  print *, &
       '         |   Ave Amplitude  | Horz Check  | Vert Check', &
       '    Det  |   Horz      Vert | 1/4    2/3  | 1/4    2/3'
endif

! Phase data.

call vxgetn ('CSRBPM PDATA', 1, 1920, data)
call vxgetn ('CSRBPM PSTAT', 1, 120, stat)

do det = 0, 120

  id = det
  if (det == 0) id = 100

  do but = 1, 4
    ix = 16*(id - 1) + 4*(but - 1) 
    horz(det)%but(but)%amp   = data(ix+1)
    horz(det)%but(but)%phase = data(ix+2) / 10000.0
    vert(det)%but(but)%amp   = data(ix+3)
    vert(det)%but(but)%phase = data(ix+4) / 10000.0
  enddo

  horz(det)%ok = (stat(id) == CBPM_GOOD_DATA) .and. (bpm_flag(id) == 1)
  vert(det)%ok = (stat(id) == CBPM_GOOD_DATA) .and. (bpm_flag(id) == 1)

  if (bpm_type(id) /= 2 .and. bpm_type(id) /= -1) then
    print *, 'ERROR: BAD CSRBPM_TYPE VALUE:', bpm_type(id), det 
  endif

  if (bpm_flag(id) /= 0 .and. bpm_flag(id) /= 1) then
    print *, 'ERROR: BAD CSRBPM_FLAGS VALUE:', bpm_flag(id), det 
  endif

  if (bpm_flag(id) == 1 .and. bpm_type(id) == -1) then
    print *, 'ERROR: CSRBPM_FLAGS INDICATES EXISTANT BPM BUT BPM_TYPE GIVES VALUE OF -1:', det
  endif

  ! type results

  if (.not. logic%auto_measurement) then

    ch1 = 180 - abs(horz(det)%but(1)%phase - horz(det)%but(4)%phase)
    cv1 = 180 - abs(vert(det)%but(1)%phase - vert(det)%but(4)%phase)
    ch2 = 180 - abs(horz(det)%but(2)%phase - horz(det)%but(3)%phase)
    cv2 = 180 - abs(vert(det)%but(2)%phase - vert(det)%but(3)%phase)

    if (bpm_flag(id) == 0 .and. det > 99) cycle

    if (bpm_flag(id) == 0) then
      print '(i8, a)', det, '  ---------------- Does Not Exist --------------'
    else
      print '(a, i4, 2i10, 4f7.0)', ' H/V', det, nint(sum(horz(det)%but(:)%amp) / 4000), &
           nint(sum(vert(det)%but(:)%amp) / 4000), ch1, ch2, cv1, cv2
    endif

  endif
  
enddo

! DC orbit data

call vxgetn ('BPM ORB DATA', 1, 840, data)
call vxgetn ('CSRBPM DSTAT', 1, 120, stat)

call poscon2 (-4, 0.0_rp, 0.0_rp, 0.0_rp, 0.0_rp, dummy, dummy)
call poscon2 (-2, 0.0_rp, 0.0_rp, 0.0_rp, 0.0_rp, dummy, dummy)

do det = 0, 120

  id = det
  if (det == 0) id = 100

  ix = 7*(id - 1)

  raw(:,id) = [data(ix+2), data(ix+3), data(ix+1), data(ix+4)]

  dc(det)%signal(1) = data(ix+2)
  dc(det)%signal(2) = data(ix+3)
  dc(det)%signal(3) = data(ix+1)
  dc(det)%signal(4) = data(ix+4)

!!  call poscon2 (id, dc(det)%signal(1), dc(det)%signal(2),  &
!!                    dc(det)%signal(3), dc(det)%signal(4), dc(det)%y, dc(det)%x)
!!  dc(det)%x = 1000.0 * dc(det)%x
!!  dc(det)%y = 1000.0 * dc(det)%y

enddo

!

call raw_butns_to_position (raw, bpm_type, butns, err_flag, &
                logic%nonlinear_calc, logic%offset_correction, logic%gain_correction)

do det = lbound(butns%det, 1), ubound(butns%det, 1)
  dc(det)%x = butns%det(det)%x_orb
  dc(det)%y = butns%det(det)%y_orb
  dc(det)%ok = butns%det(det)%ok
enddo


end subroutine all_phase_measurement

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

subroutine to_amp_phase (sine_in, cosine_in, amp, phase)

implicit none

integer sine_in, cosine_in
real(rp) amp, phase
real(rp) sine, cosine

!

sine = sine_in
cosine = cosine_in

amp = sqrt(sine**2 + cosine**2)

if (amp == 0) then
  phase = 0.0
else
  phase = atan2d (sine, cosine)
endif

end subroutine to_amp_phase

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

subroutine write_file_header (what, freq, p, ac_eta, who)

implicit none

type (cesr_freq_struct) freq
type (cesr_phase_params_struct) p
type (cesr_all_data_struct) ac_eta

integer j, ix, iu, vnumbr, ios
integer i, ival(100)

character(*) what, who
character(20) date_str
character(10) species_name
character(100) :: comment_out, route_ctl
character(40) line, route_name

integer, parameter :: db_num = 15
character(12) :: name, db_name(db_num) =  [ &
        'CSR_QUAD_CUR', 'CSR_QADD_CUR', 'CSR_SEXT_CUR', 'CSR_SQEWQUAD', &
        'CSR_SQEWSEXT', 'CSR_HORZ_CUR', 'CSR_HBND_CUR', 'CSR_VERT_CUR', &
        'CSR_HSP_VOLT', 'CSR_HSP_VVAL', 'CSR_OCTU_CUR', 'SCW_CUR_READ', &
        'UND_VERT_CUR', 'UND_CNTG_CUR', 'UND_CNTGTRIM']

! open file

open (unit = p%unit, file = p%file_name)

! get header info

if (p%species == -1) then
  species_name = 'ELECTRONS'
else
  species_name = 'POSITRONS'
endif

call date_and_time_stamp(date_str)
call make_legal_comment(p%comment, comment_out)

iu = lunget()
call fullfilename('$CESR_ONLINE/acc_control/program_info/route/route.ctl', &
     route_ctl)
open (iu, iostat = ios, file = route_ctl, status = 'old', action = 'read')
if (ios /= 0) then
  print *, 'ERROR: CANNOT OPEN: ', trim(route_ctl)
  route_name = ' '
else
  read (iu, '(a)') route_name
  read (iu, '(a)') line
  close (iu)
  route_name = trim(route_name) // '_' // line
  call string_trim (route_name, route_name, ix)
  call str_upcase (route_name, route_name)
endif

! write header info

iu = p%unit
write (iu, *) '&DATA_PARAMETERS'
write (iu, *) '  file_type  = ', "'", who, " DATA'"
write (iu, *) '  data_date  = ', "'", date_str, "'"                
write (iu, *) '  lattice    = ', "'", trim(p%lattice), "'"
write (iu, *) '  save_set   = ', p%save_set
write (iu, *) '  route_name = ', "'", trim(route_name), "'"
write (iu, *) '  comment = ', trim(comment_out)
write (iu, *) '/END'
write (iu, *) 
write (iu, *) '&PHASE_PARAMETERS'
write (iu, *) '  horiz_reflection_shake = ', freq%x%reflection
write (iu, *) '  vert_reflection_shake  = ', freq%y%reflection
write (iu, *) '  species          = ', p%species, '  ! ', species_name
write (iu, *) '  horiz_beta_freq  = ', freq%x%tune, '  ! kHz'
write (iu, *) '  vert_beta_freq   = ', freq%y%tune
write (iu, *) '/END'

if (what == 'ALL') then

  if (who == 'AC_ETA') then
    write (iu, *)
    write (iu, *) '&AC_ETA_PARAMETERS'
    write (iu, *) '  cesr_data_param%ac_z_amp_fit   = ', ac_eta%param%ac_z_amp_fit
    write (iu, *) '  cesr_data_param%ac_z_phase_fit = ', ac_eta%param%ac_z_phase_fit
    write (iu, *) '  cesr_data_param%chisq          = ', ac_eta%param%chisq
    write (iu, *) '/END'
  endif

  write (iu, *)    
  write (iu, *) '&DATA_BASE'
  do i = 1, db_num
    name = db_name(i)
    call to_node_name(name)
    ix = vnumbr(name)
    call vxgetn (name, 1, ix, ival)
    write (iu, '(3x, 2a)') db_name(i), ' = '
    write (iu, '((6x, 10i7))') (ival(j), j = 1, ix)
  enddo
  write (iu, *) '/END'
endif

write (iu, *)
                                               
end subroutine write_file_header

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

subroutine check_setup (skew_quad_on, species, who, abort_wanted)

implicit none

type (cesr_freq_struct) freq

logical H_tracker_locked, V_tracker_locked, track_2, species_ok, train_ok
logical fdbk_ok(8), abort_wanted
logical h_volt_stat        ! H_Separators off (warning if not)
logical error              ! .true. if any errors in setup
logical skew_error         ! .true. if skew quad on
logical :: skew_quad_on(0:)
logical check_v_tracker_lock ! ac_eta measurement does not use the vert tracker

integer i, ix, fdbk(8), fdbk_limit(8)
integer species, ipol
integer :: h_volts(12), h_volt_limit = 20 ! limit at 20 cu => 2kV
integer cu_skew_quad(98)

character(8) :: ans = 'Y'
character(*) who

! Check the control system to see if the conditions are right to begin a
! measurement. Eventually, will check:
!     * Horizontal Feedback Off
!     * Vertical Feedback Off
!     * Tune trackers locked/not locked as appropriate.
!     * Horiz and Vert Separators off.
!     * Buttons connected
!     * Current in CESR of correct (eg: there is positron current when ..._rp)
!     * Skew quads turned on
!
! For the logical variables which follow, True => OK (eg: fdbk off)

! Loop until everything is OK

check_v_tracker_lock = .true.
if (who == 'AC_ETA') check_v_tracker_lock = .false.

do 

  error = .false.

  ! Check feedback

  fdbk_limit = (/ 0, 0, 0, 0, 200, 200, 99999, 99999 /)
  fdbk_ok = .true.
  call vmgcmd('CSR FEED CON', 1, 8, fdbk)

  do i = 5, 8
    if (species == +1 .and. (i == 6 .or. i == 8)) cycle  ! positrons
    if (species == -1 .and. (i == 5 .or. i == 7)) cycle  ! electrons
    if (abs(fdbk(i)) > fdbk_limit(i)) then
      fdbk_ok(i) = .FALSE.
      error = .TRUE.
    endif
  enddo

  ! If shaking check that the tune tracker is locked.

  call get_tracker_lock_status (h_tracker_locked, v_tracker_locked)
  if (.not. h_tracker_locked) error = .true.

  if (check_v_tracker_lock) then
    if (.not. v_tracker_locked) error = .true.
  else
    v_tracker_locked = .true.
  endif

  ! Check the Horizontal Separators for voltage on the plates. Set
  ! h_volt_stat to false if any supplies exceed limit.

  call vxgetn('CSR HSP VVAL', 1, 12, h_volts)
  h_volt_stat = .true.
!!!!  do i = 1, 12   ! Reinstate this after undulator run....
  do i = 3, 12    ! Skip Sep 8w which is out for undulator run.
    if (h_volts(i) .gt. h_volt_limit) then
      h_volt_stat = .false.
      error = .true.
     endif
  enddo

  ! check skew quads

  skew_error = .false.
  call vxgetn ('CSR SQEWQUAD', 1, 98, cu_skew_quad)
  do i = 1, 98
    if (abs(cu_skew_quad(i)) > 2 .and. .not. skew_quad_on(i)) then
      skew_error = .true.
    endif
  enddo

  ! check that we have not locked both tune trackers to the same tune

  call get_freq (freq, 5)
  if (abs(freq%x%shake - freq%y%shake) < 1.0_rp) then
    track_2 = .false.
    error = .true.
  else
    track_2 = .true.
  endif

  ! Check detec species

  species_ok = .true.
  call vmgetn ('CSRBPM COMD ', 8, 8, ipol)
  if (ipol /= species) then
    species_ok = .false.
    error = .true.
  endif

  ! Check detec train number

  train_ok = .true.
  call vmgetn ('CSRBPM COMD ', 9, 9, ix)
  if (ix /= 1) then
    train_ok = .false.
    error = .true.
  endif

  ! Display the results of the check

  ! print *,' '
  ! print *,'         ****** Checking Setup *******'
  ! print *,' '

  ! First, type out the things that are right. . .

  ! if (fdbk_ok(5)) print *, 'OK: Horz e+ train feedback'
  ! if (fdbk_ok(6)) print *, 'OK: Horz e- train feedback'
  ! if (fdbk_ok(7)) print *, 'OK: Vert e+ train feedback'
  ! if (fdbk_ok(8)) print *, 'OK: Vert e- train feedback'
  ! if (h_tracker_locked) print *, 'OK: Horz Tracker locked'
  ! if (v_tracker_locked .and. check_v_tracker_lock) &
  !                      print *, 'OK: Vert Tracker locked'
  ! if (h_volt_stat)      print *, 'OK: Horiz separators off'
  ! if (.not. skew_error) print *, 'OK: Skew Quads off'
  ! if (species_ok)       print *, 'OK: DETEC species setting'
  ! if (train_ok)         print *, 'OK: DETEC train setting'

  ! print *

  ! . . .then the things that are wrong.

  if (.not. fdbk_ok(5)) print *, 'ERROR: HORZ E+ TRAIN FEEDBACK IS ON !'
  if (.not. fdbk_ok(6)) print *, 'ERROR: HORZ E- TRAIN FEEDBACK IS ON !'
  if (.not. fdbk_ok(7)) print *, 'ERROR: VERT E+ TRAIN FEEDBACK IS ON !'
  if (.not. fdbk_ok(8)) print *, 'ERROR: VERT E- TRAIN FEEDBACK IS ON !'
  if (.not. h_tracker_locked) print *, 'ERROR: HORZ TRACKER NOT LOCKED !!!!!!!!'
  if (.not. v_tracker_locked .and. check_v_tracker_lock) &
                     print *, 'ERROR: VERT TRACKER NOT LOCKED !!!!!!!!!'
  if (track_2 .and. .not. h_tracker_locked .and. .not. v_tracker_locked) &
                     print *, 'TO SWITCH TRACKER SET: LOGIC%DIGITAL_SHAKER_ON = T/F'
  if (.not. track_2) print *, 'ERROR: BOTH TRACKERS LOCKED ON TO THE SAME TUNE !!!!!!!!!!'
  if (.not. h_volt_stat) print *, 'WARNING: VOLTAGE ON HORIZ SEPARATORS !'
  !!! do i = 1, 98
  !!!   if (abs(cu_skew_quad(i)) > 2 .and. .not. skew_quad_on(i)) &
  !!!                 print *, 'WARNING: Skew Quad ON:', i, cu_skew_quad(i)
  !!! enddo
  if (.not. species_ok) print *, 'ERROR: DETEC SPECIES SETTING WRONG!!!!!!!!!!!!!!!!!'
  if (.not. train_ok) print *,   'ERROR: DETEC TRAIN SETTING WRONG!!!!!!!!!!!!!!!!!'

  abort_wanted = .false.
  if (.not. error) return

  if (logic%command_file_open) then
    call speech('ATTENTION! CESRV HAS DETECTED A PROBLEM. ASSISTANCE IS NECESSARY.')
  endif

  print *
  call strget('Do you want to check the setup again? Y,N,Abort: @', ans)
  call string_trim(ans, ans, ix)
  call str_upcase(ans, ans)
  if (ans(1:1) == 'A') then
    abort_wanted = .true.
    if (logic%command_file_open) then
      print *, 'TO BE SAFE I AM CLOSING THE COMMAND FILE...'
      close (logic%iu_command_file)
      logic%command_file_open = .false.
    endif
  endif
  if (ans(1:1) == 'N' .or. ans(1:1) == 'A') return

enddo  

end subroutine check_setup

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

subroutine get_freq (freq, n_read, slow_meas)

implicit none

type (cesr_freq_struct) freq

real(rp) f_avg(3), f_dev(3)
integer, optional :: n_read
logical, optional :: slow_meas
logical stat(3)


! Read the vertical and horizontal frequencies of the tune trackers

if (logic%digital_shaker_on) then
  call read_digital_tracker_frequencies (3, stat, f_avg, f_dev, integer_option(5, n_read))
  freq%x%shake = f_avg(1)
  freq%y%shake = f_avg(2)
else
  call read_analog_tracker_frequencies (freq%x%shake, freq%y%shake, &
                                  integer_option(1, n_read), logic_option(.false., slow_meas))
endif  

if (freq%x%reflection) then
  freq%x%tune = freq%rev - freq%x%shake
else
  freq%x%tune = freq%x%shake
endif

if (freq%y%reflection) then
  freq%y%tune = freq%rev - freq%y%shake
else
  freq%y%tune = freq%y%shake
endif

return

end subroutine get_freq

!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------

subroutine button_phase (det, amp_minn, hv)
    
implicit none

type (cesr_det_plane_struct) det

real(rp) sum, sum2, amp_min, target_phase, amp_minn, phase(4)
integer n_buts, but, hv


! flip phases

phase = det%but(:)%phase

if (hv == x_plane$) then
  phase(1) = phase(1) - 180
  phase(3) = phase(3) - 180
else
  phase(1) = phase(1) - 180
  phase(2) = phase(2) - 180
endif

! find which buttons are ok
! Phases of buttons can be different by multiples of 360 deg.
! This is a problem when we take an average.
! Find a OK button and move other phases accordingly by multiples of 360.

if (amp_minn .ne. 0) then
  amp_min = amp_minn
else
  amp_min = amp_min_default
endif

! for the new system assume all the buttons are ok

if (det%system_id > 0) then  ! for new system
  det%but%ok = .true.
  target_phase = phase(1)
else
  do but = 1, 4
    if (det%but(but)%amp .gt. amp_min) then
      det%but(but)%ok = .true.
      target_phase = phase(but)
    else
      det%but(but)%ok = .false.
    endif
  enddo
endif

! In any case, button is not ok if the value is NaN

do but = 1, 4
  if (isnan(det%but(but)%amp) .or. isnan(det%but(but)%phase)) then
    det%but(but)%ok = .false.
    det%but(but)%amp = 0
    det%but(but)%phase = 666.0 ! Garbage number
  endif
enddo

! Adjust phases of the buttons if off by 360 deg

do but = 1, 4
  if(det%but(but)%ok)then
    phase(but) = phase(but) + 360 *  nint((target_phase - phase(but)) / 360)
  endif
enddo

! compute average and rms

n_buts = 0
sum = 0
sum2 = 0

do but = 1, 4
  if (det%but(but)%ok) then
    n_buts = n_buts + 1
    sum  = sum  + phase(but)
    sum2 = sum2 + phase(but)**2
  endif
enddo

if (n_buts > 1) then
  det%phase_meas = sum / n_buts
  det%rms_phase_meas = sqrt(max((sum2/n_buts - (sum/n_buts)**2), 0.0_rp))
endif

det%n_buts = n_buts                         

! for the new system we already know if the data is good.

if (det%system_id == 0) det%ok = (n_buts > 1)

end subroutine button_phase

!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
! Note: This routine is not being used. DCS 2009-12-22

subroutine xy_shake_components (plane)

implicit none

type (cesr_det_plane_struct) plane(0:)
!  type (cesr_det_dc_position_struct) dc(0:)

real(rp) x_real, y_real, x_imag, y_imag, re(4), im(4)
integer det, i

! Turn off denominator sumation and zero offset in poscon.

call err_exit  ! And make sure no one uses this routine.

! offset off
call poscon2 (-1, 0._rp, 0._rp, 0._rp, 0._rp, 0._rp, 0._rp)
! denominator off
call poscon2 (-3, 0._rp, 0._rp, 0._rp, 0._rp, y_real, x_real)

! Compute the signal amplitude in the x and y planes
     
do det = 0, ubound(plane, 1)

  do i = 1, 4
    re(i) = plane(det)%but(i)%amp * cosd(plane(det)%but(i)%phase)
    im(i) = plane(det)%but(i)%amp * sind(plane(det)%but(i)%phase)
  enddo

  call poscon2 (det, re(1), re(2), re(3), re(4), y_real, x_real)
  call poscon2 (det, im(1), im(2), im(3), im(4), y_imag, x_imag)

  plane(det)%x%amp = sqrt (x_real ** 2 + x_imag ** 2)
  plane(det)%y%amp = sqrt (y_real ** 2 + y_imag ** 2)

  if (plane(det)%x%amp /= 0 .and. plane(det)%y%amp /= 0) then
    plane(det)%x%phase = atan2d (x_imag, x_real)
    plane(det)%y%phase = atan2d (y_imag, y_real)
  else
    plane(det)%x%phase = 0
    plane(det)%y%phase = 0
  endif
  
enddo    

call poscon2 (-2, 0._rp, 0._rp, 0._rp, 0._rp, 0._rp, 0._rp)       ! offset off
call poscon2 (-4, 0._rp, 0._rp, 0._rp, 0._rp, y_real, x_real)     ! denominator off

end subroutine xy_shake_components

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

subroutine phase_shift_near_design (plane)

implicit none

type (cesr_det_plane_struct) plane(0:n_det_maxx)

real(rp) offset, design(0:n_det_maxx), meas(0:n_det_maxx)
real(rp) rms, rms1

integer n

! Shift the measured phases until all of them are within +/- 180 deg of design.
! Must take into account that there is an overall arbitrary offset between
! the design and the measured. 

n = count(plane(:)%ok)
design = plane(:)%phase_design
meas = plane(:)%phase_meas

do 
  offset = sum(meas - design, mask = plane(:)%ok) / n
  if (all(meas - design - offset < 180)) exit
  meas = meas + 360 *  nint((design + offset - meas) / 360)
enddo

rms = sum((meas - design - offset)**2)
plane(:)%phase_meas = meas

! See if shifting the offset improves anything

offset = offset + 90
meas = meas + 360 *  nint((design + offset - meas) / 360)
rms1 = sum((meas - design - offset)**2)
if (rms1 < rms) then
  plane(:)%phase_meas = meas
  rms = rms1
endif

end subroutine phase_shift_near_design
     
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

subroutine poscon2 (det, b1, b2, b3, b4, r1, r2)

implicit none

real(rp) b1, b2, b3, b4, r1, r2
real rr1, rr2
integer det, d

!      

d = det
if (d == 0) d = 100
call poscon (d, real(b1), real(b2), real(b3), real(b4), rr1, rr2)
r1 = rr1
r2 = rr2

end subroutine poscon2

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

subroutine analyze_ac_eta_data (horz, freq, dc, p, ac_eta, u, do_walski)

implicit none

type (cesr_freq_struct), target :: freq
type (cesr_phase_params_struct) :: p
type (cesr_det_plane_struct), target :: horz(0:n_det_maxx)
type (cesr_det_dc_position_struct) :: dc(0:n_det_maxx)
type (universe_struct) u
type (cesr_all_data_struct) ac_eta, ac_eta_CBPM, ac_eta_old
type (db_struct) db

real(rp) vec(2)

integer det
integer i
integer ix_data_set

logical do_walski

!

do det = 0, ubound(horz, 1)
  if (horz(det)%system_id > 0) cycle  ! Already set
  horz(det)%ok = all(horz(det)%but(:)%amp > amp_min_default)
enddo

do det = 0, 120
   dc(det)%signal = 1000
enddo

if(logic%nonlinear_calc)then
  call nonlin_xy_shake_components (horz, dc)
 else
  call test_amp(horz, dc)
endif

if (logic%wolski_normal_mode_calc_on .and. do_walski) then
  do i = lbound(horz, 1), ubound(horz, 1)
    vec = matmul (u%but_to_mode(i)%mat, horz(i)%but(:)%amp)
    horz(i)%x%amp = vec(1)
    horz(i)%y%amp = vec(2)
  enddo
endif

call ac_eta_calc (u, horz, freq, ac_eta_CBPM)

do i = lbound(db%detector, 1), ubound(db%detector, 1)
  if(horz(i)%system_id < 1) cycle  ! If not new system or does not exist
  ac_eta%ac_etap_x(i)%value = ac_eta_CBPM%ac_etap_x(i)%value
  ac_eta%ac_eta_x(i)%value = ac_eta_CBPM%ac_eta_x(i)%value
  ac_eta%ac_etap_y(i)%value = ac_eta_CBPM%ac_etap_y(i)%value
  ac_eta%ac_eta_y(i)%value = ac_eta_CBPM%ac_eta_y(i)%value

  ac_eta%cbar12_a(i)%value = ac_eta_CBPM%cbar12_a(i)%value !V16
  ac_eta%cbar12_b(i)%value = ac_eta_CBPM%cbar12_b(i)%value !V36
  ac_eta%yxcos(i)%value = ac_eta_CBPM%yxcos(i)%value !yamp/xamp * cos(phi_y-phi_x)
  ac_eta%yxsin(i)%value = ac_eta_CBPM%yxsin(i)%value !yamp/xamp * sin(phi_y-phi_x)
end do

end subroutine analyze_ac_eta_data

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

subroutine analyze_ac_eta_data2 (horz, freq, dc, p, ac_eta, u, do_walski)

implicit none

type (cesr_freq_struct), target :: freq
type (cesr_phase_params_struct) :: p
type (cesr_det_plane_struct), target :: horz(0:n_det_maxx)
type (cesr_det_dc_position_struct) :: dc(0:n_det_maxx)
type (universe_struct) u
type (cesr_all_data_struct) ac_eta, ac_eta_CBPM, ac_eta_old
type (db_struct) db

real(rp) vec(2)

integer det
integer i
integer ix_data_set

logical do_walski

!

do det = 0, ubound(horz, 1)
  if (horz(det)%system_id > 0) cycle  ! Already set
  horz(det)%ok = all(horz(det)%but(:)%amp > amp_min_default)
enddo

do det = 0, 120
   dc(det)%signal = 1000
enddo

if(logic%nonlinear_calc)then
  call nonlin_xy_shake_components (horz, dc)
 else
  call test_amp(horz, dc)
endif

if (logic%wolski_normal_mode_calc_on .and. do_walski) then
  do i = lbound(horz, 1), ubound(horz, 1)
    vec = matmul (u%but_to_mode(i)%mat, horz(i)%but(:)%amp)
    horz(i)%x%amp = vec(1)
    horz(i)%y%amp = vec(2)
  enddo
endif

call ac_eta_calc2 (u, horz, freq, ac_eta_CBPM)

do i = lbound(db%detector, 1), ubound(db%detector, 1)
  if(horz(i)%system_id < 1) cycle  ! If not new system or does not exist
  ac_eta%ac_etap_x(i)%value = ac_eta_CBPM%ac_etap_x(i)%value
  ac_eta%ac_eta_x(i)%value = ac_eta_CBPM%ac_eta_x(i)%value
  ac_eta%ac_etap_y(i)%value = ac_eta_CBPM%ac_etap_y(i)%value
  ac_eta%ac_eta_y(i)%value = ac_eta_CBPM%ac_eta_y(i)%value
end do

end subroutine analyze_ac_eta_data2

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

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

implicit none

type (cesr_freq_struct), target :: freq
type (a_cesr_freq_struct), pointer :: a_freq, b_freq
type (cesr_phase_params_struct) :: p
type (cesr_det_plane_struct), target :: horz(0:n_det_maxx), vert(0:n_det_maxx)
type (cesr_det_plane_struct), pointer :: a_mode, b_mode
type (cesr_det_dc_position_struct) :: dc(0:n_det_maxx)
type (lat_struct), pointer :: ring
type (twiss_struct), pointer :: h_mode, v_mode
type (universe_struct), target :: u

real(rp) rel_amp, rel_phase, beta_ratio
real(rp) tune_x, tune_y
real(rp) beta_a, beta_b, cbar_max
real(rp) a_x_design_amp, a_x_design_phase, a_y_design_amp, a_y_design_phase
real(rp) a_y_in_phase_amp, a_y_out_phase_amp
real(rp) b_x_design_amp, b_x_design_phase, b_y_design_amp, b_y_design_phase
real(rp) b_x_in_phase_amp, b_x_out_phase_amp
real(rp) x_phase_net, y_phase_net
real(rp) x_wgt, y_wgt, amp_ratio, gamma_c 
real(rp) cbar(2,2)
real(rp) horz_cos(1:4), horz_sin(1:4), vert_cos(1:4), vert_sin(1:4), ccos, ssin, deg_to_rad
real(rp) h_scale, v_scale
real(rp) betah(120), betav(120)

integer i, idet, ix, ir

!---------------------------------------------------------------
! calculate the Cbar's and phases
! NOTE: Since cos(phi + phi0) = cos(phi_0) cos(phi) - sin(phi_0) sin(phi)
! there is an extra minus sign for the cbar12 (sin like) component:
! x = x_amp * cos(omega*t)
! y = y_amp * cos(omega*t + del_phase)   ! del_phase = y_phase - x_phase
!   = y_amp * (cos(omega*t) * cos(del_phase) - sin(omega*t) * sin(del_phase))
!   = -cbar22 * cos(omega*t) - cbar12 * sin(omega*t)

! first some initial calculations

ring => u%design

call nonlin_xy_shake_components (horz, dc)
call nonlin_xy_shake_components (vert, dc)

do i = 0, 120
  ix = u%db%detector(i)%ix_lat
  if (ix /= 0) then
    if (ring%ele(ix)%mode_flip) then
      h_mode => ring%ele(ix)%b
      v_mode => ring%ele(ix)%a
    else
      h_mode => ring%ele(ix)%a
      v_mode => ring%ele(ix)%b
    endif
    horz(i)%beta_design = h_mode%beta
    vert(i)%beta_design = v_mode%beta
    horz(i)%phase_design = h_mode%phi * 180 / pi
    vert(i)%phase_design = v_mode%phi * 180 / pi
  endif
enddo
     
do idet = 0, 120

  ir = u%db%detector(idet)%ix_lat
  if (ir == 0) then
    horz(idet)%cbar12 = 0
    horz(idet)%cbar22 = 0
    vert(idet)%cbar12 = 0
    vert(idet)%cbar11 = 0

  else
    call c_to_cbar (ring%ele(ir), cbar)

    !-----------------------------------------------------------------
    ! Mobius Cbar and phases

    cbar_max = max(abs(cbar(1,1)), abs(cbar(1,2)), abs(cbar(2,2)))
    if (ring%ele(ir)%mode_flip .or. cbar_max .gt. 0.2_rp) then

      gamma_c = ring%ele(ir)%gamma_c

      horz(idet)%ok = .true.
      horz(idet)%but%ok = .true.
      horz(idet)%n_buts = 4
      vert(idet)%ok = .true.
      vert(idet)%but%ok = .true.
      vert(idet)%n_buts = 4

      beta_a = ring%ele(ir)%a%beta 
      beta_b = ring%ele(ir)%b%beta

      if (ring%ele(ir)%mode_flip) then
        a_mode => vert(idet)
        a_freq => freq%y
        b_mode => horz(idet)
        b_freq => freq%x
      else
        a_mode => horz(idet)
        a_freq => freq%x
        b_mode => vert(idet)
        b_freq => freq%y
      endif

      ! a mode anal

      beta_ratio = beta_a / beta_b
      a_x_design_amp = gamma_c * sqrt(beta_a)
      a_x_design_phase = 0
      a_y_in_phase_amp  = -cbar(2,2) * sqrt(beta_b)
      a_y_out_phase_amp = -cbar(1,2) * sqrt(beta_b)
      a_y_design_amp = sqrt (a_y_in_phase_amp**2 + a_y_out_phase_amp**2)
      a_y_design_phase = -atan2(a_y_out_phase_amp, a_y_in_phase_amp)
      if (a_freq%reflection) a_y_design_phase = -a_y_design_phase

      x_phase_net = a_mode%x%phase - a_x_design_phase * 180 / pi
      y_phase_net = a_mode%y%phase - a_y_design_phase * 180 / pi

      ! adjust by factors of 360

      if (a_x_design_amp .gt. a_y_design_amp) then
        y_phase_net = y_phase_net - nint((y_phase_net - x_phase_net)/360) * 360
      else
        x_phase_net = x_phase_net - nint((x_phase_net - y_phase_net)/360) * 360
      endif

      !

      if (a_mode%x%amp == 0) then
        a_mode%phase_meas = 0
        a_mode%cbar12 = 0
        horz(idet)%cbar22 = 0

      else
        x_wgt = (a_x_design_amp * a_mode%x%amp)
        y_wgt = (a_y_design_amp * a_mode%y%amp)
        a_mode%phase_meas = (x_phase_net * x_wgt + y_phase_net * y_wgt) / (x_wgt + y_wgt)

        rel_phase = a_mode%y%phase - a_mode%x%phase
        if (freq%x%reflection) rel_phase = -rel_phase
        if (p%species == -1) rel_phase = -rel_phase
      
        amp_ratio = a_mode%y%amp / a_mode%x%amp
        a_mode%cbar12 = amp_ratio * sqrt(beta_ratio) * sind(rel_phase) * gamma_c
        horz(idet)%cbar22 = -amp_ratio * sqrt(beta_ratio) * cosd(rel_phase) * gamma_c
      endif
               
      ! b mode anal

      beta_ratio = beta_b / beta_a
      b_y_design_amp = gamma_c * sqrt(beta_b)          
      b_y_design_phase = 0
      b_x_in_phase_amp  =  cbar(1,1) * sqrt(beta_a) 
      b_x_out_phase_amp = -cbar(1,2) * sqrt(beta_a)
      b_x_design_amp = sqrt (b_x_in_phase_amp**2 + b_x_out_phase_amp**2)
      b_x_design_phase = -atan2(b_x_out_phase_amp, b_x_in_phase_amp)
      if (b_freq%reflection) b_x_design_phase = -b_x_design_phase

      x_phase_net = b_mode%x%phase - b_x_design_phase * 180 / pi
      y_phase_net = b_mode%y%phase - b_y_design_phase * 180 / pi

      ! adjust by factors of twopi

      if (b_x_design_amp .gt. b_y_design_amp) then
        y_phase_net = y_phase_net - nint((y_phase_net-x_phase_net)/360) * 360
      else
        x_phase_net = x_phase_net - nint((x_phase_net-y_phase_net)/360) * 360
      endif

      !

      if (b_mode%y%amp == 0) then
        b_mode%phase_meas = 0
        b_mode%cbar12 = 0
        vert(idet)%cbar11 = 0

      else

        x_wgt = (b_x_design_amp * b_mode%x%amp)
        y_wgt = (b_y_design_amp * b_mode%y%amp)
        b_mode%phase_meas = (x_phase_net * x_wgt + y_phase_net * y_wgt) / (x_wgt + y_wgt)

        rel_phase = b_mode%x%phase - b_mode%y%phase
        if (freq%x%reflection) rel_phase = -rel_phase
        if (p%species == -1) rel_phase = -rel_phase
        amp_ratio = b_mode%x%amp / b_mode%y%amp
        b_mode%cbar12 =  -amp_ratio * sqrt(beta_ratio) * sind(rel_phase) * gamma_c
        vert(idet)%cbar11 =  amp_ratio * sqrt(beta_ratio) * cosd(rel_phase) * gamma_c
      endif
        
    !---------------------------------------------------------
    ! Non-Mobius Cbar and phases

    else  

      if (horz(idet)%x%amp /= 0) then
        rel_amp = horz(idet)%y%amp / horz(idet)%x%amp
        rel_phase = horz(idet)%y%phase - horz(idet)%x%phase
        if (freq%x%reflection) rel_phase = -rel_phase
        if (p%species == -1) rel_phase = -rel_phase
        horz(idet)%cbar22 = -rel_amp * cosd(rel_phase)
        horz(idet)%cbar12 = rel_amp * sind(rel_phase)
      else
        horz(idet)%cbar22 = 0
        horz(idet)%cbar12 = 0
      endif

      if (vert(idet)%y%amp /= 0) then
        rel_amp = vert(idet)%x%amp / vert(idet)%y%amp
        rel_phase = vert(idet)%x%phase - vert(idet)%y%phase
        if (freq%y%reflection) rel_phase = -rel_phase
        if (p%species == -1) rel_phase = -rel_phase
        vert(idet)%cbar11 = rel_amp * cosd(rel_phase)
        vert(idet)%cbar12 = rel_amp * sind(rel_phase)
      else
        vert(idet)%cbar11 = 0
        vert(idet)%cbar12 = 0
      endif

      ! get the throry beta's and convert C's to Cbar's

      beta_ratio = vert(idet)%beta_design / horz(idet)%beta_design
      horz(idet)%cbar12 = horz(idet)%cbar12 / sqrt(beta_ratio)
      horz(idet)%cbar22 = horz(idet)%cbar22 / sqrt(beta_ratio)
      vert(idet)%cbar12 = vert(idet)%cbar12 * sqrt(beta_ratio)
      vert(idet)%cbar11 = vert(idet)%cbar11 * sqrt(beta_ratio)

      call button_phase (horz(idet), 0.0_rp, x_plane$)
      call button_phase (vert(idet), 0.0_rp, y_plane$)
      ! nonlinear correction
      horz(idet)%phase_meas = horz(idet)%x%phase
      vert(idet)%phase_meas = vert(idet)%y%phase


    endif ! non-Mobius cbar and phases

  endif ! idet exists

enddo ! idet

!---------------------------------------------------------------
! limit cbar range so values do not overflow the output file

do idet = 0, 120
  horz(idet)%cbar12 = max(min(horz(idet)%cbar12, 9.0_rp), -9.0_rp)
  horz(idet)%cbar22 = max(min(horz(idet)%cbar22, 9.0_rp), -9.0_rp)
  vert(idet)%cbar12 = max(min(vert(idet)%cbar12, 9.0_rp), -9.0_rp)
  vert(idet)%cbar11 = max(min(vert(idet)%cbar11, 9.0_rp), -9.0_rp)
enddo

! reverse the phases for electron data

tune_x = ring%ele(ring%n_ele_track)%a%phi * 180 / pi
tune_y = ring%ele(ring%n_ele_track)%b%phi * 180 / pi

if (p%species == -1) then
  do idet = 0, 120
    horz(idet)%phase_design = tune_x - horz(idet)%phase_design
    vert(idet)%phase_design = tune_y - vert(idet)%phase_design
  enddo
endif

! compute x-y amplitudes for measurement of beta

deg_to_rad = twopi/360.
i=0
do idet = 0,120 !cbpm II
  if(horz(idet)%system_id < 1)cycle    
  ir = u%db%detector(idet)%ix_lat
  if(ir == 0)cycle
  if(any(horz(idet)%but(1:4)%amp < 1.e-10))cycle 
  i=i+1
  horz_cos(1:4) = horz(idet)%but(1:4)%amp * cos(horz(idet)%but(1:4)%phase * deg_to_rad)
  horz_sin(1:4) = horz(idet)%but(1:4)%amp * sin(horz(idet)%but(1:4)%phase * deg_to_rad)
  ccos = horz_cos(4)-horz_cos(3)+horz_cos(2)-horz_cos(1)
  ssin = horz_sin(4)-horz_sin(3)+horz_sin(2)-horz_sin(1)
  horz(idet)%beta_meas = (ccos**2 + ssin**2)
  betah(i) = horz(idet)%beta_meas/ring%ele(ir)%a%beta 

  vert_cos(1:4) = vert(idet)%but(1:4)%amp * cos(vert(idet)%but(1:4)%phase * deg_to_rad)
  vert_sin(1:4) = vert(idet)%but(1:4)%amp * sin(vert(idet)%but(1:4)%phase * deg_to_rad)
  ccos = vert_cos(4)+vert_cos(3)-vert_cos(2)-vert_cos(1)
  ssin = vert_sin(4)+vert_sin(3)-vert_sin(2)-vert_sin(1)
  vert(idet)%beta_meas = (ccos**2 + ssin**2)
  betav(i) = vert(idet)%beta_meas/ring%ele(ir)%b%beta

end do
call sort(betah(1:i))
h_scale = 1/betah(i/2)
horz(0:120)%beta_meas = horz(0:120)%beta_meas * h_scale

call sort(betav(1:i))
v_scale = 1/betav(i/2)
vert(0:120)%beta_meas = vert(0:120)%beta_meas * v_scale

!  open(unit=11)
!  do idet = 0,120
!    ir = u%db%detector(idet)%ix_lat
!    print '(1x,i,4e12.4)', idet, horz(idet)%beta_meas, vert(idet)%beta_meas, &
!                        ring%ele(ir)%a%beta, ring%ele(ir)%b%beta 
!    write(11, '(1x,i,4e12.4)') idet, horz(idet)%beta_meas, vert(idet)%beta_meas, &
!                        ring%ele(ir)%a%beta, ring%ele(ir)%b%beta 
!  end do
!  close(unit=11)

! Correct for:  
!     Sign if shaking at reflected.

do idet = 0, 120
  if (freq%x%reflection) horz(idet)%phase_meas = -horz(idet)%phase_meas
  if (freq%y%reflection) vert(idet)%phase_meas = -vert(idet)%phase_meas
enddo

! To be able to compare design with experiment we arbitrarily shift all the
! measured phases so that there is no phase difference at the first OK
! detector that has an rms less than 1.0
! Also: Phase can be off by multiples of 360 so shift phase to
! roughly correspond to the theoretical values.

call phase_shift_near_design (horz)
call phase_shift_near_design (vert)

end subroutine analyze_phase_data

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

subroutine write_shaking_data (horz, vert, freq, dc, p, ac_eta, who)

implicit none

type (cesr_freq_struct), target :: freq
type (cesr_phase_params_struct) p
type (cesr_det_plane_struct), target :: horz(0:n_det_maxx), vert(0:n_det_maxx)
type (cesr_det_dc_position_struct) dc(0:n_det_maxx)
type (cesr_all_data_struct) ac_eta

real(rp) phase(4)
integer i, det

character(*) who
character(40) fmt1
character(48) line

!----------------------------------------------------------------------------
! First write the raw data.

call write_file_header ('ALL', freq, p, ac_eta, who)

write (p%unit, *)
write (p%unit, '(10x, a)') 'Phase_1 Phase_2 Phase_3 Phase_4 | ' //  &
                             '    Amp_1       Amp_2       Amp_3       Amp_4'
write (p%unit, '(10x, a)') '             (deg)              | ' 
write (p%unit, '(a)') '&RAWDATA'
fmt1 = '(1x, a, i3, a, 4F8.2, a48, i4)'
do det = 0, n_det_maxx
  line = ''
  do i = 1, 4
    call amp_to_string (horz(det)%but(i)%amp, line(i*12-11:))
  enddo
  phase = max(min(999.99_rp, horz(det)%but(:)%phase), -999.99_rp)
  write (p%unit, fmt1) 'H_(', det, ')=', phase, line, horz(det)%system_id
  if (who /= 'PHASE') cycle
  line = ''
  do i = 1, 4
    call amp_to_string (vert(det)%but(i)%amp, line(i*12-11:))
  enddo
  phase = max(min(999.99_rp, vert(det)%but(:)%phase), -999.99_rp)
  write (p%unit, fmt1) 'V_(', det, ')=', phase, line, vert(det)%system_id
enddo
write (p%unit, '(a)') '/END'
write (p%unit, *) 


write (p%unit, '(9x, a)') &
        'det    x_orbit  y_orbit     But_1     But_2     But_3     But_4'
write (p%unit, *) '&RAWORBIT'
do det = 0, n_det_maxx
  if (any(dc(det)%signal < 0) .or. any(dc(det)%signal > 1e9)) then
    print *, 'ERROR: BPM BUTTON AMPLITUDE OUT OF RANGE FOR BPM:', det
    print *, '       SETTING ALL AMPLITUDES TO ZERO'
    dc(det)%signal = 0
    dc(det)%x = 0; dc(det)%y = 0
  endif
  write (p%unit, '(2x, a, i3, a, 2f9.2, 4i10)') 'ORBIT_(', det, ')=', &
                                 dc(det)%x, dc(det)%y,  &
                                 (nint(dc(det)%signal(i)), i = 1, 4)
enddo
write (p%unit, *) '/END'
write (p%unit, *)

! AC_ETA data

if (who == 'AC_ETA') then

  write (p%unit, '(a)') &
  '            |      Horizontal   |        Vertical     |', &
  '   Det      |    Eta       Etap |      Eta       Etap | OK'

  write (p%unit, '(a)') '&AC_ETA_DATA'
  do det = 0, 120
    write (p%unit, '(2x, a, i3, a, 4(f11.6), 1x, l3)') &
        'AC_ETA(', det, ') =', &
        ac_eta%ac_eta_x(det)%value, ac_eta%ac_etap_x(det)%value, &
        ac_eta%ac_eta_y(det)%value, ac_eta%ac_etap_y(det)%value, &
        horz(det)%ok
  enddo

  write (p%unit, '(a)') '/END'

! write phase and cbar data

else

  write (p%unit, '(a)') &
    '           |       Horizontal         |         Vertical         | OK', &
    '   Det     |  Phase   Cbar22   Cbar12 |  Phase   Cbar12   Cbar11 | x  y', &
    '           |  (deg)                   |  (deg)                   |'

  write (p%unit, '(a)') '&PHASE_CBAR_DATA'
  do det = 0, 120
    write (p%unit, '(2x, a, i3, a, 2(f9.2, 2f9.5), 1x, 2l3)') &
        'PC_(', det, ') =', &
        horz(det)%phase_meas, horz(det)%cbar22, horz(det)%cbar12,  &
        vert(det)%phase_meas, vert(det)%cbar12, vert(det)%cbar11,  &
        horz(det)%ok, vert(det)%ok     
  enddo

  write (p%unit, '(a)') '/END'

  ! beta data

  write (p%unit, *)
  write (p%unit, '(2a)') '   Det         |    beta_a      beta_b  OK_a  OK_b '
  
  write (p%unit, '(a)') '&BETA_DATA'
  do i = 0, 120
    write (p%unit, '(2x, a, i3, a, 2es12.4, 2l6)') &
        'beta_meas(', i, ') =', &
        horz(i)%beta_meas, vert(i)%beta_meas, horz(i)%ok, vert(i)%ok
  enddo

  write (p%unit, '(a)') '/END'

  ! raw amp/phase data

  write (p%unit, *)
  write (p%unit, '(2a)') &
    '               |                      a-mode                  ', &
                   '|                      b-mode                  ', &
    '   Det         |     x_amp   x_phase       y_amp   y_phase OK ', &
                   '|     x_amp   x_phase       y_amp   y_phase OK '
  
  write (p%unit, '(a)') '&SHAKE_DATA'
  do i = 0, 120
    write (p%unit, '(2x, a, i3, a, 2(2(es12.4, f10.3), l3))') &
        'shake(', i, ') =', &
        horz(i)%x%amp, horz(i)%x%phase, horz(i)%y%amp, horz(i)%y%phase, horz(i)%ok, &
        vert(i)%x%amp, vert(i)%x%phase, vert(i)%y%amp, vert(i)%y%phase, vert(i)%ok     
  enddo

  write (p%unit, '(a)') '/END'

endif

close (unit = p%unit)

end subroutine write_shaking_data
                       
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

subroutine amp_to_string (amp, str)

implicit none

real(rp) amp
character(*) str

!

if (amp < 2e9) then
  write (str, '(i12)') nint(amp)
else
  write (str, '(es12.3)') amp
endif

end subroutine amp_to_string

end module