!........................................................................
!
! Subroutine :
!
! Description:
!
! Arguments  :          
!
! Mod/Commons:    
!
! Calls      :      
!
! Author     : 
!
! Modified   :     
! 
!
!    field_file = 'w7-20_15_fid.table'
!
!  write(12,'(6e12.4)') x(i0(3),i0(1),i0(2)),y(i0(3),i0(1),i0(2)),z(i0(3),i0(1),i0(2)), &
!          position
!  write(13,'(6e12.4)')position, &
!                     Bx(i0(1),i0(2),i0(3)),By(i0(1),i0(2),i0(3)),Bz(i0(1),i0(2),i0(3))              
!........................................................................
!
!
! $Log$
! Revision 1.7  2007/01/30 16:15:13  dcs
! merged with branch_bmad_1.
!
! Revision 1.4  2006/01/20 23:37:07  mjf7
! Formatting and linux compatibility changes
!
! Revision 1.3  2003/06/05 18:33:27  cesrulib
! synch with bmad union removal
!
! Revision 1.2  2003/04/30 17:14:50  cesrulib
! dlr's changes since last import
!
! Revision 1.1.1.1  2002/12/13 19:23:28  cesrulib
! import bmadz
!
!
!........................................................................
!

subroutine field_rk_custom (ele,param, s, here, field, field_type)

  use bmad
  use em_field_mod, only: field_interpolate_3d

  implicit none

  type (ele_struct), intent(in) :: ele
  type (coord_struct), intent(in) :: here
  type (lat_param_struct) param

  real(rp), intent(in) ::s
  real(rp), intent(out) ::  field(3)
  real(rp)  position(3), position0(3)
  real(rp)  Z(-500:500,0:20,0:20),X(-500:500,0:20,0:20)
  real(rp)  Y(-500:500,0:20,0:20) 
  real(rp)  BBz(-500:500,0:20,0:20), BBx(-500:500,0:20,0:20) 
  real(rp)  BBy(-500:500,0:20,0:20)
  real(rp)  Bx(-20:20,-20:20,-500:500)
  real(rp)  By(-20:20,-20:20,-500:500)
  real(rp)  Bz(-20:20,-20:20,-500:500)
  real(rp)  half_z
  real(rp)  deltas(3)
  real(rp) r(3), r1(3)

  integer, intent(out) :: field_type
  integer i, j, k
  integer Nx, Ny, Nz
  integer i0(3), i01(3)

  logical :: init_needed = .true.
  character* 70 field_file

  save field_file, deltas, half_z, Bx, By, Bz

!

  field = 0


  if (init_needed) then
    field_file = ele%descrip
!    field_file = 'w7-20_15_fid.table'
    field_type = B_FIELD$ 
    OPEN (1, FILE = field_file, STATUS='OLD')
    read(1,*)Nx, Ny, Nz

    if(Nx <= 0 .or. Ny <= 0 .or. Nz <= 0)then
      print *,' Nx = ',nx,'  ny=',ny,'  Nz = ',nz
      print *,'FIELD_RK_CUSTOM: Bfield data file is not right'
      stop
    endif
    
    do j = 0, Nx - 1
      do k = 0, Ny -1
        do i = -(Nz - 1)/2 , (Nz -1)/2
          read (1,*) Z(i,j,k), X(i,j,k), &
               Y(i,j,k), BBz(i,j,k), BBx(i,j,k), BBy(i,j,k)
        end do
      end do
    end do
    close (1)
    
    X = 0.01 * X !cm to meter
    Y = 0.01 * Y
    Z = 0.01 * Z
    
    half_z = z((Nz-1)/2,0,0) - z(0,0,0)

    do j = 0,Nx-1
      do k = 0, Ny-1
        do i = -(Nz-1)/2, (Nz-1)/2
          Bx(j, k, i) =    BBx(i,j,k)/10. !kg to Tesla and 2->x 3->y 1->z
          Bx(-j, k, i) =  -BBx(i,j,k)/10.  
          Bx(j, -k, i) =  -BBx(i,j,k)/10. ! reverse sign 9-august 2002 to correct symmetry error
          Bx(-j, -k, i) =  BBx(i,j,k)/10. ! reverse sign 9-august 2002 to correct symmetry error
          By(j,k, i) =    BBy(i,j,k)/10.
          By(-j,k, i) =   BBy(i,j,k)/10.
          By(j,-k, i) =   BBy(i,j,k)/10.
          By(-j,-k, i) =  BBy(i,j,k)/10.
          Bz(j,k, i) =    BBz(i,j,k)/10.
          Bz(-j,k, i) =   BBz(i,j,k)/10.
          Bz(j,-k, i) =   -BBz(i,j,k)/10.
          Bz(-j,-k, i) =  -BBz(i,j,k)/10.
          
        end do
      end do
    end do
    
    
    deltas(3) = Z(1,0,0)-Z(0,0,0)
    deltas(1) = X(0,1,0)-X(0,0,0)
    deltas(2) = Y(0,0,1)-Y(0,0,0)
    init_needed = .false.
  endif
  
  field_type = B_FIELD$ 
  !
  
  
  position0(1) = -20*deltas(1)
  position0(2) = -20*deltas(2)
  position0(3) = -500*deltas(3)
  
  position(1) = here%vec(1)
  position(2) = here%vec(3)
  position(3) = s - half_z
  
  r1 = (position-position0)/deltas
  r = (position)/deltas
  i0 = int(r)
  i01 = int(r1)
  !  write(12,'(6e12.4)') x(i0(3),i0(1),i0(2)),y(i0(3),i0(1),i0(2)),z(i0(3),i0(1),i0(2)), &
  !          position
  !  write(13,'(6e12.4)')position, &
  !                     Bx(i0(1),i0(2),i0(3)),By(i0(1),i0(2),i0(3)),Bz(i0(1),i0(2),i0(3))  
  
  field(1) = field_interpolate_3d (position, Bx, deltas, position0) * ele%value(polarity$)
  field(2) = field_interpolate_3d (position, By, deltas, position0) * ele%value(polarity$)
  field(3) = field_interpolate_3d (position, Bz, deltas, position0) * ele%value(polarity$)
  
end subroutine field_rk_custom