module sim_ac_and_cbar_mod

  use bmad
  use sim_bpm_mod
  use ele_noise_mod
  use mode3_mod
  use tt_tools_mod
  use bsim_cesr_interface
  use nonlin_bpm_mod

contains

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Subroutine sim_cbar_data(ring, bpm, n_turns, track6x6, n_TTs, tt_params, &
!                     n_damping, init_vec, phi, cbar, A, ele_noise)
!
! Simulate a phase/coupling measurement. Lattice and BPM errors 
! are taken into account.
!
! Input:
!   ring  -- lat_struct
!   bpm(:) -- det_struct(:), holds information about a single BPM
!   n_turns -- integer, total number of turns to sample in the measurement
!   track6x6 -- logical. If .true., use full 6x6 If .false., offset particle 
!                 and track as 4x4.
!   n_TTs    -- number of tune trackers. if zero, no damping/excitation used
!   tt_params -- type(tt_param_struct), for tune trackers.
!   n_damping -- integer. If damping = .true., how many turns to damp before
!                recording orbit data?
!   init_vec -- real(rp)(1:6), in meters. How far to offset 
!                     particles before tracking. 
!   ele_noise(:) -- ele_noise_struct(:), optional, noise to apply to elements
!                   during tracking
!
! Output:
!   phi(:,2,2)  -- real(rp). First index runs over all BPMs. Second and 
!                  third indices represent phase measurements as:
!                  ((ax,ay),(bx,by)), where first letter is the mode
!                  and second letter is the plane of motion.
!   cbar(:,2,2) -- real(rp). First index runs over all BPMs. Second and 
!                  third indices are "measured" cbar matrix.
!   A(:,2,2)    -- real(rp). First index runs over all BPMs. Second and 
!                  third indices are amplitudes of the oscillations, 
!                  represented as: 
!                  ((Ax,Ay),(Bx,By))
!                  where first letter is the mode (a or b), and the second 
!                  is which plane the motion is in (x or y). 
!--------------------------------------------------------------------------
  subroutine sim_cbar_data(ring, bpm, n_turns, track6x6, n_TTs, tt_params, id, n_damping, init_vec, phi, cbar, A, ele_noise)

  ! Note: the returned phases will have a constant offset WRT phases calculated by Bmad. 
  ! This can be corrected by adding a constant shift to this "measured" data such that the average difference
  ! between Bmad and "measured" is zero.

    implicit none
    type(lat_struct) :: ring
    type(det_struct) :: bpm(:)
    type(det_data_struct), allocatable :: bpm_save(:,:)
    type(tt_param_struct) :: tt_params(:)
    type(ele_noise_struct), optional :: ele_noise(:)
    integer, intent(in) :: n_TTs, id(:)
    
    logical, intent(in) :: track6x6
    integer, intent(in) :: n_turns
    real(rp), intent(in) :: init_vec(6)
    integer, intent(in) :: n_damping
    real(rp) :: harvest
    real(rp), allocatable :: phi(:,:,:), cbar(:,:,:), A(:,:,:)
    
    ! 2011.05.25 - incorporate Suntao's button-based phase/amplitude notation:
    type (cesr_det_plane_struct) :: horz(0:n_det_maxx), vert(0:n_det_maxx)
    type (cesr_det_dc_position_struct) :: dc(0:n_det_maxx)
    
    integer :: i_turn = 0, i, ix, jx, n_bpms, idet, ix_db, ix_lat
    real(rp) rel_amp, rel_phase
    logical err
    real(rp), allocatable :: tt_tunes(:,:)
       
 
    n_bpms = ubound(bpm,1)
    
    if (.not. allocated(bpm_save)) allocate(bpm_save(n_turns,n_bpms))

    if (n_TTs .gt. 0) then
       if (.not. allocated(tt_tunes)) allocate(tt_tunes(n_turns,n_TTs))
    endif
    
    phi(:,:,:) = 0 ! = (n_bpms,2,2) = (bpm%ix_db,(ax,ay),(bx,by))
    A(:,:,:) = 0   ! = (n_bpms,2,2) = (bpm%ix_db,(ax,ay),(bx,by))
    cbar(:,:,:) = 0! = (n_bpms,2,2) = (bpm%ix_db,(11,12),(21,22))
        
    if (present(ele_noise)) then
       call sim_tbt_data(ring, bpm, n_turns, track6x6, n_TTs, tt_params, id, n_damping, init_vec, bpm_save, tt_tunes, ele_noise)
    else
       call sim_tbt_data(ring, bpm, n_turns, track6x6, n_TTs, tt_params, id, n_damping, init_vec, bpm_save, tt_tunes)
    end if

    call cal_butn_amp_phase(ring,bpm_save,tt_tunes, horz, vert) 

    write(*,*) " "
    write(*,*) "========================================================================"
    write(*,*) "RMS (measured - Bmad reported):"

    call calc_phase_cbar(ring, bpm_save, horz, vert, phi, cbar, A)

    write(*,*) "========================================================================"
    write(*,*) " " 

  end subroutine sim_cbar_data
  

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!Subroutine sim_ac_eta_data(ring, bpm, n_turns, n_TTs, tt_params, id, n_damping, init_vec, use_new_method, ac_eta, ele_noise)
! Simulate a ac_eta (dispersion) measurement. Lattice and BPM errors 
! are taken into account.
!
! Input:
!   ring  -- lat_struct
!   bpm(:) -- det_struct(:), holds information about a single BPM
!   n_turns -- integer, total number of turns to sample in the measurement
!   n_TTs    -- number of tune trackers. if zero, no damping/excitation used
!   tt_params -- type(tt_param_struct), for tune trackers.
!   n_damping -- integer. If damping = .true., how many turns to damp 
!                recording orbit data? 
!   init_vec -- real(rp)(1:6), in meters. How far to offset 
!                     particles before tracking. 
!   use_new_method --logical, .true. to use new method to calculate ac_eta
!                             .false. to use old method
!   ele_noise(:) -- ele_noise_struct(:), optional, noise to apply to elements
!                   when tracking
!
! Output:
!   ac_eta   -- cesr_all_data_struct, calculated ac disperion
!--------------------------------------------------------------------------
!
Subroutine sim_ac_eta_data(ring, bpm, n_turns, n_TTs, tt_params, id, n_damping, init_vec, use_new_method, ac_eta, ele_noise)

    implicit none
    type(lat_struct) :: ring
    type(det_struct) :: bpm(:)
    integer :: n_turns
    integer :: n_TTs
    type(tt_param_struct) :: tt_params(:)
    integer :: id(:)
    integer :: n_damping
    real(rp) :: init_vec(6)
    real(rp) :: harvest
    logical :: use_new_method
    type (cesr_all_data_struct) :: ac_eta 
    type(ele_noise_struct), optional :: ele_noise(:)

    type (det_data_struct), allocatable :: bpm_save(:,:)
    real(rp), allocatable :: tt_tunes(:,:)
    integer n_bpms,i, ix, ix_db
    real(rp) zFracTune

! parameters for ac_eta calculation
    type (cesr_det_dc_position_struct) dc(0:n_det_maxx)
    type (cesr_freq_struct), target :: freq
    type (cesr_det_plane_struct) :: long(0:n_det_maxx)
    real(rp) ,allocatable :: phiz(:), Az(:)

  n_bpms = ubound(bpm,1)
  zFracTune = MOD(ring%ele(ring%n_ele_track)%mode3%c%phi,2.0_rp*pi)/2.0_rp / pi

  if (.not. allocated(bpm_save)) allocate(bpm_save(n_turns,n_bpms))
  if (n_TTs .gt. 0) then
     if (.not. allocated(tt_tunes)) allocate(tt_tunes(n_turns,n_TTs))
  endif
  allocate(phiz(n_bpms),Az(n_bpms))

! Tracking the particles
  if (present(ele_noise)) then
     call sim_tbt_data(ring, bpm, n_turns, .true., n_TTs, tt_params, id, n_damping, init_vec, bpm_save, tt_tunes, ele_noise)
  else
     call sim_tbt_data(ring, bpm, n_turns, .true., n_TTs, tt_params, id, n_damping, init_vec, bpm_save, tt_tunes)
  end if
      
! Calculate the amplitude and phase 
  call cal_butn_amp_phase(ring, bpm_save,tt_tunes, long)

  dc(:)%x = 0.
  dc(:)%y = 0.
  do ix=1, n_bpms
     ix_db = bpm(ix)%ix_db
     dc(ix_db)%signal(1) = sum(bpm_save(:,ix)%amp(1))/size(bpm_save(:,ix)%amp(1))
     dc(ix_db)%signal(2) = sum(bpm_save(:,ix)%amp(2))/size(bpm_save(:,ix)%amp(2))
     dc(ix_db)%signal(3) = sum(bpm_save(:,ix)%amp(3))/size(bpm_save(:,ix)%amp(3))
     dc(ix_db)%signal(4) = sum(bpm_save(:,ix)%amp(4))/size(bpm_save(:,ix)%amp(4))
  enddo


  call nonlin_xy_shake_components_cesrvidx (long, dc)

  freq%rev = 390.14
  freq%x%reflection = .false.
  freq%x%tune=abs(zFracTune)*390.14

! calculate the ac_eta 
  if(use_new_method) then
     ! new method to calculate the ac_eta
     call z_amp_phase_eta_calc(ring, bpm, long, freq, ac_eta, phiz, Az) 
  else
     ! fit the eta_x function to get the phase
     call ac_eta_calc(ring, bpm, n_bpms,long, freq, ac_eta) 
!     do i=1,n_bpms
!        phiz(i)=long(i)%z%phase*twopi/360
!        Az(i)=long(i)%z%amp
!     enddo
  endif

end subroutine sim_ac_eta_data

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

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!Subroutine sim_ac_and_cbar_data(ring, bpm, n_turns, n_TTs, tt_params, id, n_damping, init_vec, use_new_method, ac_eta, ele_noise)
! Simulate simultaneous ac_eta (dispersion) and phase/cbar measurement, using the same set of turns data. 
!
! Input:
!   ring  -- lat_struct
!   bpm(:) -- det_struct(:), holds information about a single BPM
!   n_turns -- integer, total number of turns to sample in the measurement
!   n_TTs    -- number of tune trackers. if zero, no damping/excitation used
!   tt_params -- type(tt_param_struct), for tune trackers.
!   n_damping -- integer. If damping = .true., how many turns to damp 
!                recording orbit data? 
!   init_vec -- real(rp)(1:6), in meters. How far to offset 
!                     particles before tracking. 
!   use_new_method --logical, .true. to use new method to calculate ac_eta
!                             .false. to use old method
!   ele_noise(:) -- ele_noise_struct(:), optional.  Noise to apply to elements
!                   when tracking.
!
! Output:
!   ac_eta   -- cesr_all_data_struct, calculated ac disperion
!--------------------------------------------------------------------------
!
Subroutine sim_ac_and_cbar_data(ring, bpm, n_turns, n_TTs, tt_params, id, n_damping, init_vec, use_new_method, ac_eta, phi, cbar, A, ele_noise)

    implicit none
    type(lat_struct) :: ring
    type(det_struct) :: bpm(:)
    integer, intent(in) :: n_turns
    integer, intent(in) :: n_TTs
    type(tt_param_struct) :: tt_params(:)
    integer, intent(in) :: id(:)
    integer, intent(in) :: n_damping
    real(rp), intent(in) :: init_vec(6)
    real(rp) :: harvest
    logical, intent(in) :: use_new_method
    type (cesr_all_data_struct) :: ac_eta 
    type(ele_noise_struct), optional :: ele_noise(:)

    type (det_data_struct), allocatable :: bpm_save(:,:)
    real(rp), allocatable :: tt_tunes(:,:), on_off_save(:)
    integer n_bpms,i
    real(rp) zFracTune
    type(coord_struct), allocatable :: orb(:), co(:)

! parameters for ac_eta calculation
    type (cesr_det_dc_position_struct) dc(0:n_det_maxx)
    type (cesr_freq_struct), target :: freq
   
    real(rp) ,allocatable :: phiz(:), Az(:)

! parameters for phase/cbar measurement
  real(rp), allocatable :: phi(:,:,:), cbar(:,:,:), A(:,:,:)
    
! 2011.05.25 - incorporate Suntao's button-based phase/amplitude notation:
  type (cesr_det_plane_struct) :: horz(0:n_det_maxx), vert(0:n_det_maxx)
  type (cesr_det_plane_struct) :: long(0:n_det_maxx)
  
  integer :: ix, jx, idet, ix_db, ix_lat
  real(rp) rel_amp, rel_phase, rms_x, rms_y
  logical err
  
 
  if (n_TTs .lt. 3) then
     write(*,*) "*****sim_ac_and_cbar_data subroutine requires three tune trackers (h,v,l)*****"
     write(*,'(a,i0,a)') "**********Only ", n_TTs, " provided. Consider yourself warned!****************"
  endif

  n_bpms = ubound(bpm,1)
  zFracTune = MOD(ring%ele(ring%n_ele_track)%mode3%c%phi,2.0_rp*pi)/2.0_rp / pi
  
  if (.not. allocated(A)) then
     allocate(A(n_bpms,2,2), cbar(n_bpms,2,2), phi(n_bpms,2,2))
  endif
  
  if (.not. allocated(bpm_save)) allocate(bpm_save(n_turns,n_bpms))
  
  if (n_TTs .gt. 0) then
     if (.not. allocated(tt_tunes)) allocate(tt_tunes(n_turns,n_TTs))
  endif
  
  phi(:,:,:) = 0 ! = (n_bpms,2,2) = (bpm%ix_db,(ax,ay),(bx,by))
  A(:,:,:) = 0   ! = (n_bpms,2,2) = (bpm%ix_db,(ax,ay),(bx,by))
  cbar(:,:,:) = 0! = (n_bpms,2,2) = (bpm%ix_db,(11,12),(21,22))
    
! Tracking the particles
  if (present(ele_noise)) then
     call sim_tbt_data(ring, bpm, n_turns, .true., n_TTs, tt_params, id, n_damping, init_vec, bpm_save, tt_tunes, ele_noise)
  else
     call sim_tbt_data(ring, bpm, n_turns, .true., n_TTs, tt_params, id, n_damping, init_vec, bpm_save, tt_tunes)
  end if


  ! calculate phase/amplitude for buttons:
  call cal_butn_amp_phase(ring, bpm_save,tt_tunes, horz, vert, long)

  if (.not. allocated(orb)) allocate(orb(0:ring%n_ele_track), co(0:ring%n_ele_track))
  call set_on_off(rfcavity$, ring, off_and_save$, saved_values = on_off_save)
  call closed_orbit_calc(ring, co, 4)
  call set_on_off(rfcavity$, ring, restore_state$, saved_values = on_off_save)

  dc(:)%x = 0.
  dc(:)%y = 0.
  rms_x   = 0.
  rms_y   = 0.

  do ix=1, n_bpms
     ix_db = bpm(ix)%ix_db
     ix_lat = bpm(ix)%ix_lat

     dc(ix_db)%signal(1) = sum(bpm_save(:,ix)%amp(1))/size(bpm_save(:,ix)%amp(1))
     dc(ix_db)%signal(2) = sum(bpm_save(:,ix)%amp(2))/size(bpm_save(:,ix)%amp(2))
     dc(ix_db)%signal(3) = sum(bpm_save(:,ix)%amp(3))/size(bpm_save(:,ix)%amp(3))
     dc(ix_db)%signal(4) = sum(bpm_save(:,ix)%amp(4))/size(bpm_save(:,ix)%amp(4))

     orb(ix_lat)%vec(1) =  sum(bpm_save(:,ix)%vec(1))/size(bpm_save(:,ix)%vec(1))
     orb(ix_lat)%vec(3) =  sum(bpm_save(:,ix)%vec(3))/size(bpm_save(:,ix)%vec(3))

     rms_x = rms_x + (orb(ix_lat)%vec(1) - co(ix_lat)%vec(1))**2/n_bpms
     rms_y = rms_y + (orb(ix_lat)%vec(3) - co(ix_lat)%vec(3))**2/n_bpms

  enddo

  rms_x = sqrt(rms_x)
  rms_y = sqrt(rms_y)

  write(*,*) " " 
  write(*,*) "========================================================================"
  write(*,*) "RMS (measured - Bmad reported):" 

  ! for ring_ma2, technically we use different data for measuring the orbit. this only gives a sort of 
  ! sanity check, to make sure we aren't way off.
  write(*,'(a,es10.3,a,es10.3)') " RMS delta(orbit.x)   = ", rms_x, "    RMS delta(orbit.y)   = ", rms_y


!  do ix = 1, size(bpm_save(1,:))
!     ix_db = bpm(ix)%ix_db
     !rms_x = rms_x + (bpm_save(,)**2
!  enddo

  
!------------------------------
! AC_ETA CALCULATIONS: 

  call nonlin_xy_shake_components_cesrvidx (long, dc)

  freq%rev = 390.14
  freq%x%reflection = .false.
  freq%x%tune=abs(zFracTune)*390.14

  if(use_new_method) then
     ! new method to calculate the ac_eta
     call z_amp_phase_eta_calc(ring, bpm, long, freq, ac_eta, phiz, Az) 
  else
     ! fit the eta_x function to get the phase
     call ac_eta_calc(ring, bpm, n_bpms,long, freq, ac_eta) 
     !do i=1,n_bpms
     !   phiz(i)=long(i)%z%phase*twopi/360
     !   Az(i)=long(i)%z%amp
     !enddo
  endif

! END AC_ETA CALCULATIONS
!------------------------------
! PHASE/COUPLING CALCULATIONS:
  
  call calc_phase_cbar(ring, bpm_save, horz, vert, phi, cbar, A)

  write(*,*) "========================================================================"
  write(*,*) " "

end subroutine sim_ac_and_cbar_data


!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Subroutine calc_phase_cbar(ring, bpm_save, horz, vert, phi, cbar, A)
!
! Subroutine to compute phase / cbar data from provided turn-by-turn data.
! To be used internally by sim_cbar_data and sim_ac_and_cbar_data. 
!
! Input:
!   ring     - lat_struct
!   bpm_save - det_struct, dimensions (1:n_turns, 1:n_bpms); holds TBT data
!   horz, vert - cesr_det_plane_struct, dimeensions (0:n_det_maxx). holds 
!                pre-processed data from cal_butn_amp_phase
! 
! Returns:
!   phi      - real, dimensions (1:n_bpms, 1:2, 1:2); holds phase data
!   cbar     - real, dimensions (1:n_bpms, 1:2, 1:2); holds cbar data
!   A        - real, dimensions (1:n_bpms, 1:2, 1:2); holds amplitude data 
!              from phase calculations
!--------------------------------------------------------------------------

subroutine calc_phase_cbar(ring, bpm_save, horz, vert, phi, cbar, A)


  implicit none
  type(lat_struct) :: ring
  type(det_data_struct) :: bpm_save(:,:)

  real(rp), allocatable :: phi(:,:,:), cbar(:,:,:), A(:,:,:), cbar_bmad(:,:,:)
  
  ! 2011.05.25 - incorporate Suntao's button-based phase/amplitude notation:
  type (cesr_det_plane_struct) :: horz(0:n_det_maxx), vert(0:n_det_maxx)
  type (cesr_det_dc_position_struct) :: dc(0:n_det_maxx)
  
  integer :: i, ix, jx, n_bpms, idet, ix_db, ix_lat
  real(rp) rel_amp, rel_phase, dphi_ax_avg, dphi_by_avg
  real(rp), allocatable :: delta_phi_ax(:), delta_phi_by(:)
  real(rp) :: rms_phia, rms_phib, rms_ac_x, rms_ac_y, rms_cbar(2,2)

  n_bpms = size(bpm_save(1,:))
 
  if (.not. allocated(phi)) then
     allocate(A(n_bpms,2,2), cbar(n_bpms,2,2), phi(n_bpms,2,2))
  endif

  if (.not. allocated(cbar_bmad)) allocate(cbar_bmad(n_bpms,2,2))

  dc(:)%x = 0.
  dc(:)%y = 0.
  do ix=1, n_bpms
     ix_db = bpm_save(1,ix)%ix_db
     dc(ix_db)%signal(1) = sum(bpm_save(:,ix)%amp(1))/size(bpm_save(:,ix)%amp(1))
     dc(ix_db)%signal(2) = sum(bpm_save(:,ix)%amp(2))/size(bpm_save(:,ix)%amp(2))
     dc(ix_db)%signal(3) = sum(bpm_save(:,ix)%amp(3))/size(bpm_save(:,ix)%amp(3))
     dc(ix_db)%signal(4) = sum(bpm_save(:,ix)%amp(4))/size(bpm_save(:,ix)%amp(4))
  enddo

  ! convert b(4) phase/amp data to h/v * a/b data
  call nonlin_xy_shake_components_cesrvidx (horz, dc)
  call nonlin_xy_shake_components_cesrvidx (vert, dc)
  
  ! from CESRv / Suntao's code:
  do idet=1, n_bpms

     ix_db  = bpm_save(1,idet)%ix_db
     ix_lat = bpm_save(1,idet)%ix_lat

     if (horz(ix_db)%x%amp /= 0) then
        rel_amp = horz(ix_db)%y%amp / horz(ix_db)%x%amp
        rel_phase = (horz(ix_db)%y%phase - horz(ix_db)%x%phase)*pi/180.
        horz(ix_db)%cbar22 = -rel_amp * cos(rel_phase)
        horz(ix_db)%cbar12 = rel_amp * sin(rel_phase)
     else
        horz(ix_db)%cbar22 = 0
        horz(ix_db)%cbar12 = 0
     endif
     
     if (vert(ix_db)%y%amp /= 0) then
        rel_amp = vert(ix_db)%x%amp / vert(ix_db)%y%amp
        rel_phase = (vert(ix_db)%x%phase - vert(ix_db)%y%phase)*pi/180.
        vert(ix_db)%cbar11 = rel_amp * cos(rel_phase)
        vert(ix_db)%cbar12 = rel_amp * sin(rel_phase)
     else
        vert(ix_db)%cbar11 = 0
        vert(ix_db)%cbar12 = 0
     endif
     
     vert(ix_db)%cbar12 = vert(ix_db)%cbar12 / sqrt(ring%ele(ix_lat)%a%beta / ring%ele(ix_lat)%b%beta)
     horz(ix_db)%cbar22 = horz(ix_db)%cbar22 * sqrt(ring%ele(ix_lat)%a%beta / ring%ele(ix_lat)%b%beta)
     vert(ix_db)%cbar11 = vert(ix_db)%cbar11 / sqrt(ring%ele(ix_lat)%a%beta / ring%ele(ix_lat)%b%beta)
     
     ! straight from cesrv...
     call button_phase (horz(ix_db), 0.0_rp, x_plane$)
     call button_phase (vert(ix_db), 0.0_rp, y_plane$)
     
  enddo

  ! Load results into bsim_cesr structures:
  A(:,1,1)    = horz(bpm_save(1,:)%ix_db)%x%amp
  A(:,1,2)    = vert(bpm_save(1,:)%ix_db)%x%amp
  A(:,2,1)    = horz(bpm_save(1,:)%ix_db)%y%amp
  A(:,2,2)    = vert(bpm_save(1,:)%ix_db)%y%amp

  do idet=1,n_bpms
     ix_db = bpm_save(1,idet)%ix_db
     cbar(idet,1,2) = vert(ix_db)%cbar12 
     cbar(idet,2,2) = horz(ix_db)%cbar22 
     cbar(idet,1,1) = vert(ix_db)%cbar11 
  enddo
  
  
  ! correct the phase advance for when atan2 'wraps around' from pi --> -pi

  do ix = 1, n_bpms
     ix_db = bpm_save(1,ix)%ix_db
     ix_lat = bpm_save(1,ix)%ix_lat
     horz(ix_db)%phase_design = ring%ele(ix_lat)%a%phi * 180./pi
     vert(ix_db)%phase_design = ring%ele(ix_lat)%b%phi * 180./pi
     horz(ix_db)%phase_meas   = horz(ix_db)%x%phase 
     vert(ix_db)%phase_meas   = vert(ix_db)%y%phase 

  enddo

  call phase_shift_near_design(horz)
  call phase_shift_near_design(vert)


  phi(:,1,1) = horz(bpm_save(1,:)%ix_db)%phase_meas * pi/180.
  phi(:,1,2) = horz(bpm_save(1,:)%ix_db)%y%phase    * pi/180.
  phi(:,2,1) = vert(bpm_save(1,:)%ix_db)%x%phase    * pi/180.
  phi(:,2,2) = vert(bpm_save(1,:)%ix_db)%phase_meas * pi/180.


  if (.not. allocated(delta_phi_ax)) allocate(delta_phi_ax(n_bpms),delta_phi_by(n_bpms))
  delta_phi_ax(:) = 0
  delta_phi_by(:) = 0
      
  delta_phi_ax(:) = ring%ele(bpm_save(1,:)%ix_lat)%a%phi - phi(:,1,1)
  delta_phi_by(:) = ring%ele(bpm_save(1,:)%ix_lat)%b%phi - phi(:,2,2)
        
  dphi_ax_avg = (sum(delta_phi_ax(:))) / (size(delta_phi_ax))
  dphi_by_avg = (sum(delta_phi_by(:))) / (size(delta_phi_by))

  phi(:,1,1)  = phi(:,1,1) + dphi_ax_avg
  phi(:,2,2)  = phi(:,2,2) + dphi_by_avg

  ! useful feedback for the user


  do ix = 1, n_bpms
     call c_to_cbar(ring%ele(bpm_save(1,ix)%ix_lat),cbar_bmad(ix,:,:))
  enddo

  rms_cbar(:,:) = 0
  rms_phia = 0
  rms_phib = 0

  do ix=1,n_bpms
     ix_db = bpm_save(1,ix)%ix_db
     if ((horz(ix_db)%ok .neqv. .true.) .or. (vert(ix_db)%ok .neqv. .true.)) cycle
     rms_phia = rms_phia + (phi(ix,1,1) - ring%ele(bpm_save(1,ix)%ix_lat)%a%phi)**2
     rms_phib = rms_phib + (phi(ix,2,2) - ring%ele(bpm_save(1,ix)%ix_lat)%b%phi)**2
     rms_cbar(1,1) = rms_cbar(1,1) + (cbar_bmad(ix,1,1)-cbar(ix,1,1))**2
     rms_cbar(1,2) = rms_cbar(1,2) + (cbar_bmad(ix,1,2)-cbar(ix,1,2))**2
     rms_cbar(2,2) = rms_cbar(2,2) + (cbar_bmad(ix,2,2)-cbar(ix,2,2))**2
  enddo

  n_bpms = min(count(horz(:)%ok),count(vert(:)%ok))

  rms_phia = sqrt(rms_phia/(n_bpms))
  rms_phib = sqrt(rms_phib/(n_bpms))
  rms_cbar(1,1) = sqrt(rms_cbar(1,1)/(n_bpms))
  rms_cbar(1,2) = sqrt(rms_cbar(1,2)/(n_bpms))
  rms_cbar(2,2) = sqrt(rms_cbar(2,2)/(n_bpms))

  write(*,'(a,es10.3,a,es10.3)')          " RMS delta(phi.a)     = ",      rms_phia, "    RMS delta(phi.b)     = ", rms_phib
  write(*,'(a,es10.3,a,es10.3,a,es10.3)') " RMS delta(cbar12)    = ", rms_cbar(1,2)
  write(*,'(a,es10.3,a,es10.3,a,es10.3)') " RMS delta(cbar22)    = ", rms_cbar(2,2), "    RMS delta(cbar11)    = ", rms_cbar(1,1)

end subroutine calc_phase_cbar


! ================================================

subroutine button_phase(det, amp_minn, hv)

implicit none

type (cesr_det_plane_struct) det
real(rp), save :: amp_min_default = 1000

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


! flip phases

phase = det%but(:)%phase

! Depending on horizontal or vertical shaking, two buttons will be 
! 180deg out-of-phase. Correct for this here.
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

! =======================================================================

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, n_iter, n_iter_tot

! 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

n_iter = 0

n_iter_tot = mod(maxval(design(:)), 360._rp)

do
  n_iter = n_iter + 1
  offset = sum(meas - design, mask = plane(:)%ok) / n
  if (all(meas - design - offset < 180._rp)) exit
  meas = meas + 360._rp *  nint((design + offset - meas) / 360._rp) - offset
  if (n_iter .gt. n_iter_tot * 2) then ! not converging-- stop here
     write(*,*) "phase_shift_near_design: not converging!"
     if (global_com%exit_on_error) then
        stop
     else
        return
     endif
  endif
enddo

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

! See if shifting the offset improves anything

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



end subroutine phase_shift_near_design


!++The following subroutines are for sim_ac_eta_data++++++++++++++++++++++++++++++++++++++++++++++

!====================================================================================================
subroutine z_amp_phase_eta_calc(ring, bpm, long, freq, ac_eta, phiz, Az)

  type (lat_struct) :: ring
  type (det_struct) :: bpm(:)
  type (cesr_freq_struct), target :: freq
  
  type (cesr_det_plane_struct) long(0:n_det_maxx)
  type (cesr_all_data_struct) ac_eta
  real(rp), allocatable :: phiz(:), Az(:)
  
  integer n_bpms
  integer i, ix, ix_db
  real(rp) eta_x,etap_x, eta_y, etap_y,eta_x_meas,etap_x_meas,eta_y_meas,etap_y_meas, beta_c,alpha_c
  real(rp) temp, temp_phi
  real(rp) rms_ac_x, rms_ac_y

  n_bpms=size(bpm)

  if (.not. allocated(phiz)) allocate(phiz(1:n_bpms),Az(1:n_bpms))

  phiz(:)=0
  Az(:)=0

  do i = 1, n_bpms

     ix = bpm(i)%ix_lat
     ix_db = bpm(i)%ix_db

     if (ix < 1) cycle
     if (.not. long(ix_db)%ok ) cycle
     
     eta_x = ring%ele(ix)%mode3%x%eta
     etap_x = ring%ele(ix)%mode3%x%etap
     eta_y = ring%ele(ix)%mode3%y%eta
     etap_y = ring%ele(ix)%mode3%y%etap

     beta_c=ring%ele(ix)%mode3%c%beta
     alpha_c=ring%ele(ix)%mode3%c%alpha

     temp=beta_c*etap_x/eta_x-alpha_c
     temp_phi=atan(1/temp) !radians between -pi and pi
     phiz(i)=long(ix_db)%x%phase*twopi/360-temp_phi
     Az(i)=beta_c*long(ix_db)%x%amp*sin(temp_phi)/eta_x

     eta_x_meas=beta_c*long(ix_db)%x%amp*sin(temp_phi)/Az(i)
     etap_x_meas=(cos(temp_phi)+alpha_c*sin(temp_phi))*long(ix_db)%x%amp/Az(i)
       
     temp_phi=long(ix_db)%y%phase*twopi/360-phiz(i)
     eta_y_meas=beta_c*long(ix_db)%y%amp*sin(temp_phi)/Az(i)
     etap_y_meas=(cos(temp_phi)+alpha_c*sin(temp_phi))*long(ix_db)%y%amp/Az(i)
  
     ac_eta%ac_etap_x(i)%value = etap_x_meas
     ac_eta%ac_eta_x(i)%value = eta_x_meas
     ac_eta%ac_etap_y(i)%value = etap_y_meas
     ac_eta%ac_eta_y(i)%value = eta_y_meas

! correct the phase advance for when atan 'wraps around' from pi --> -pi
     do ix=2,n_bpms     
        do while (phiz(ix) .lt. phiz(ix-1))
           phiz(ix) = phiz(ix) + twopi
        enddo
     enddo

  enddo


  ! useful feedback for the user: 
  rms_ac_x = 0
  rms_ac_y = 0

  do ix=1,n_bpms
     rms_ac_x = rms_ac_x + (ac_eta%ac_eta_x(ix)%value - ring%ele(bpm(ix)%ix_lat)%mode3%x%eta)**2
     rms_ac_y = rms_ac_y + (ac_eta%ac_eta_y(ix)%value - ring%ele(bpm(ix)%ix_lat)%mode3%y%eta)**2
  enddo

  rms_ac_x = sqrt(rms_ac_x / n_bpms)
  rms_ac_y = sqrt(rms_ac_y / n_bpms)

  write(*,'(a,es10.3,a,es10.3)') " RMS delta(ac_eta.x)  = ", rms_ac_x, "    RMS delta(ac_eta.y)  = ", rms_ac_y



end subroutine z_amp_phase_eta_calc

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!  Subroutine  ac_eta_calc (ring, bpm, n_bpms,long, freq, ac_eta)
!
! Fit amplitude and overall phase of longitudinal-horizontal data to model 
! longitudinal-horizontal. Then use fitted values along with measured 
! longitudinal-vertical phase and amplitude data to get longitudinal-vertical coupling
!
! Input:
!   ring  -- lat_struct
!   long  -- cesr_det_plane_struct: Measured horizontal and vertical amp and phase
!   freq  -- cesr_freq_struct: Shaking frequency.
!   bpm   -- det_struct:  cesr bpm structure (sim_bpm_mod)
!   n_bpms -- number of bpms
!
! Output:
!   ac_eta   -- cesr_all_data_struct - ac eta, etap, x, y
!                                      Fitted Amplitude and phase
!                                      chis2 of fit
!                                      data point that contributes most to figure of merit
!-

subroutine ac_eta_calc (ring, bpm, n_bpms, long, freq, ac_eta)

use cesr_basic_mod
use sim_bpm_mod
use mode3_mod

implicit none

type (lat_struct) ring
type (cesr_freq_struct), target :: freq
type (cesr_det_plane_struct) :: long(0:n_det_maxx)
type (cesr_all_data_struct) ac_eta
type (det_struct) :: bpm(:)
integer n_bpms

real(rp) del_phi
real(rp) pp(1:2), gradient(2), chis_0, chis, grad_mag, b
real(rp) pp_min(1:2), chis_min
logical error
real(rp) rms_ac_x, rms_ac_y

integer i, j, ix

! mode3 calc using the measured tune

!call set_on_off (rfcavity$, ring, on$) !assume this is done during initialization in top level
call set_z_tune (ring%branch(0), -freq%x%tune * twopi / 390.14)

call twiss3_at_start (ring, error)
if (error) call err_exit
call twiss3_propagate_all (ring)

call calc_all_c(ring, bpm, n_bpms, long, ac_eta)

chis_min = 999.
do i = 1, 1000

  call get_gradient (ring, bpm, n_bpms, long, ac_eta, gradient)

  grad_mag = (gradient(1)**2)
  chis_0 = ac_eta%param%chisq
  b = -chis_0/grad_mag/i

  pp(1) = ac_eta%param%ac_z_phase_fit
  pp(1) =  gradient(1)*b + pp(1)

  ac_eta%param%ac_z_phase_fit = pp(1)

  call calc_all_c(ring, bpm, n_bpms, long, ac_eta)

  chis = ac_eta%param%chisq
  pp(1) = ac_eta%param%ac_z_phase_fit

  if (chis < chis_min) then
    chis_min = chis
    pp_min(1) = pp(1)
  endif

  if (abs(chis_0-chis) < 0.0001)exit
  chis_0 = chis
end do

pp(1) = pp_min(1)

ac_eta%param%ac_z_phase_fit = pp(1)

call calc_all_c(ring, bpm, n_bpms, long,  ac_eta)

  ! useful feedback for the user: 
  rms_ac_x = 0
  rms_ac_y = 0

  do ix=1,n_bpms
     rms_ac_x = rms_ac_x + (ac_eta%ac_eta_x(ix)%value - ring%ele(bpm(ix)%ix_lat)%mode3%x%eta)**2
     rms_ac_y = rms_ac_y + (ac_eta%ac_eta_y(ix)%value - ring%ele(bpm(ix)%ix_lat)%mode3%y%eta)**2
  enddo

  rms_ac_x = sqrt(rms_ac_x / n_bpms)
  rms_ac_y = sqrt(rms_ac_y / n_bpms)

  write(*,*) " "
  write(*,'(a,es10.3,a,es10.3)') " RMS delta(ac_eta.x)  = ", rms_ac_x, "    RMS delta(ac_eta.y)  = ", rms_ac_y

end subroutine ac_eta_calc

!------------------------------------------------------
!+
! Convert measured amplitude and phase to V-matrix elements
! Input:
!   long -- Cesr_det_plane_struct: x, y, z phase and amplitude.
!   ele  -- Ele_struct: Element containing normal mode beta values.
!
! Output:
!   v(6,6) -- real(rp) V-matrix with coupling x-y, x-z, y-z.
!-

subroutine amp_phase_to_cbar(long, ele, V)

use cesr_basic_mod
use sim_bpm_mod

implicit none

type (cesr_det_plane_struct) :: long
type g_struct
  real(rp) mat(1:2,1:2)
endtype
type (g_struct) C(3)
type (ele_struct) ele

real(rp) rbeta(3)
real(rp) cbar(2,2)
real(rp) rel_amp, rel_phase
real(rp) DegToRad
real(rp) V(6,6)

DegToRad = twopi/360.
rbeta(1) = sqrt(ele%mode3%a%beta)
rbeta(2) = sqrt(ele%mode3%b%beta)
rbeta(3) = sqrt(ele%mode3%c%beta)

cbar(1:2,1:2)=0.

rel_amp = (long%x%amp/long%z%amp)
rel_phase = (long%x%phase - long%z%phase) *DegToRad
Cbar(1,1) = rel_amp*cos(rel_phase)*rbeta(3)/rbeta(1)
Cbar(1,2) = rel_amp*sin(rel_phase)*rbeta(3)/rbeta(1)
call cbar_to_c (Cbar, ele%mode3%a, ele%mode3%c, C(2)%mat)

rel_amp = (long%y%amp/long%z%amp)
rel_phase = (long%y%phase - long%z%phase)*DegToRad
Cbar(1,1) = rel_amp*cos(rel_phase)*rbeta(3)/rbeta(2)
Cbar(1,2) = rel_amp*sin(rel_phase)*rbeta(3)/rbeta(2)
call cbar_to_c(Cbar, ele%mode3%b, ele%mode3%c, C(3)%mat)

V(1:6,1:6)=0.
V(1:2,3:4) = C(1)%mat(1:2,1:2)
V(1:2,5:6) = C(2)%mat(1:2,1:2)
V(3:4,5:6) = C(3)%mat(1:2,1:2)

end subroutine amp_phase_to_cbar

!------------------------------------------------------------------
subroutine calc_all_c(lat, bpm, n_bpms, long, ac_eta)

use cesr_basic_mod
use mode3_mod
use sim_bpm_mod
implicit none

type (lat_struct) lat
type (cesr_det_plane_struct) :: long(0:n_det_maxx)
type (cesr_all_data_struct) ac_eta
type (twiss_struct) twiss
type (det_struct) :: bpm(:)
integer n_bpms

real(rp) A, dphi, chis
real(rp) c2_12_calc
real(rp) V(6,6)
real(rp) contrib, contrib_0
real(rp) sum1, sum2
real(rp) B
real(rp) eta_calc(0:120), eta_meas(0:120)

integer i, ix, ix_db
integer j, n
integer isave

logical error

A = ac_eta%param%ac_z_amp_fit
A=0.5
!A= 1.
dphi = ac_eta%param%ac_z_phase_fit

contrib=0.
chis=0.
n=0

! Loop over all detectors...

sum1 = 0
sum2 = 0
chis = 0

do i = 1, n_bpms

  ix = bpm(i)%ix_lat
  ix_db = bpm(i)%ix_db

  if (ix < 1) cycle
  if (.not. long(i)%ok ) cycle

  twiss = lat%ele(ix)%mode3%c
 
  long(ix_db)%z%amp = A * sqrt(twiss%beta)
  long(ix_db)%z%phase = (twiss%phi - dphi) *360/twopi

  call amp_phase_to_cbar(long(ix_db), lat%ele(ix), V)

  ac_eta%ac_etap_x(i)%value = V(1,5)
  ac_eta%ac_eta_x(i)%value = V(1,6)
  ac_eta%ac_etap_y(i)%value = V(3,5)
  ac_eta%ac_eta_y(i)%value = V(3,6)

  c2_12_calc = lat%ele(ix)%mode3%v(1,6)
  eta_calc(i) = c2_12_calc
  eta_meas(i) = ac_eta%ac_eta_x(i)%value

  if (long(i)%ok  ) then
    n = n+1
    sum1 = sum1 + ac_eta%ac_eta_x(i)%value * c2_12_calc
    sum2 = sum2 + ac_eta%ac_eta_x(i)%value**2
    chis = chis + (ac_eta%ac_eta_x(i)%value - c2_12_calc)**2
  endif
end do

B = sum1/sum2
chis=0.
n=0
contrib_0 = 0

do i=1, n_bpms
  ix_db = bpm(i)%ix_db

  ac_eta%ac_etap_x(i)%value = ac_eta%ac_etap_x(i)%value*B
  ac_eta%ac_eta_x(i)%value = ac_eta%ac_eta_x(i)%value*B
  ac_eta%ac_etap_y(i)%value = ac_eta%ac_etap_y(i)%value*B
  ac_eta%ac_eta_y(i)%value = ac_eta%ac_eta_y(i)%value*B 

  if (long(ix_db)%ok ) then
    n=n+1
    contrib = (B*eta_meas(i) - eta_calc(i))**2

    if (contrib > contrib_0) then
      isave = i
      contrib_0 = contrib
    endif

    chis = chis + (B*eta_meas(i) - eta_calc(i))**2
  endif
end do

A= 1./B

chis = chis/float(n)
!print '(a13,i,a7,i,2e12.4)',' system_id = ',ac_eta%param%ix_data_set,' n,B,A ',n,B,A
ac_eta%param%ac_z_amp_fit = A
ac_eta%param%ac_z_phase_fit = dphi
!ac_eta%param%ix_data_set = isave
ac_eta%param%chisq = chis

end subroutine calc_all_c

!--------------------------------------------------------
subroutine get_gradient(lat, bpm, n_bpms, long, ac_eta, gradient)

use cesr_basic_mod
use mode3_mod
use sim_bpm_mod

implicit none

type (cesr_det_plane_struct) :: long(0:n_det_maxx)
type (cesr_all_data_struct) ac_eta
type (lat_struct) lat
 type (det_struct):: bpm(:)
integer n_bpms

real(rp) del_phi, A, dphi
real(rp)  pp(1:2), gradient(2), chis_0, chis
real(rp) deltap(2)

integer i
integer j


call calc_all_c(lat, bpm, n_bpms, long, ac_eta)

pp(1) = ac_eta%param%ac_z_phase_fit
chis_0 = ac_eta%param%chisq
 

i = 1
deltap(i) = 0.001
pp(i) = pp(i)+deltap(i)

ac_eta%param%ac_z_phase_fit = pp(1)

call calc_all_c(lat, bpm, n_bpms, long, ac_eta)

chis = ac_eta%param%chisq

pp(i) = pp(i)-deltap(i)
gradient(i) = (chis-chis_0)/deltap(i)

end subroutine get_gradient

!======================================================================================
subroutine cal_butn_amp_phase(ring,bpm_save,tt_tunes, horz, vert, long)

  type(lat_struct) :: ring
  type(det_data_struct) :: bpm_save(:,:)
  real(rp) :: tt_tunes(:,:)

  type(cesr_det_plane_struct) :: horz(0:n_det_maxx)
  type(cesr_det_plane_struct), optional :: vert(0:n_det_maxx), long(0:n_det_maxx)

  integer n_turns, n_bpms
  real(rp), allocatable :: cax(:,:),sax(:,:),  cby(:,:), sby(:,:), ccz(:,:), scz(:,:)
  integer ix,jx, ix_db
  logical err

  n_turns = size(bpm_save(:,1))
  n_bpms  = size(bpm_save(1,:))

  if (.not. allocated(cax)) then
     allocate(cax(4,n_bpms),sax(4,n_bpms),cby(4,n_bpms),sby(4,n_bpms),ccz(4,n_bpms),scz(4,n_bpms))
  endif
  
  ! Notation: c/s   = cosine / sine component
  !           a/b/c = mode
  !           x/y/z = physical axes
  
  cax(:,:) = 0
  sax(:,:) = 0
  cby(:,:) = 0
  sby(:,:) = 0
  ccz(:,:) = 0
  scz(:,:) = 0
  
  turns_loop1:do ix=1,n_turns
     buttns_loop1:do jx=1,4
        ! in the end, we take ratios of these terms, so factors of n_turns will drop out.
        cax(jx,:) = cax(jx,:) + bpm_save(ix,:)%amp(jx)*cos(sum(tt_tunes(1:ix,1)))  ! fractional tune is stored, in radians 
        sax(jx,:) = sax(jx,:) + bpm_save(ix,:)%amp(jx)*sin(sum(tt_tunes(1:ix,1)))
        if (size(tt_tunes(1,:)) .eq. 1) cycle buttns_loop1 ! go on to next iteration unless n_TTs > 1
        cby(jx,:) = cby(jx,:) + bpm_save(ix,:)%amp(jx)*cos(sum(tt_tunes(1:ix,2)))
        sby(jx,:) = sby(jx,:) + bpm_save(ix,:)%amp(jx)*sin(sum(tt_tunes(1:ix,2)))
        if (size(tt_tunes(1,:)) .eq. 2) cycle buttns_loop1 ! go on to next iteration unless n_TTs > 2
        ccz(jx,:) = ccz(jx,:) + bpm_save(ix,:)%amp(jx)*cos(sum(tt_tunes(1:ix,3)))
        scz(jx,:) = scz(jx,:) + bpm_save(ix,:)%amp(jx)*sin(sum(tt_tunes(1:ix,3)))
     enddo buttns_loop1

  enddo turns_loop1
    
  sax(:,:)=-sax(:,:)  ! this is correct! need to flip sign
  sby(:,:)=-sby(:,:)
  scz(:,:)=-scz(:,:)

  ! set default to "not ok"; flag good data points individually
  horz(:)%ok = .false.
  if (size(tt_tunes(1,:)) .gt. 1) then
     vert(:)%ok = .false.
     vert(:)%system_id = 0 ! default to 'old' system
  endif
  if (size(tt_tunes(1,:)) .gt. 2) then
     long(:)%ok = .false.
     long(:)%system_id = 0
  endif
 
  ! cax(:) ~ 0.5 * n_turns * cos(phi(:)); same idea for other terms
  !
  turns_loop2:do ix=1,n_bpms
     ix_db = bpm_save(1,ix)%ix_db
     horz(ix_db)%ok = .true.
     if (size(tt_tunes(1,:)) .gt. 1) then
        vert(ix_db)%ok = .true.
        vert(ix_db)%system_id = 1
     endif
     if (size(tt_tunes(1,:)) .gt. 2) then
        long(ix_db)%ok = .true.
        long(ix_db)%system_id = 1
     endif
     buttns_loop2:do jx=1,4 
        horz(ix_db)%but(jx)%phase = atan2(sax(jx,ix),cax(jx,ix))*360._rp/twopi
        horz(ix_db)%but(jx)%amp =  (2./n_turns)*sqrt((sax(jx,ix))**2 + (cax(jx,ix))**2) 
        if (size(tt_tunes(1,:)) .eq.1 ) cycle buttns_loop2 ! go on to next iteration unless n_TTs > 1
        vert(ix_db)%but(jx)%phase = atan2(sby(jx,ix),cby(jx,ix))*360._rp/twopi
        vert(ix_db)%but(jx)%amp =  (2./n_turns)*sqrt((sby(jx,ix))**2 + (cby(jx,ix))**2)
        if (size(tt_tunes(1,:)) .eq. 2) cycle buttns_loop2 ! go on to next iteration unless n_TTs > 2
        long(ix_db)%but(jx)%phase = atan2(scz(jx,ix),ccz(jx,ix))*360._rp/twopi
        long(ix_db)%but(jx)%amp =  (2./n_turns)*sqrt((scz(jx,ix))**2 + (ccz(jx,ix))**2)
     enddo buttns_loop2
  enddo turns_loop2


end subroutine cal_butn_amp_phase



end module sim_ac_and_cbar_mod