! Read the file inf_field_alone.dat that was created by Wuzeng and write in bmad grid field element format
program read_inf_field_write_grid_element
use cesr_utils
use magfield
use magfield_interface, dummy => get_g2_fields
use parameters_bmad

implicit none

   type(magfield_struct), allocatable,save :: map_inj(:,:,:),map_inf(:,:,:)
   type(magfield_struct), allocatable,save :: map(:,:,:)
   type(magfield_struct), save :: xmin_inj, xmin_inf
   type(em_field_struct) field

integer i,n, io
integer ix, ix2, iy
integer nargs
integer ixx, iyy,  ixx_max/0/, iyy_min/0/, iyy_max/0/
integer iz
integer multiplier, fringe/1/, inflector/2/, inflector_or_fringe
integer ind(3)
integer  ind_min(1:3), ind_max(1:3)

character*100 new_file, new_file2, diff_file
character*120 line
character*60 inf_map/'inf_field_alone.dat'/
character*60 fringe_map/'ring_inject_bfield3d_wuzheng_extended.dat'/ !,'injec_fld.dat'
character*60 field_type
character*120 bmad_file

real(rp) deltax/0.001/,deltay/0.001/
real(rp) x0,y0,z0,dx(1:3),r0(1:3),dy,dz
real(rp) Ex,Ey,Ez, B(1:3)
real(rp) xmin(1:3),xmax(1:3)
real(rp) xoffset/718.9/, yoffset/0.0/, zoffset/-430./ !xoffset [cm] of centerline of old inflector, zoffset to end of inflector

field_file(fringe)%name = 'ring_inject_bfield3d_wuzheng_extended.dat'
field_file(fringe)%grid_spacing=0.5
field_file(fringe)%dGdx = 0.
field_file(fringe)%type = 2
field_file(fringe)%flag= 0

field_file(inflector)%name = 'inf_field_alone.dat'
field_file(inflector)%grid_spacing =0.5
field_file(inflector)%dGdx=0.
field_file(inflector)%type =2
field_file(inflector)%flag = 0

  nargs = command_argument_count()
  if (nargs == 1)then
     call get_command_argument(1, field_type)
     if(field_type == 'fringe')inflector_or_fringe = fringe
     if(field_type == 'inflector')inflector_or_fringe = inflector
     print *, 'Using ', field_file(inflector_or_fringe)%name
  endif

!    inflector_element = 'inflector.bmad'

!inflector_or_fringe = fringe

call read_field(inflector_or_fringe,map,xmin_inf)

!write bmad grid element
 ix = index(field_file(inflector_or_fringe)%name,'.dat')
 bmad_file = field_file(inflector_or_fringe)%name(1:ix-1)//'.bmad'
print *,' Write file = ', bmad_file
 open(unit=12, file = bmad_file)
 dx(1:3)=field_file(inflector)%grid_spacing/100.  ! cm to m

r0(1) = -xoffset/100.
r0(2) = yoffset/100.
r0(3) = -zoffset/100.
 B(1:3)=0.

if(inflector_or_fringe == inflector) write(12,'(a)')'INFLECTOR: em_field,  field_calc=fieldmap, l=4.3, &' 
if(inflector_or_fringe == fringe) write(12,'(a)')'FRINGE : em_field, field_calc=fieldmap, l= 4.3, &' 
 write(12, '(a1,a)')'!', field_file(inflector_or_fringe)%name
 write(12, '(a)')'grid_field= {'
 write(12, '(a)')'geometry= xyz,'
 write(12,'(a)')'field_type = magnetic, '
 if(inflector_or_fringe == fringe)write(12,'(a)')'field_scale = fringe_field_scale, '
 if(inflector_or_fringe == inflector)write(12,'(a)')'field_scale = inf_field_scale, '
 write(12,'(a)')'curved_ref_frame = .false.,'
 write(12, '(a6,es16.8,a1,es16.8,a1,es16.8,a2)')'r0 = (',r0(1),',',r0(2),',',r0(3),'),'
 write(12,'(a6,es16.8,a1,es16.8,a1,es16.8,a2)')'dr = (',dx(1),',',dx(2),',',dx(3),'),'
 write(12,'(a)')'ele_anchor_pt = BEGINNING,'

! Here is how the map is laid out
!    x(1:3) = (ix -1) * dx(1:3) + r0(1:3)
!  r0(1:3) = xmin%x(1:3)
   ind(1:3) = map(1,1,1)%x(1:3)/.5
   ind_min(1:3) = map(1,1,1)%x(1:3)/.5
   ind_max(1:3) = map(1,1,1)%x(1:3)/.5
   xmin(1:3) = map(1,1,1)%x(1:3)
   xmax(1:3) = map(1,1,1)%x(1:3)
   do iz=1, size(map, 3)
    do ix = 1,size(map,1)
     do iy = 1,size(map,2)
      ind(1:3) = map(ix,iy,iz)%x(1:3)/.5
      ind_min(1:3) = min(ind_min(1:3), ind(1:3))
      ind_max(1:3) = max(ind_max(1:3), ind(1:3))
      xmin(1:3) = min(xmin(1:3), map(ix,iy,iz)%x(1:3))
      xmax(1:3) = max(xmax(1:3), map(ix,iy,iz)%x(1:3))
      B = map(ix,iy,iz)%B /1.e4  !gauss to tesla
    end do
   end do
  end do
write(12,'(a1,1x,a,1x,3i10,1x,a,1x,3es12.4)')'!',' ind_min = ', ind_min(1:3),' xmin(1:3) = ', xmin(1:3)
write(12,'(a1,1x,a,1x,3i10,1x,a,1x,3es12.4)')'!',' ind_max = ', ind_max(1:3),' xmax(1:3) = ', xmax(1:3)

   do iz=1, size(map,3)
    do ix = 1,size(map,1)
     do iyy = -size(map,2),size(map,2)
      iy = abs(iyy)
      if(iyy == 0 .or. iyy ==-1)cycle
      ind(1:3) = map(ix,iy,iz)%x(1:3)/.5
      B = map(ix,iy,iz)%B /1.e4  !gauss to tesla
      if(iyy < 0)then
       ind(2) = -map(ix,iy,iz)%x(2)/.5
       B(1) = -map(ix,iy,iz)%B(1)/1.e4
      endif
!      B = map(ind(1),ind(2),ind(3))%B /1.e4  !gauss to tesla
      if(iyy==size(map,2) .and. ix == size(map,1) .and. iz == size(map,3))then
!        write(12,'(a3,i5,a1,i5,a1,i5,a5,1x,2(es22.14,a1),es22.14,a2)') 'pt(',ix,',',iy,',',iz,') = (', B(1),',',B(2),',',B(3),')}'
!       else
!        write(12,'(a3,i5,a1,i5,a1,i5,a5,1x,2(es22.14,a1),es22.14,a2)') 'pt(',ix,',',iy,',',iz,') = (', B(1),',',B(2),',',B(3),'),'
        write(12,'(a3,i5,a1,i5,a1,i5,a5,1x,2(es22.14,a1),es22.14,a2)') 'pt(',ind(1),',',ind(2),',',ind(3),') = (', B(1),',',B(2),',',B(3),')}'
       else
        write(12,'(a3,i5,a1,i5,a1,i5,a5,1x,2(es22.14,a1),es22.14,a2)') 'pt(',ind(1),',',ind(2),',',ind(3),') = (', B(1),',',B(2),',',B(3),'),'
      endif
    end do
   end do
  end do

end