      ! File: lookup_table_5d_mod.f90
      !
      ! Copyright (c) 2019, LEONARDO COMPANY S.p.A.  All rights reserved.
      !
      ! LEONARDO COMPANY S.p.A. and its licensors retain all intellectual property
      ! and proprietary rights in and to this software, related documentation
      ! and any modifications thereto.  Any use, reproduction, disclosure or 
      ! distribution of this software and related documentation,via any medium,
      ! without an express license agreement from LEONARDO COMPANY S.p.A. is strictly 
      ! prohibited.
      !
      ! Written by: Fabrizio Petri <fabrizio.petri@leonardocompany.com>
      ! Creation date:  
      ! 
      module LookupTable5dClass_mod
        !Class Dependency
        use LookupTableClass_mod
        implicit none
        
        !Class Visibility
        private
        public :: LookupTable5dClass
        
        integer , private, parameter :: TABLE_SIZE = 5
        
        !Class
        type , extends(LookupTableClass) :: LookupTable5dClass
           
           private
           
           !Public Fields
           ! N.B: Must be declared as public :: <FieldName>
           
           !Protected Friend Fields
           ! N.B: PROTECTED statement is not supported in Type Bound fields. Protected field must be declared as public :: m_<FieldName>
           !      Programmer MUST NOT modificate these fields outside the ancestor class, or outside inherited child class
           integer                       :: m_NPoints(TABLE_SIZE)  !< The size of the table
           double precision              :: m_XMin(TABLE_SIZE)     !< The minimum lookup coordinate
           double precision              :: m_Dx(TABLE_SIZE)       !< The x-spacing in the lookup coordinate
           double precision              :: m_InvDx(TABLE_SIZE)    !< The inverse x-spacing
           double precision, pointer     :: m_y(:,:,:,:,:)    
           
           !Private Fields
           ! N.B: Protected field must be declared as m_<FieldName>
           
        contains
        
          !Public Methods
          procedure , public  :: GetValue => LT5dGetValue
          procedure , private :: GetCol
          procedure , private :: GetLoc
          procedure , private :: GetColAtLoc
          procedure  :: Dispose => DisposeLookupTable3d        ! Default destructor  method for the current class
          procedure, non_overridable  :: DisposeLookupTable3d    ! Destructor  method for inherited classes
          
          
        end type LookupTable5dClass
         
        !> Type to indicate a location in a 3D lookup table
        type LT5_loc_t
           private
           !> The x-value lies between low_ix and low_ix+1
           integer  :: low_ix(TABLE_SIZE)
           !> The distance from low_ix (up to low_ix+1), given as a real number
           !> between 0 and 1.
           double precision :: low_frac(TABLE_SIZE)
        end type LT5_loc_t        
        
        !Class Constructor Interface  
        public :: New_LookupTable4d
        interface New_LookupTable4d
          procedure LT5d_CreateEmpty
          procedure LT5d_CreateFromData
        end interface
        
        !Class Method Interfaces  
         
      contains
       
        !Class Constructor
        function LT5d_CreateEmpty(xMin,xMax,nPoints) result(NewLT)
          implicit none
          class(LookupTable5dClass) , pointer :: NewLT
          double precision, intent(in) :: xMin(TABLE_SIZE)    !< Minimum coordinate
          double precision, intent(in) :: xMax(TABLE_SIZE)    !< Maximum coordinate
          integer, intent(in)  :: nPoints(TABLE_SIZE) !< How many values to store
          
          allocate(NewLT)
          
          if (any(xMax <= xMin)) 
     &      stop "FATAL.LOOKUP.TABLE.4D.1 create error: x_max <= x_min"
          if (any(nPoints <= 1)) 
     &      stop "FATAL.LOOKUP.TABLE.4D.2 Create error: n_points <= 1"
          
          NewLT%m_nPoints = nPoints
          NewLT%m_xMin    = xMin
          NewLT%m_dx      = (xMax - xMin) / (nPoints - 1)
          NewLT%m_InvDx   = 1.0d0 / NewLT%m_Dx
          
          allocate(NewLT%m_y(
     &          nPoints(1),nPoints(2),nPoints(3),nPoints(4),nPoints(5)))
          NewLT%m_Y = 0.0d0

        end function


        !> This function returns a new lookup table from regularly spaced data
        function LT5d_CreateFromData(xMin,xMax,SpacedData) result(NewLT)
          implicit none
          class(LookupTable5dClass) , pointer :: NewLT          
          double precision, intent(in) :: xMin(TABLE_SIZE)      !< Minimum coordinate
          double precision, intent(in) :: xMax(TABLE_SIZE)      !< Maximum coordinate
          double precision, intent(in), target :: SpacedData(:,:,:,:,:)!< Input data of shape (nx, ny, n_cols)
          integer              :: nPoints(TABLE_SIZE)

        
          nPoints(1) = size(SpacedData, 1)
          nPoints(2) = size(SpacedData, 2)
          nPoints(3) = size(SpacedData, 3)
          nPoints(4) = size(SpacedData, 4)
          nPoints(5) = size(SpacedData, 5)
          
         if (any(xMax <= xMin)) 
     &      stop "FATAL.LOOKUP.TABLE.5D.3 create error: x_max <= x_min"
          if (any(nPoints <= 1)) 
     &      stop "FATAL.LOOKUP.TABLE.5D.4 Create error: n_points <= 1"
        
          NewLT%m_nPoints = nPoints
          NewLT%m_xMin    = xMin
          NewLT%m_Dx      = (xMax - xMin) / (nPoints - 1)
          NewLT%m_InvDx   = 1.0d0 / NewLT%m_Dx

          NewLT%m_Y => SpacedData
        end function        
        
        
        !Class Destructor
        subroutine DisposeLookupTable3d(this)
          implicit none
          class(LookupTable5dClass) :: this
          
          This%m_NPoints = 0
          This%m_XMin    = 0
          This%m_Dx      = 0
          This%m_InvDx   = 0
          if (associated(This%m_Y)) This%m_Y => null()
        end subroutine
         
        !Class Private Methods
        function LT5dGetValue(This,x) result(y)
          implicit none
          class(LookupTable5dClass)   :: This 
          double precision, intent(in) :: x(:)
          double precision :: y
          
          integer n
          
          if (size(x) /= TABLE_SIZE) 
     &      stop "FATAL.LOOKUP.TABLE.5D.5 Get Value: x size /= 4"

          y =  This%GetCol(x(1),x(2),x(3),x(4),x(5))
          
        end function      
         
        !> Get a location in the lookup table
        elemental function GetLoc(This,x1,x2,x3,x4,x5) result(my_loc)
          implicit none
          class(LookupTable5dClass), intent(in) :: This   
          double precision, intent(in) :: x1,x2,x3,x4,x5
          type(LT5_loc_t)         :: my_loc
          double precision        :: frac(TABLE_SIZE)
        
          frac            = ([x1,x2,x3,x4,x5]-This%m_xMin)*This%m_InvDx
          my_loc%Low_ix   = ceiling(frac)
          my_loc%Low_frac = my_loc%Low_ix - frac
        
          ! Check bounds
          where (my_loc%Low_ix < 1)
             my_loc%Low_ix   = 1
             my_loc%Low_frac = 1
          end where
        
          where (my_loc%Low_ix >= This%m_nPoints)
             my_loc%Low_ix   = This%m_nPoints - 1
             my_loc%Low_frac = 0
          end where
        end function    

        !> Get the value of a single column at x
        pure function GetCol(This,x1,x2,x3,x4,x5) result(value)
          implicit none
          class(LookupTable5dClass), intent(in) :: This  
          double precision, intent(in)    :: x1,x2,x3,x4,x5
          double precision                :: value
          type(LT5_loc_t)         :: loc
        
          loc   = This%GetLoc(x1,x2,x3,x4,x5)
          value = This%GetColAtLoc(loc)
          
        end function
        
        !> Get the value of a single column at a location
        pure function GetColAtLoc(This,loc) result(value)
          implicit none
          class(LookupTable5dClass), intent(in) :: This  
          type(LT5_loc_t), intent(in) :: loc
          double precision            :: w(2,2,2)
          double precision            :: y3d(2,2,2)
          double precision            :: y00,y01,y10,y11,value
          double precision            :: W2d(2,2),y2d(2,2) 
          integer i,j,k,m,l

          !Trilinear interpolation coefficients
          w(1, 1, 1) = loc%low_frac(1) * loc%low_frac(2)
          w(2, 1, 1) = (1.0d0 - loc%low_frac(1)) * loc%low_frac(2)
          w(1, 2, 1) = loc%low_frac(1) * (1.0d0 - loc%low_frac(2))
          w(2, 2, 1) = (1.0d0-loc%low_frac(1))*(1.0d0-loc%low_frac(2))
          w(:, :, 2) = w(:, :, 1) * (1.0d0 - loc%low_frac(3))
          w(:, :, 1) = w(:, :, 1) * loc%low_frac(3)
          
          i = loc%low_ix(1)
          j = loc%low_ix(2)
          k = loc%low_ix(3)
          m = loc%low_ix(4)
          l = loc%low_ix(5)
          
          !Can I get better performance prefetching 
          !   w
          !   This%m_Y(i:i+1,j:j+1,k:k+1,m,l)
          !   This%m_Y(i:i+1,j:j+1,k:k+1,m+1,l)
          !   This%m_Y(i:i+1,j:j+1,k:k+1,m,l+1)
          !   This%m_Y(i:i+1,j:j+1,k:k+1,m+1,l+1)
          !
          ! In L2 Cache Through MM_PREFETCH intrinsic subroutine?
          ! Have you got any other suggestion to increase performance of this function? 
          
          !Trilinear interpolation in (/x(1),x(2),x(3),x4Low,x5Low)
          y2d(1,1) = sum(w*This%m_Y(i:i+1,j:j+1,k:k+1,m,l))
          
          !Trilinear interpolation in (/x(1),x(2),x(3),x4High,x5Low)
          y2d(2,1) = sum(w*This%m_Y(i:i+1,j:j+1,k:k+1,m+1,l))
          
          !Trilinear interpolation in (/x(1),x(2),x(3),x4Low,x5High)
          y2d(1,2) = sum(w*This%m_Y(i:i+1,j:j+1,k:k+1,m,l+1))
          
          !Trilinear interpolation in (/x(1),x(2),x(3),x4High,x5High)
          y2d(2,2) = sum(w*This%m_Y(i:i+1,j:j+1,k:k+1,m+1,l+1))         
          
          !Bilinear intepolation
          w2d(1, 1)= loc%Low_frac(4) * loc%Low_frac(5)
          w2d(2, 1)= (1.0d0 - loc%Low_frac(4)) * loc%Low_frac(5)
          w2d(1, 2)= loc%Low_frac(4) * (1.0d0 - loc%Low_frac(5))
          w2d(2, 2)= (1.0d0 - loc%Low_frac(4)) * (1.0d0-loc%Low_frac(5))

          value = sum(w2d*y2d)
        end function
        

      end module LookupTable5dClass_mod
      