subroutine wave_phase (phase_z, tune_z, wave)

use cesrv_struct
use cesrv_interface

implicit none

type (d1_data_struct), target :: phase_z, tune_z
type (wave_struct), target :: wave

type (p2_plot_struct), target :: tune_plot
type (p1_plot_struct), pointer :: plot1

real(rp) phi, cos_phi(250), sin_phi(250)
real(rp) coef_a(4), coef_b(4), const_(250), dummy_(250)
real(rp) phi_(250), coef_ba(4), amp_ba, amp_a, amp_b, phi_0
real(rp) rms_a(5), rms_b(5), rms_ba(5), var_tune, c
real(rp) phi_ring, x

integer i_a1, i_a2, i_b1, i_b2, n_a, n_b, m_min, m_max, nc
integer i, ix, j, n_pts_min / 4 /, n_data_use

character fmt*80

! Load wave arrays
! Remember to take out the conversion factor put in by plotting_data_calc

plot1 => wave%p2%plot1
call plotting_data_calc (plot1, phase_z, tune_z)
plot1%y = plot1%y / plot1%conversion_factor

! calculate the extra phase to add to elements with index over 100

tune_plot%plot1%p2 => tune_plot
tune_plot%plot1%d1 => tune_z
tune_plot%d2 => tune_z%d2
tune_plot%plot1%conversion_factor = 1
tune_plot%plot1%y_limit = -1
tune_plot%plot_data = wave%p2%plot_data
tune_plot%base      = wave%p2%base
tune_plot%plot1%normalize = wave%p2%plot1%normalize
call plotting_data_calc (tune_plot%plot1, phase_z, tune_z)
var_tune = tune_plot%plot1%y(1)

!

phi_ring = 2 * tune_z%d(1)%model

n_data_use = plot1%n_use
if (n_data_use == 0) then
  print *, 'ERROR: NO DATA FOR THE WAVE ANALYSIS.'
  return
endif

plot1%n_use = 0

do i = 1, n_data_use

  ix = wave%p2%ix(i)
  x = wave%p2%x(i)

  phi = 2 * phase_z%d(ix)%model

  phi_(i) = phi
  cos_phi(i) = cos(phi)
  sin_phi(i) = sin(phi)
  const_(i) = 1.0

  plot1%y(i+n_data_use) = plot1%y(i) + var_tune
  wave%p2%ix(i+n_data_use) = wave%p2%ix(i) + 100
  wave%p2%x(i+n_data_use)  = wave%p2%x(i) + 100
  phi_(i+n_data_use) = phi + phi_ring
  cos_phi(i+n_data_use) = cos(phi + phi_ring)
  sin_phi(i+n_data_use) = sin(phi + phi_ring)
  const_(i+n_data_use) = 1.0

  if (x > 50 .and. plot1%n_use == 0) plot1%n_use = i-1+n_data_use

enddo

! For the init: Compute the least squares model in a sliding window...

!-----------------------------------------
! Fit to regions...
! First: Find how many points are in the regions and where they start.
! Second: Find model coeficients COEF_A, and COEF_B.
! Third: Compute residuals WAVE%P2%a, and VAR_B.

call wave_region_setup (wave, i_a1, i_b1, n_a, n_b)

if (n_a < n_pts_min) then
  print *, 'ERROR: REGION A HAS LESS THAN', n_pts_min, ' POINTS'
  return
elseif (n_b < n_pts_min .and. (wave%ix_b1 /= 0 .or. wave%ix_b2 /= 0)) then
  print *, 'ERROR: REGION B HAS LESS THAN', n_pts_min, ' POINTS'
  return
endif

call wave_anal (plot1%y(i_a1:), sin_phi(i_a1:), cos_phi(i_a1:), &
                     const_(i_a1:), dummy_(i_a1:), 3, n_a, coef_a, rms_a)
call wave_anal (plot1%y(i_b1:), sin_phi(i_b1:), cos_phi(i_b1:), &
                     const_(i_b1:), dummy_(i_b1:), 3, n_b, coef_b, rms_b)

wave%p2%plot2%n_use = plot1%n_use
wave%p2%plot3%n_use = plot1%n_use

do i = 1, plot1%n_use
  wave%p2%plot2%y(i) = (plot1%y(i) - coef_a(1) * sin_phi(i) - &
            coef_a(2) * cos_phi(i) - coef_a(3) * const_(i)) * 180 / pi
  wave%p2%plot3%y(i) = (plot1%y(i) - coef_b(1) * sin_phi(i) - &
            coef_b(2) * cos_phi(i) - coef_b(3) * const_(i)) * 180 / pi
  if (n_b < n_pts_min) wave%p2%plot3%y(i) = 0
enddo

! compute the possible places between the regions where there could be a kicker
! This is when A_wave - B_wave = 0

if (n_b < n_pts_min) then
  wave%n_cross = 0
  return
endif

do i = 1, 3
  coef_ba(i) = coef_b(i) - coef_a(i)
  rms_ba(i) = sqrt(rms_a(i)**2 + rms_b(i)**2)
enddo

amp_a = sqrt(coef_a(1)**2 + coef_a(2)**2)
amp_b = sqrt(coef_b(1)**2 + coef_b(2)**2)
amp_ba = sqrt(coef_ba(1)**2 + coef_ba(2)**2)

wave%kick(1) = -(sign(amp_ba, coef_ba(3)) + coef_ba(3))

! k is k of quad so positive k means horizontal focusing.

if (wave%plane == x_plane$) wave%kick(1) = -wave%kick(1)

if (coef_ba(1) == 0 .and. coef_ba(2) == 0) then
  phi_0 = 0
else
  phi_0 = atan2 (-coef_ba(1), -coef_ba(2))
endif

if (coef_ba(3) < 0) phi_0 = phi_0 + pi

i_a2 = i_a1 + n_a - 1          ! end of region A

m_min = int((phi_(i_a2) - phi_0)/twopi) + 1
m_max = int((phi_(i_b1) - phi_0)/twopi)

nc = 0
do i = m_min, m_max
  nc = nc + 1
  wave%phi_kick(nc) = twopi * i + phi_0
  do j = i_a2, i_b1
    if (phi_(j) < wave%phi_kick(nc)) wave%ix_cross(nc) = wave%p2%x(j)
  enddo
  if (wave%phi_kick(nc) > phi_ring) &
                           wave%phi_kick(nc) = wave%phi_kick(nc) - phi_ring
  wave%phi_kick(nc) = wave%phi_kick(nc) / 2      ! back to normal phase
enddo
wave%n_cross = nc

phi_ring = phi_ring / 2

! rms to amplitude ratios

if (amp_a == 0 .or. amp_b == 0 .or. amp_ba == 0) then
  wave%chi_a = 0
  wave%rms_a = 0
  wave%rms_b = 0
  wave%rms_k = 0
  wave%rms_phi = 0
else
  wave%chi_a = abs(abs(coef_ba(3)) - amp_ba) / abs(wave%kick(1))
  wave%rms_a = sqrt((rms_a(1)**2 + rms_a(2)**2)) / amp_a
  wave%rms_b = sqrt((rms_b(1)**2 + rms_b(2)**2)) / amp_b
  wave%rms_k = sqrt(rms_ba(1)**2 * coef_ba(1)**2 + &
      rms_ba(2)**2 * coef_ba(2)**2 + rms_ba(3)**2 * coef_ba(3)**2) /  &
      wave%kick(1)**2
  wave%rms_phi = sqrt(rms_ba(2)**2 * coef_ba(1)**2 +  &
                        rms_ba(1)**2 * coef_ba(2)**2) / (2 * amp_ba**2)
endif

! typeout info

c = 180 / pi
print *, 'NUMBER OF POSSIBLE KICKER LOCATIONS:', wave%n_cross
fmt = '(a11, 2i6, 7f12.5)'
print *, '  (Deg)     Ix1   Ix2      Lambda  RMS_Lambda         Rho     RMS_Rho           C       RMS_C         Amp' 
print fmt, 'A Region:', wave%ix_a1, wave%ix_a2, c*coef_a(1), c*rms_a(1), &
                                                c*coef_a(2), c*rms_a(2), c*coef_a(3), c*rms_a(3), c*amp_a
print fmt, 'B Region:', wave%ix_b1, wave%ix_b2, c*coef_b(1), c*rms_b(1), &
                                                c*coef_b(2), c*rms_b(2), c*coef_b(3), c*rms_b(3), c*amp_b
fmt = '(a11, 12x, 5f12.5)'
print fmt, 'B-A:', c*coef_ba(1), c*rms_ba(1), c*coef_ba(2), c*rms_ba(2), c*coef_ba(3), c*rms_ba(3), c*amp_ba

! convert to degrees

do i = 1, plot1%n_use
  plot1%y(i) = plot1%y(i) * 180 / pi
enddo

end subroutine