! Read efield data file generated by Wanwei using Opera as 2-D slice of quads
! The 2-D file the wanwei generated is 'output_designed_27k2_all_cylinder_middle.txt'
! write as bmad grid map style element
program read_write_curvilinear_grid
use precision_def
use cesr_utils
implicit none
integer i,n, io
integer nargs
integer ix, iy, ixx_min/0/, ixx_max/0/, iyy_min/0/, iyy_max/0/
integer iz
character*100 efield_file, label, voltage_file2/''/
character*100 new_file, new_file2, diff_file
character*120 line
character*120 prefix
real(rp), allocatable :: x(:),r(:),z(:)
real(rp), allocatable :: E_radial(:),E_theta(:),E_z(:), Mag_E(:)
real(rp), allocatable :: voltage(:)
real(rp) r0,z_0,E_radial0, E_theta0, E_z0, Mag_E0,voltage0
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, efield_file)
  else
    efield_file = 'output_designed_27k2_all_cylinder_middle.txt'
  endif
  print *, 'Using ', trim(efield_file)


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

io=0
do while(io == 0)
 i=i+1
 read(11,*,iostat=io)r0,z_0,E_radial0, E_theta0, E_z0, Mag_E0,voltage0
 print '(i10,7es12.4)',i,r0,z_0,E_radial0, E_theta0, E_z0, Mag_E0,voltage0
end do
n=i
allocate(r(n),z(n),E_radial(n),E_theta(n),E_z(n), Mag_E(n), voltage(n))
close(11)

open(unit=11, file = efield_file)
print *,' number of lines = ', n
i=0
do i=1,n-1
 read(11,*)r(i),z(i),E_radial(i),E_theta(i),E_z(i), Mag_E(i), voltage(i)
 print '(i10,7es12.4)',i,r(i),z(i),E_radial(i), E_theta(i), E_z(i), Mag_E(i),voltage(i)
end do
close(unit=11)
 
 prefix = efield_file(1:index(efield_file,'.'))
 open(unit=12, file = trim(prefix)//'bmad')
  print *, 'write to ', trim(prefix)//'bmad'
 dx=0.001
 dy=0.001
 x0=0.0
 r0=7.112

 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)')'!', trim(efield_file)
 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 i=1,n-1
!(x, y, z) = dr * (ix, iy, iz) + r0 + r_anchor
      ix = (r(i)/100.-r0)/dx + sign(0.5_rp,r(i)/100.-r0)
      iy = (z(i)/100.-y0)/dy + sign(0.5_rp, z(i)/100.-y0)
      if(i == n-1 .and. iz==1)then
        write(12,'(a3,i5,a1,i5,a1,i5,a5,1x,2(es22.14,a1),es22.14,a2)') 'pt(',ix,',',iy,',',iz,') = (', E_radial(i),',',E_z(i),',',Ez,')}'
       else
        write(12,'(a3,i5,a1,i5,a1,i5,a5,1x,2(es22.14,a1),es22.14,a2)') 'pt(',ix,',',iy,',',iz,') = (', E_radial(i),',',E_z(i),',',Ez,'),'
      endif

   end do
  end do

end