!#########################################################################
!		
!    Copyright (C) 2003-2012 Department of Physics and Astronomy,
!                            University of Rochester,
!                            Rochester, NY
!
!    profiles.f90 is part of AstroBEAR.
!
!    AstroBEAR is free software: you can redistribute it and/or modify	  
!    it under the terms of the GNU General Public License as published by 
!    the Free Software Foundation, either version 3 of the License, or    
!    (at your option) any later version.
!
!    AstroBEAR is distributed in the hope that it will be useful, 
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with AstroBEAR.  If not, see <http://www.gnu.org/licenses/>.
!
!#########################################################################
!> tables module

MODULE Tables
   USE GlobalDeclarations
   IMPLICIT NONE

   INTEGER, PARAMETER :: EXTEND_TABLE=0, EXTRAPOLATE_TABLE=1, ZERO_TABLE=2
   TYPE TableDef
      INTEGER :: nPoints(3)=1
      REAL(KIND=qPREC) :: xmin(3), xmax(3), dx(3)
      REAL(KIND=qPREC), POINTER, DIMENSION(:,:,:) :: data
      INTEGER :: bc(3,2) = EXTEND_TABLE
      INTEGER :: dim=1
   END type TableDef

   INTERFACE GetTableValue
      MODULE PROCEDURE GetTableValueScalar, GetTableValueArray
   END INTERFACE

   INTERFACE CreateTable
      MODULE PROCEDURE CreateTableVector, CreateTableScalar
   END INTERFACE

   INTERFACE InTableRange
      MODULE PROCEDURE InTableRangeScalar, InTableRangeVector
   END INTERFACE InTableRange

CONTAINS
   SUBROUTINE CreateTableVector(table, xmin, xmax, dx, minbc, maxbc)
      TYPE(TableDef), POINTER :: table
      REAL(KIND=qPREC), DIMENSION(:) :: xmin, xmax, dx
      INTEGER, DIMENSIOn(:) :: minbc, maxbc
      ALLOCATE(table)
      table%dim=size(xmin)
      table%xmin(1:table%dim)=xmin(1:table%dim)
      table%xmax(1:table%dim)=xmax(1:table%dim)
      table%dx(1:table%dim)=dx(1:table%dim)
      table%bc(1:table%dim,1)=minbc(1:table%dim)
      table%bc(1:table%dim,2)=maxbc(1:table%dim)
      table%nPoints(1:table%dim)=nint((table%xmax(1:table%dim)-table%xmin(1:table%dim))/table%dx(1:table%dim))+1
      ALLOCATE(table%data(table%npoints(1),table%npoints(2),table%npoints(3)))
      table%data=0d0
   END SUBROUTINE CreateTableVector

   SUBROUTINE CreateTableScalar(table, xmin, xmax, dx, minbc, maxbc)
      TYPE(TableDef), POINTER :: Table
      REAL(KIND=qPREC) :: xmin, xmax, dx
      INTEGER :: minbc, maxbc
      ALLOCATE(table)
      table%dim=1
      table%xmin(1)=xmin
      table%xmax(1)=xmax
      table%dx(1)=dx
      table%bc(1,1)=minbc
      table%bc(1,2)=maxbc
      table%nPoints(1)=nint((table%xmax(1)-table%xmin(1))/table%dx(1))+1
      ALLOCATE(table%data(table%npoints(1),1,1))
      table%data=0d0
   END SUBROUTINE CreateTableScalar
   
   

   SUBROUTINE ReadTable(filename, table, iErr)
      TYPE(TableDef), POINTER :: table
      CHARACTER(LEN=*) :: filename
      INTEGER :: i, j, iErr, iErr2
      CHARACTER(LEN=17) :: X
      CHARACTER(LEN=10) :: Y
      IF (ASSOCIATED(table)) RETURN
      ALLOCATE(table)
      OPEN(UNIT=TABLE_DATA_HANDLE, FILE=filename, STATUS='old', IOSTAT=iErr2)
      IF (iErr2 == 0) THEN
         READ(TABLE_DATA_HANDLE, *) table%dim
         IF (table%dim > 3) THEN
            IF (MPI_ID == 0) PRINT*, 'Error tables only supported up to 3 dimensions'
            STOP
         END IF
         DO i=1,table%dim
            READ(TABLE_DATA_HANDLE, *) table%xmin(i), table%xmax(i), table%dx(i), table%bc(i,1), table%bc(i,2)
            table%nPoints(i)=nint((table%xmax(i)-table%xmin(i))/table%dx(i))+1
         END DO
         ALLOCATE(table%data(table%npoints(1),table%npoints(2),table%npoints(3)))
         DO i=1, 9
            IF (table%npoints(3) < 10**i) EXIT
         END DO
         WRITE(Y,'(A5,I1,A4)') '(A1,I',i,',A2)'
         WRITE(X,Y) '(',table%npoints(3),'E)'
!         write(*,*) Y
!         WRITE(*,*) X
         DO i=1, table%npoints(1)
            DO j=1,table%npoints(2)
               READ(TABLE_DATA_HANDLE, *) table%data(i,j,:)
            END DO
         END DO
         CLOSE(TABLE_DATA_HANDLE)
      ELSE
         IF (MPI_ID == 0) WRITE(*,*) 'Error loading table from ', filename
         iErr=iErr2
      END IF

   END SUBROUTINE ReadTable


   SUBROUTINE WriteTable(filename, table)
      TYPE(TableDef) :: table
      CHARACTER(LEN=*) :: filename
      INTEGER :: i, j
      CHARACTER(LEN=17) :: X
      CHARACTER(LEN=10) :: Y
      OPEN(UNIT=TABLE_DATA_HANDLE, FILE=filename, STATUS='unknown')
      WRITE(TABLE_DATA_HANDLE, *) table%dim
      DO i=1,table%dim
         WRITE(TABLE_DATA_HANDLE, '(3E26.15,2I5)') table%xmin(i), table%xmax(i), table%dx(i), table%bc(i,1), table%bc(i,2)
      END DO
      DO i=1, 9
         IF (table%npoints(3) < 10**i) EXIT
      END DO
      WRITE(Y,'(A5,I1,A4)') '(A1,I',i,',A2)'
      WRITE(X,Y) '(',table%npoints(3),'E)'
      DO i=1, table%npoints(1)
         DO j=1,table%npoints(2)
            WRITE(TABLE_DATA_HANDLE, X) table%data(i,j,:)
         END DO
      END DO
      CLOSE(TABLE_DATA_HANDLE)
   END SUBROUTINE WriteTable


   FUNCTION InTableRangeScalar(table, x)
      TYPE(TableDef) :: table
      REAL(KIND=qPREC) :: x
      LOGICAL :: InTableRangeScalar
      InTableRangeScalar=x >= table%xmin(1) .AND. x <= table%xmax(1)
   END FUNCTION InTableRangeScalar

   FUNCTION InTableRangeVector(table, x)
      TYPE(TableDef) :: table
      REAL(KIND=qPREC), DIMENSION(:) :: x
      LOGICAL :: InTableRangeVector
      InTableRangeVector=ALL(x >= table%xmin(1:size(x)) .AND. x <= table%xmax(1:size(x)))
   END FUNCTION InTableRangeVector

   FUNCTION GetTableValueScalar(table, x)
      TYPE(TableDef) :: table
      REAL(KIND=qPREC) :: x, GetTableValueScalar
      IF (table%dim == 1) THEN
         GetTableValueScalar=GetTableValueArray(table, (/x/))
      ELSE
         IF (MPI_ID == 0) PRINT*, 'Error - GetTableValue requires multiple components for multi-dimensional arrays'
         STOP
      END IF
   END FUNCTION GetTableValueScalar

   FUNCTION GetTableValueArray(table, x)
      TYPE(TableDef) :: table
      REAL(KIND=qPREC) :: x(:), GetTableValueArray
      INTEGER :: i1(1:3,0:1), l(3), m(3), dim, i, j, k
      REAL(KIND=qPrec) :: ix(3), values(0:1,0:1,0:1)
      dim=table%dim
      ix=0d0
      ix(1:dim)=(x(1:dim)-table%xmin(1:dim))/table%dx(1:dim)+1d0
      !Goal is to build coefficients for (a1*x+b1)(a2*x+b2)...(an*x+b2)
      !or to successively to linear interpolation of data points...
      
      i1=1
      i1(1:dim,0)=min(max(floor(ix(1:dim)), 1), table%nPoints(1:dim)-1)
      i1(1:dim,1)=i1(1:dim,0)+1 !max(min(ceiling(ix(1:dim)), table%npoints(1:dim)), i1(1:dim,0)+1)
      l=0
      DO i=0,2**dim
         DO j=1,dim
            l(j)=mod(i/2**(j-1),2)
         END DO
         values(l(1),l(2),l(3))=table%data(i1(1,l(1)),i1(2,l(2)),i1(3,l(3)))
      END DO
      DO k=dim,1,-1 !do z on 4 points, y on 2 points, x on 1 point
         DO i=0,2**(k-1)-1
            l=0
            DO j=1,k-1
               l(j)=mod(i/2**(j-1),2)
            END DO
            l(k)=0
            m=l
            m(k)=1
            CALL interp(values(l(1),l(2),l(3)),values(m(1),m(2),m(3)), ix(k)-REAL(i1(k,0)), table%bc(k,:)) !replaces 1st value with interpolated value
         END DO
      END DO
      GetTableValueArray=values(0,0,0)

   END FUNCTION GetTableValueArray


   SUBROUTINE Interp(y1, y2, dx, bc)
      REAL(KIND=qPREC) :: y1, y2, dx
      INTEGER :: bc(2)
      IF (dx >= 0d0 .AND. dx <=1d0) THEN
         y1=y1+(y2-y1)*dx
      ELSEIF (dx < 0d0) THEN
         SELECT CASE(bc(1))
         CASE(ZERO_TABLE)
            y1=0d0
         CASE(EXTEND_TABLE)
            y1=y1
         CASE(EXTRAPOLATE_TABLE)
            y1=y1+(y2-y1)*dx
         END SELECT
      ELSEIF (dx > 1d0) THEN
         SELECT CASE(bc(2))
         CASE(ZERO_TABLE)
            y1=0d0
         CASE(EXTEND_TABLE)
            y1=y2
         CASE(EXTRAPOLATE_TABLE)
            y1=y1+(y2-y1)*dx
         END SELECT
      END IF
   END SUBROUTINE Interp

   SUBROUTINE DestroyTable
      TYPE(TableDef) :: Table
      IF (ASSOCIATED(table%data)) DEALLOCATE(table%data)
   END SUBROUTINE DestroyTable
END MODULE Tables
