program read_quad_grid
use precision_def
use cesr_utils
implicit none
integer i,n, io
integer ix, ix2
integer nargs
integer ixx, iyy, ixx_min/0/, ixx_max/0/, iyy_min/0/, iyy_max/0/
integer iz
character*100 voltage_file, label, voltage_file2/''/
character*100 new_file, new_file2, diff_file
character*120 line
real(rp), allocatable :: x1(:),y1(:),V1(:)
real(rp), allocatable :: x2(:),y2(:),V2(:), delta_v(:)
real(rp), allocatable :: DV(:,:)
real(rp) a,b,c
real(rp) deltax/0.001/,deltay/0.001/
real(rp) dv_dx, dv_dy
real(rp) x0,y0,z0,dx,dy,dz
real(rp) Ex,Ey,Ez,Bx,By,Bz


  nargs = command_argument_count()
  if (nargs == 1)then
     call get_command_argument(1, voltage_file)
     print *, 'Using ', trim(voltage_file)
  elseif(nargs == 2)then
     call get_command_argument(1, voltage_file)
     call get_command_argument(2, voltage_file2)
     print *, 'Using ', trim(voltage_file),' and ', trim(voltage_file2)

  else

    voltage_file = 'Quad_field_32kV.dat'
    print '(a,$)',' Lattice file name ? (default= Quad_field_32kV.dat) '
    read(5,'(a)') line
    call string_trim(line, line, ix)
    voltage_file = line
    if(ix == 0) voltage_file = 'Quad_field_32kV.dat'
    print *, ' voltage_file = ', voltage_file
  endif



open(unit=11, file = voltage_file)
i=0

read(11,'(a)')label
print *,' label = ', label
io=0
do while(io == 0)
 i=i+1
 read(11,*,iostat=io)a,b,c
! print '(i10,3es12.4)',i,a,b,c
end do
n=i
allocate(x1(n),y1(n),V1(n), delta_v(n))
close(11)

open(unit=11, file = voltage_file)
read(11,'(a)')label
print *,' label = ', label
i=0
do i=1,n-1
 read(11,*)x1(i),y1(i),V1(i)
! print '(i10,3es12.4)',i,x1(i),y1(i),V1(i)
end do
close(unit=11)

ix=index(voltage_file,'.dat')
new_file = voltage_file(1:ix-1)//'_efield.dat'
open(unit=12, file=new_file)
write(12,'(a)')label
do i=1,n-1
  delta_v(i) = v1(i)
  write(12,'(3es12.4)')x1(i),y1(i),V1(i)
end do
print *, new_file,' written'
close(unit=12)


if(voltage_file2 /= '')then
 allocate(x2(n),y2(n),V2(n))
 open(unit=11, file = voltage_file2)
 read(11,'(a)')label
 print *,' label = ', label
 i=0
 do i=1,n-1
  read(11,*)x2(i),y2(i),V2(i)
!  print '(i10,3es12.4)',i,x2(i),y2(i),V2(i)
 end do

 ix2=index(voltage_file2,'.dat')
 new_file = voltage_file2(1:ix2-1)//'-efield.dat'
 open(unit=12, file=new_file)
 write(12,'(a)')label
 do i=1,n-1
   write(12,'(3es12.4)')x2(i),y2(i),V2(i)
 end do
 print *,new_file,' written'
 close(unit=12)

 open(unit=12, file=voltage_file(1:ix-1)//'-'//voltage_file2(1:ix2-1)//'.dat')
 write(12,'(a)')label
 do i=1,n-1
   delta_v(i) = V1(i)-V2(i)
   write(12,'(5es12.4,2i10)')x1(i),x2(i),y1(i),y2(i),V1(i)-V2(i)
 end do
 print *, voltage_file(1:ix-1)//'-'//voltage_file2(1:ix2-1)//'.dat', ' written'
 close(unit=12)


endif

! take the gradient

do i = 1,n-1
 ixx = x1(i)*10+sign(0.5_rp,x1(i))
 iyy = y1(i)*10+sign(0.5_rp,y1(i))
 ixx_min = min(ixx,ixx_min)
 ixx_max = max(ixx,ixx_max)
 iyy_min = min(iyy,iyy_min)
 iyy_max = max(iyy,iyy_max)
 print '(2i10,2es12.4)',ixx,iyy,x1(i)*10, y1(i)*10
end do
print *,'ixx_min = ',ixx_min,' ixx_max = ',ixx_max, ' iyy_min = ', iyy_min,' iyy_max = ', iyy_max
allocate(DV(ixx_min:ixx_max,iyy_min:iyy_max))

do i = 1,n-1
 ixx = x1(i)*10+sign(0.5_rp,x1(i))
 iyy = y1(i)*10+sign(0.5_rp,y1(i))
 DV(ixx,iyy) = delta_v(i)
end do

 new_file=voltage_file(1:ix-1)//'-efield.dat'
 if(voltage_file2 /= '')new_file=voltage_file(1:ix-1)//'-'//voltage_file2(1:ix2-1)//'-efield.dat'
print *, new_file
open(unit=12, file=new_file)
  do ixx = ixx_min+1, ixx_max-1
   do iyy = iyy_min+1, iyy_max-1
    dv_dx = -(DV(ixx+1,iyy)-DV(ixx-1,iyy))/2/deltax
    dv_dy = -(DV(ixx,iyy+1)-DV(ixx,iyy-1))/2/deltay
    write(12,'(2i10,2es12.4)') ixx,iyy,dv_dx,dv_dy
  end do
end do

print *, new_file, ' written'
close(unit=12)


!write bmad grid element
 open(unit=12, file = 'quad_one.dat')
 dx=0.001
 dy=0.001
 x0=0.0

 y0=0.
 z0=0.
 Ez=0.
 Bx=0.
 By=0.
 Bz=0.
! write(12,'(a)')'QUAD_LONG: sbend, rho=7.112, angle=dThetaQlong, field_calc=fieldmap, &' 
 write(12, '(a1,a,a3,a)')'!', trim(voltage_file),' - ',trim(voltage_file2)
 write(12, '(a)')'grid_field= {'
 write(12, '(a)')'geometry= xyz,'
 write(12,'(a)')'field_type = electric, '
 write(12,'(a)')'field_scale = -1., '
 write(12,'(a)')'curved_ref_frame = .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
    do ixx = ixx_min+1, ixx_max-1
     do iyy = iyy_min+1, iyy_max-1
      Ex = -(DV(ixx+1,iyy)-DV(ixx-1,iyy))/2/deltax
      Ey = -(DV(ixx,iyy+1)-DV(ixx,iyy-1))/2/deltay
 
      if(iyy==iyy_max-1 .and. ixx == ixx_max-1 .and. iz == 1)then
        write(12,'(a3,i5,a1,i5,a1,i5,a5,1x,2(es22.14,a1),es22.14,a2)') 'pt(',ixx,',',iyy,',',iz,') = (', Ex,',',Ey,',',Ez,')}'
       else
        write(12,'(a3,i5,a1,i5,a1,i5,a5,1x,2(es22.14,a1),es22.14,a2)') 'pt(',ixx,',',iyy,',',iz,') = (', Ex,',',Ey,',',Ez,'),'
      endif
    end do
   end do
  end do

end