program bpm_test

  use precision_def
  use gain_anal_mod
  use nonlin_bpm_mod

  implicit none

  integer np/500/
  integer nx/10/,ny/100/
  real(rp)  x,y !this is were the bunch is
  real(rp) buttons(4)
  real(rp) dx
  real(rp) sum, xx, yy
  real(rp) kx/25.75/,ky/19.17/
  real(rp) r0/0.1/
  real(rp) d(1:4,0:2,0:2), intensity/1./
  real(rp) nonlin_buttons(4)
  real(rp) z(3)
  real(rp) deltax, deltay, xspan/40./,yspan/20./
  real(rp) x_nonlin(3)

  integer i,j

  nonlin_bpm_use_offsets = .false.
  call nonlin_bpm_set_pointers(11)

  np=600
  dx = 1.0
  sum = 0.
  np=701

 deltax = xspan/nx
 deltay = yspan/ny
 x = -xspan-deltax
 y = -yspan-deltay
 do while(x <= xspan)
  x= x + deltax
y=0
!  y = -yspan-deltay
!   do while(y <= yspan)
!     y=y + deltay

    call bpm_beampipe(x,y,np, r0,buttons(1:4))
!   buttons = 1.
   sum = buttons(1)+buttons(2)+buttons(3)+buttons(4)
   xx= buttons(4)-buttons(3) + buttons(2)-buttons(1)
   yy= buttons(4)-buttons(2) + buttons(3)-buttons(1)
   xx = kx * xx/sum
   yy = ky * yy/sum

    z(1) = x/1000.
    z(2) = y/1000.
    z(3) = 1.
    call nonlin_bpm_interpolate(z,d)
    nonlin_buttons(1:4) = d(1:4,0,0)

    call nonlin_orbit(11,abs(buttons),x_nonlin)

   write(14,'(i,9e12.4)')np, x,y,buttons, xx, yy, r0
   print '(i,8e12.4)' , np, x,y,buttons, xx, yy
   print '(12e12.4)',x,y,buttons, nonlin_buttons, x_nonlin(1:2)
   write(20,'(i,10e12.4)')np, x,y,buttons, nonlin_buttons
   write(21,'(i,12e12.4)')np, x,y,buttons, nonlin_buttons, x_nonlin(1:2)
!   end do
  end do
  

 end