program phase_amp

   use bmad
   use gain_fit

   implicit none

   integer i
   integer ix
   integer det(4)
   integer j
   integer index
   integer m
   integer detectors(100)
   integer status
 
   logical read_current_from_butns_file
   character*100 name, file, file_fit
   character*1 answer
   real(rp) current
   real(rp) ampx, ampy, phix, phiy

   logical just_one/.false./

  print '(a,$)', ' Real or fake data ? '
  accept *, answer
  call str_upcase(answer,answer)

  index = 0  
 if(answer(1:1) == 'F')then
 
   call get_fake_data

 else  !real data

    call get_turn_by_turn

    just_one = .false.
    if(bpm > 0) just_one = .true.

  endif ! real or fake data

  print '(a)','     bpm       ampx        phix         ampy        phiy'
  write(11, '(a)')'     bpm       ampx        phix         ampy        phiy'
 do m=1,100

 if(.not. just_one) then
     bpm = m
 endif

    if(all(orbit_data(1)%detector(bpm)%button(1:4) == 0))cycle

   call compute_phase_amp(2,0,phix,ampx,phiy,ampy)
   print '(i,4e12.4)',bpm,ampx,phix,ampy,phiy
   write(11,'(i,4e12.4)')bpm,ampx,phix,ampy,phiy

  if(just_one) exit

 end do
 
end

   subroutine compute_phase_amp(nx,ny,phix,ampx,phiy,ampy)

   use gain_fit

   implicit none
   
   real(rp) phix, phiy, ampx, ampy
   real(rp) phase_h, phase_v, x,y,phi,sumx_cos, sumx_sin, sumy_cos, sumy_sin
   real(rp) cx,cy,sx,sy
   real(rp) my_pi/3.141592/

   integer nx, ny
   integer i

   sumx_cos=0.
   sumx_sin=0.
   sumy_cos=0.
   sumy_sin=0.
   do i=1,norbits
     phase_h = orbit_data(i)%detector(bpm)%phase_h
     phase_v = orbit_data(i)%detector(bpm)%phase_h
     x = orbit_data(i)%detector(bpm)%x
     y = orbit_data(i)%detector(bpm)%y
     phi = nx*phase_h + ny*phase_v
     sumx_cos = sumx_cos + x*cos(phi)
     sumx_sin = sumx_sin + x*sin(phi)
     sumy_cos = sumy_cos + y*cos(phi)
     sumy_sin = sumy_sin + y*sin(phi)

!    print '(i,7es12.4)',i,phase_h,phase_v,x,y,phi, sumx_cos, sumx_sin
   end do

   cx = sumx_cos/norbits
   sx = sumx_sin/norbits
   cy = sumy_cos/norbits
   sy = sumy_sin/norbits

   phix = atan2(sx,cx)*180./my_pi
   ampx = 2*sqrt(cx**2+sx**2 )/norbits

   phiy = atan2(sy,cy)*180./my_pi
   ampy = 2*sqrt(cy**2+sy**2 )/norbits

   return

  end