MODULE fieldMaps

USE nrtype
IMPLICIT NONE 

TYPE fieldMap2D_B ! Two-dimensional magnetic field, B = ( Bx(x,y), By(x,y) )
  CHARACTER(LEN=72)                     :: filename   ! name of field-map file
  INTEGER                               :: ngx,  ngy  ! number of gridpoints
  REAL(dp), DIMENSION(:), ALLOCATABLE   :: x,    y    ! coordinates
  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Bx,   By   ! field values
  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Bx2,  By2  ! 2nd derivatives (for cubic spline)
END TYPE fieldMap2D_B

TYPE fieldMap2D_E ! Two-dimensional electric field, E = ( Ex(x,y), Ey(x,y) )
  CHARACTER(LEN=72)                     :: filename   ! name of field-map file
  INTEGER                               :: ngx,  ngy  ! number of gridpoints
  REAL(dp), DIMENSION(:), ALLOCATABLE   :: x,    y    ! coordinates
  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ex,   Ey   ! field values
  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ex2,  Ey2  ! 2nd derivatives (for cubic spline)
END TYPE fieldMap2D_E

TYPE fieldMap3D_B ! TODO: 3D Ring+Inflector map
  CHARACTER(LEN=72)                     :: filename         ! name of field-map file
  INTEGER                               :: ngx,  ngy,  ngz  ! number of gridpoints
  REAL(dp), DIMENSION(:), ALLOCATABLE   :: x,    y,    z    ! coordinates
  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Bx,   By,   Bz   ! field values
  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Bx2,  By2,  Bz2  ! 2nd derivatives (for cubic spline) -- more splines?
END TYPE fieldMap3D_B

TYPE fieldMap3D_E ! TODO: 3D Quad maps (i.e. include fringe fields)
  CHARACTER(LEN=72)                     :: filename         ! name of field-map file
  INTEGER                               :: ngx,  ngy,  ngz  ! number of gridpoints
  REAL(dp), DIMENSION(:), ALLOCATABLE   :: x,    y,    z    ! coordinates
  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ex,   Ey,   Ez   ! field values
  REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Ex2,  Ey2,  Ez2  ! 2nd derivatives (for cubic spline) -- more spines?
END TYPE fieldMap3D_E


TYPE (fieldMap2D_B) Kicker
TYPE (fieldMap2D_E) Q1S
TYPE (fieldMap2D_E) Q1L
TYPE (fieldMap2D_E) QStorage


CONTAINS
SUBROUTINE readFieldMap2D( filename, ngx, ngy, x, y, fx, fy, fx2, fy2 )
  IMPLICIT NONE

  CHARACTER(LEN=*), INTENT(IN)                         :: filename  ! field-map filename
  INTEGER,  INTENT(INOUT)                              :: ngx, ngy  ! number of gridpoints
  REAL(dp), INTENT(INOUT), DIMENSION(:),   ALLOCATABLE :: x,   y    ! coordinates
  REAL(dp), INTENT(INOUT), DIMENSION(:,:), ALLOCATABLE :: fx,  fy   ! field values
  REAL(dp), INTENT(INOUT), DIMENSION(:,:), ALLOCATABLE :: fx2, fy2  ! 2nd derivatives

  INTEGER :: i, j ! iterators

  OPEN (UNIT=99, FILE=filename, ACTION='read', STATUS='old')
    READ(99,*) ngx, ngy
    ALLOCATE( x(ngx) )
    ALLOCATE( y(ngy) )
    ALLOCATE( fx(ngx,ngy) )
    ALLOCATE( fy(ngx,ngy) )
    ALLOCATE( fx2(ngx,ngy) )
    ALLOCATE( fy2(ngx,ngy) )
    ! Be mindful of the loop structure here -- make sure it matches your input file!
    DO i=1, ngx
      DO j=1, ngy
        READ(99,*) x(i), y(j), fx(i,j), fy(i,j)
      END DO
    END DO
  CLOSE(UNIT=99)
END SUBROUTINE readFieldMap2D

END MODULE fieldMaps