      ! File: lookup_table_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 LookupTableClass_mod
        !Class Dependency
      
        implicit none
        
        !Class Visibility
        private
        public :: LookupTableClass
        
        !Static Methods
        public :: LookUpTable_GetSpacedData
        public :: LookUpTable_LinInterpList
        public :: LookUpTable_FindIndexAdaptive
        public :: LookUpTable_FindIndexLinear
        public :: LookUpTable_FindIndexBsearch
        public :: LookUpTable_GetXData
        
        !Class
        type , abstract :: LookupTableClass
           
           private
           
           !Public Field
           
           ! 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
           
           !Private Fields
           ! N.B: Protected field must be declared as m_<FieldName>
           
           !Private Field
        contains
        
          !Public Methods
          procedure(ILookUpGetValue) , deferred :: GetValue
          procedure, non_overridable    :: DisposeLookupTable    ! Destructor  method for inherited classes
           
          
        end type LookupTableClass
         
        !Class Constructor Interface (No Constructor Abstract Class) 
        
        !Class Method Interfaces  
        abstract interface
          function ILookUpGetValue(This,x) result(y)
            import LookupTableClass
            implicit none
            class(LookupTableClass) :: this
            double precision, intent(in) :: x(:)
            double precision :: y
          end function          
        end interface
      contains
       
        !Class Constructor (No Constructor Abstract Class)
      
      !Class Destructor
      subroutine DisposeLookupTable(this)
        implicit none
        class(LookupTableClass) :: this
         
      end subroutine
       
      !Class Private Methods
        !> Linearly interpolate the (x, y) input data to the new_x coordinates
        pure function LookUpTable_GetSpacedData(in_x, in_y, new_x) 
     &                                                   result(out_yy)
          implicit none
          double precision, intent(in) :: in_x(:), in_y(:), new_x(:)
          double precision            :: out_yy(size(new_x))
          integer              :: ix
        
          do ix = 1, size(new_x)
             call LookUpTable_LinInterpList
     &                              (in_x, in_y, new_x(ix), out_yy(ix))
          end do
        end function
     
        !> Compute by use of linear interpolation the value in the middle
        ! of a domain D = [x_list(1) , x_list(size(x_list))].
        ! If x_value is left of domain  D,
        ! then the value becomes the value at the left side of D,
        ! if x_value is right of domain D,
        ! then the value becomes the value at the rigth side of D
        pure subroutine LookUpTable_LinInterpList
     &                             (x_list, y_list, x_value, y_value)
          double precision, intent(in)  :: x_list(:), y_list(:)
          double precision, intent(in)  :: x_value
          double precision, intent(out) :: y_value
        
          integer               :: ix, iMin, iMax
          double precision      :: temp
        
          iMin = 1
          iMax = size(x_list)
        
          if (x_value <= x_list(iMin)) then
             y_value = y_list(iMin)
          else if (x_value >= x_list(iMax)) then
             y_value = y_list(iMax)
          else
            ix = LookUpTable_FindIndexAdaptive(x_list, x_value)
            temp = (x_value - x_list(ix-1))/(x_list(ix) - x_list(ix-1))
            y_value = (1.0d0 - temp) * y_list(ix-1) + temp * y_list(ix)
          end if
        end subroutine 
     
        !> Adaptive search (combination of linear and binary search) of sorted list
        !> for the smallest ix such that list(ix) >= val. On failure, returns
        !> size(list)+1
        pure function LookUpTable_FindIndexAdaptive(list, val) 
     &                                                      result(ix)
          double precision, intent(in) :: list(:) !< Sorted list
          double precision, intent(in) :: val     !< Value to search for
          integer              :: ix
          integer, parameter   :: binary_search_limit = 40
        
          if (size(list) < binary_search_limit) then
             ix = LookUpTable_FindIndexLinear(list, val)
          else
             ix = LookUpTable_FindIndexBsearch(list, val)
          end if
        end function  
     
        !> Linear search of sorted list for the smallest ix such that list(ix) >= val.
        !> On failure, returns size(list)+1
        pure function LookUpTable_FindIndexLinear(list, val) result(ix)
          double precision, intent(in) :: list(:) !< Sorted list
          double precision, intent(in) :: val     !< Value to search for
          integer              :: ix
        
          do ix = 1, size(list)
             if (list(ix) >= val) exit
          end do
        end function          
     
        !> Binary search of sorted list for the smallest ix such that list(ix) >= val.
        !> On failure, returns size(list)+1
        pure function LookUpTable_FindIndexBsearch(list, val) result(ix)
          double precision, intent(in) :: list(:) !< Sorted list
          double precision, intent(in) :: val     !< Value to search for
          integer              :: ix, i_min, i_max, i_middle
        
          i_min = 1
          i_max = size(list)
        
          do while (i_min < i_max)
             ! This safely performs: i_middle = (i_max + i_min) / 2
             i_middle = i_min + ishft(i_max - i_min, -1)
        
             if (list(i_middle) >= val) then
                i_max = i_middle
             else
                i_min = i_middle + 1
             end if
          end do
        
          ix = i_min
          if (val > list(ix)) ix = size(list) + 1
        end function
          
        !> Returns the x-coordinates of the lookup table
        pure function LookUpTable_GetXData(x_min, dx, n_points) 
     &                                                   result(xdata)
          double precision, intent(in) :: x_min, dx
          integer, intent(in)  :: n_points
          double precision     :: xdata(n_points)
          integer              :: ix
        
          do ix = 1, n_points
             xdata(ix) = x_min + (ix-1) * dx
          end do
        end function             
       
      end module LookupTableClass_mod
      