program potential_to_field

use precision_def
implicit none
real(rp) r,y,x
real(rp) x0,y0,z0,dx,dy,dz
real(rp), allocatable:: Efield_x(:,:), Efield_y(:,:)
real(rp), allocatable :: V(:,:), xx(:), yy(:)
real(rp) volt(-1:1), xvec(-1:1), E
real(rp) V0
real(rp) Ex,Ey,Ez,Bx,By,Bz
integer ix, iy, iz
integer ix_max/0/, ix_min/0/, iy_max/0/, iy_min/0/
integer i,j,k
character*120 string

dx=0.001
dy=0.001
x0=7.112
!dz = dThetaQlong * x0
y0=0.
z0=0.
Ez=0.
Bx=0.
By=0.
Bz=0.
write(12,'(a)')'QUAD_LONG: FULLRING, angle=dThetaQlong, field_calc=grid, &' 
write(12, '(a)')'field = {mode = {grid = {'
write(12, '(a)')'type = xyz,'
write(12,'(a)')'curved_coords = .true.,'
write(12, '(a6,es16.8,a1,es16.8,a1,es16.8,a2)')'r0 = (',x0,',',y0,',',z0,'),'
write(12,'(a6,es16.8,a1,es16.8,a1,a16,a2)')'dr = (',dx,',',dy,',','dThetaQlong * x0','),'
write(12,'(a)')'ele_anchor_pt = BEGINNING,'
do iz = 0,1
 open(unit=11, file='Quadrupole_field_31kV_roundup.dat')
 read(11,'(a)')string
 print *,string
 read(11,'(a)')string
 print *,string
 do while(.true.)
  read(11,'(a)', end=99)string
  call string_trim(string,string,i)
  read(string(1:i),*)x
  call string_trim(string(i+1:),string,i)
  read(string(1:i),*)y
  ix = x/100./dx+sign(0.5_rp,x)
  iy = y/100./dy+sign(0.5_rp,y)
  call string_trim(string(i+1:),string,i)
  read(string(1:i),*)V0
  if(allocated(V))then
      V(ix,iy) = V0
      xx(ix) = x/100.
      yy(iy) = y/100.
  endif
  write(14, '(5es12.4)')x/100.,y/100.,ix*dx, iy*dy,V0
 ix_max = max(ix,ix_max)
 ix_min = min(ix,ix_min)
 iy_max = max(iy,iy_max)
 iy_min = min(iy,iy_min)
 end do
99 close(unit=11)
print '(4i10)',ix_min,ix_max,iy_min,iy_max
if(.not. allocated(V))then
   allocate(V(ix_min:ix_max,iy_min:iy_max))
   allocate(Efield_x(ix_min:ix_max,iy_min:iy_max))
   allocate(Efield_y(ix_min:ix_max,iy_min:iy_max))
   allocate(xx(ix_min:ix_max), yy(iy_min:iy_max))
endif
end do


do iz=0,1
do i=ix_min+1, ix_max-1
  do j=iy_min+1, iy_max-1
   volt(-1:1) = V(i-1:i+1,j)
   call gradient(volt,dx,E)
   Efield_x(i,j) = -E
   Ex=E
   volt(-1:1) = V(i,j-1:j+1)
   call gradient(volt,dx,E)
   Efield_y(i,j) = -E
   Ey=E
   write(13, '(2i10,1x,2es12.4,1x,es12.4)')i,j,i*dx, j*dy, v(i,j)
   write(16,'(2i10,1x,2es12.4,1x,2es12.4)')i,j,i*dx,j*dy, Efield_x(i,j), Efield_y(i,j)
!  write(*,'(a3,1x,3i10,a5,1x,6es22.14,a2)') 'pt(',ix,',',iy,',',iz,') = (', Ex,',',Ey,',',Ez,',',Br,',',By,',',Bz,'),'
if(j==iy_max-1 .and. i == ix_max-1 .and. iz == 1)then
   write(12,'(a3,i5,a1,i5,a1,i5,a5,1x,5(es22.14,a1),es22.14,a4)') 'pt(',i,',',j,',',iz,') = (', Ex,',',Ey,',',Ez,',',Bx,',',By,',',Bz,')}}}'
  else
   write(12,'(a3,i5,a1,i5,a1,i5,a5,1x,5(es22.14,a1),es22.14,a2)') 'pt(',i,',',j,',',iz,') = (', Ex,',',Ey,',',Ez,',',Bx,',',By,',',Bz,'),'
endif
 end do
end do
end do

!write(12,'(2a)')'}}'
end

subroutine gradient(volt,dx,g)
use precision_def
implicit none
real(rp) volt(-1:1),  E
real(rp) a,b, dx,g

b = (volt(1) - volt(-1))/2/dx
a = (volt(1)+volt(-1)-2*volt(0))/2/dx**2

g = b

g=(volt(1)-volt(-1))/2/dx
return
end