!= 内挿関連 プログラム
!
!= interpolation module
!
! Authors::   Masanori Onishi
! Version::   $Id: interpol.f90,v 1.00 onishi$
! Tag Name::  $Name: dcrtm-20170126 $
! Copyright:: Copyright (C) GFD Dennou Club, 2016. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]


module interpol
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !

  ! モジュール引用 ; USE statements
  !


  implicit none

  private

  public :: func_LogLinInterPol
  public :: func_LinearInterPol
  public :: func_BoundGridNum
  public :: PointInterPol
  public :: ArrInterPol

contains 

  function func_zValue( &
    & r_Press1,r_Press2,&
    & r_Temp1, r_Temp2, &
    & r_Val1,  r_Val2)  result(z_Val)
    !
    !calculatez_value
    !
    real(8)::r_Press1,r_Press2
    real(8)::r_Temp1,r_Temp2
    real(8)::r_Val1,r_Val2
    real(8)::z_Val

    z_Val=(r_Press1*r_Val1/r_Temp1+r_Press2*r_Val2/r_Temp2)/&
       &(r_Press1/r_Temp1+r_Press2/r_Temp2)

  endfunction func_zValue

  function func_LogLinInterPol( &
    & Flag, X1, X2, Y1, Y2, XX ) result(YY)

    integer :: Flag
    real(8) :: X1, X2, Y1, Y2, XX, YY

    if(Flag == 1) then ! X-linear, Y-liner

      YY = func_LinearInterPol( X1, X2, Y1, Y2, XX )

    elseif(Flag == 2) then ! X-log, Y-liner

      if( (X1 <= 0) .or. (X2 <= 0) ) then
        print *, 'X1 or X2 <= 0:', X1, X2
        print *, 'stop'
        stop
      end if

      if( XX <= 0 ) then
        print *, 'XX <= 0:', XX
        print *, 'stop'
        stop
      end if

      YY = func_LinearInterPol( log(X1), log(X2), Y1, Y2, log(XX) )

    elseif(Flag == 3) then ! X-linear, Y-log

      if( (Y1 <= 0) .or. (Y2 <= 0) ) then
        print *, 'Y1 or Y2 <= 0:', Y1, Y2
        print *, 'stop'
        stop
      end if

      YY = func_LinearInterPol( X1, X2, log(Y1), log(Y2), XX )
      YY = exp(YY)

    elseif(Flag == 4) then ! X-log, Y-log

      if( (X1 <= 0) .or. (X2 <= 0) ) then
        print *, 'X1 or X2 <= 0:', X1, X2
        print *, 'stop'
        stop
      end if

      if( (Y1 <= 0) .or. (Y2 <= 0) ) then
        print *, 'Y1 or Y2 <= 0:', Y1, Y2
        print *, 'stop'
        stop
      end if

      if( XX <= 0 ) then
        print *, 'XX <= 0:', XX
        print *, 'stop'
        stop
      end if

      YY = func_LinearInterPol( log(X1), log(X2), log(Y1), log(Y2), log(XX))
      YY = exp(YY)

    elseif(Flag == 5) then ! X-linear, Y-log

      if( (Y1 <= 0) .and. (Y2 <= 0) ) then
        YY = 0.0_8
      elseif( (Y1 <= 0) .and. (Y2 > 0) ) then
        YY = Y2
      elseif( (Y1 > 0) .and. (Y2 <= 0) ) then
        YY = Y1
      else

        YY = func_LinearInterPol( X1, X2, log(Y1), log(Y2), XX )
        YY = exp(YY)
      end if

    elseif(Flag == 6) then ! X-log, Y-log

      if( (X1 <= 0) .or. (X2 <= 0) ) then
        print *, 'X1 or X2 <= 0:', X1, X2
        print *, 'stop'
        stop
      end if

      if( XX <= 0 ) then
        print *, 'XX <= 0:', XX
        print *, 'stop'
        stop
      end if

      if( (Y1 <= 0) .and. (Y2 <= 0) ) then
        YY = 0.0_8
      elseif( (Y1 <= 0) .and. (Y2 > 0) ) then
        YY = Y2
      elseif( (Y1 > 0) .and. (Y2 <= 0) ) then
        YY = Y1
      else

        YY = func_LinearInterPol( log(X1), log(X2), log(Y1), log(Y2), log(XX))
        YY = exp(YY)

      end if

    end if

  end function func_LogLinInterPol

  function func_LinearInterPol( x1, x2, y1, y2, x0 ) result(a)

    real(8) :: x1, x2, y1, y2, x0, a

    if ( (x2-x1) .ne. 0.0_8 ) then
      a = ( y2 - y1 ) * ( x0 - x1 ) / ( x2 - x1 ) + y1
    else
      print *, "値が不適切です"
      stop
    end if

  end function func_LinearInterPol

  subroutine PointInterPol(           & ! 2017/01/26 modify
    & Flag, nmax, n_ArrX, n_ArrY, XX, & ! (in )
    & YY                              ) ! (out)
    !
    integer, intent(in ) :: Flag
    integer, intent(in ) :: nmax
    real(8), intent(in ) :: n_ArrX(1:nmax)
    real(8), intent(in ) :: n_ArrY(1:nmax)
    real(8), intent(in ) :: XX
    real(8), intent(out) :: YY
    integer :: ks
    integer :: kt

    ! Flag=1: ascending sequence
    ! Flag=2: descending sequence

    if(nmax < 2) then
      print *, 'dimension of array is smaller than 2:', nmax
      print *, 'stop'
      stop
    end if


    if(Flag == 1) then

      if( XX < n_ArrX(1) ) then
        kt = 2
      elseif( XX > n_ArrX(nmax) ) then
        kt = nmax
      else
        do ks = 2, nmax
          if( (XX >= n_ArrX(ks-1)) .and. (XX <= n_ArrX(ks))) then
            kt = ks
            exit
          end if
        end do
      end if

    elseif(Flag == 2) then

      if( XX > n_ArrX(1) ) then
        kt = 2
      elseif( XX < n_ArrX(nmax) ) then
        kt = nmax
      else
        do ks = 2, nmax
          if( (XX <= n_ArrX(ks-1)) .and. (XX >= n_ArrX(ks))) then
            kt = ks
            exit
          end if
        end do
      end if

    end if

    YY = n_ArrY(ks-1) + &
      & (n_ArrY(ks  ) - n_ArrY(ks-1)) * (XX - n_ArrX(ks-1))/&
      & (n_ArrX(ks  ) - n_ArrX(ks-1))

  end subroutine PointInterPol

  subroutine ArrInterPol(   & ! 2016/11/10
    & Flag,                 & ! (in )
    & amax, a_ArrX, a_ArrY, & ! (in )
    & bmax, b_ArrX,         & ! (in )
    &               b_ArrY  ) ! (out)
    !
    ! Flag=1: ascending sequence
    ! Flag=2: descending sequence
    !
    integer, intent(in ) :: Flag
    integer, intent(in ) :: amax
    real(8), intent(in ) :: a_ArrX(1:amax)
    real(8), intent(in ) :: a_ArrY(1:amax)
    integer, intent(in ) :: bmax
    real(8), intent(in ) :: b_ArrX(1:bmax)
    real(8), intent(out) :: b_ArrY(1:bmax)
    integer :: k_min
    integer :: k_max
    integer :: k_a, k_b
    integer :: num1, num2

    if( amax < 1) then
      print *, 'amax is invalid value (amax >= 1), amax:', amax
      stop
    end if
    if( bmax < 1) then
      print *, 'bmax is invalid value (bmax >= 1), bmax:', bmax
      stop
    end if

    if(Flag == 1) then

      if(a_ArrX(1) <= b_ArrX(1)) then
        k_min = 1
      else

        if(a_ArrX(1) == b_ArrX(bmax)) then
          k_min = bmax
        else
          k_min = func_BoundGridNum(1, bmax, b_ArrX, a_ArrX(1))
        end if

        do k_b = 1, k_min-1
          b_ArrY(k_b) = 0.0_8
        end do
      end if

      if(a_ArrX(amax) > b_ArrX(bmax)) then
        k_max = bmax
      else
        k_max = func_BoundGridNum(1, bmax, b_ArrX, a_ArrX(amax)) - 1
        do k_b = k_max+1, bmax
          b_ArrY(k_b) = 0.0_8
        end do
      end if

      if( (k_max - k_min) > amax ) then

        num1 = k_min
        do k_a = 2, amax
          num2 = func_BoundGridNum(1, bmax, b_ArrX, a_ArrX(k_a)) - 1
          do k_b = num1, num2
            b_ArrY(k_b) = func_LinearInterPol(        &
              &           a_ArrX(k_a-1), a_ArrX(k_a), &
              &           a_ArrY(k_a-1), a_ArrY(k_a), &
              &                          b_ArrX(k_b)  )
          end do
          num1 = num2 + 1
        end do

      else

        do k_b = k_min, k_max
          num1 = func_BoundGridNum(1, amax, a_ArrX, b_ArrX(k_b))
          b_ArrY(k_b) = func_LinearInterPol(          &
            &           a_ArrX(num1-1), a_ArrX(num1), &
            &           a_ArrY(num1-1), a_ArrY(num1), &
            &                           b_ArrX(k_b)  )
        end do

      end if

    elseif(Flag == 2) then

      if(a_ArrX(1) >= b_ArrX(1)) then
        k_min = 1
      else

        if(a_ArrX(1) == b_ArrX(bmax)) then
          k_min = bmax
        else
          k_min = func_BoundGridNum(2, bmax, b_ArrX, a_ArrX(1))
        end if

        do k_b = 1, k_min-1
          b_ArrY(k_b) = 0.0_8
        end do
      end if

      if(a_ArrX(amax) < b_ArrX(bmax)) then
        k_max = bmax
      else
        k_max = func_BoundGridNum(2, bmax, b_ArrX, a_ArrX(amax)) - 1
        do k_b = k_max+1, bmax
          b_ArrY(k_b) = 0.0_8
        end do
      end if

      if( (k_max - k_min) > amax ) then

        num1 = k_min
        do k_a = 2, amax
          num2 = func_BoundGridNum(2, bmax, b_ArrX, a_ArrX(k_a)) - 1
          do k_b = num1, num2
            b_ArrY(k_b) = func_LinearInterPol(        &
              &           a_ArrX(k_a-1), a_ArrX(k_a), &
              &           a_ArrY(k_a-1), a_ArrY(k_a), &
              &                          b_ArrX(k_b)  )
          end do
          num1 = num2 + 1
        end do

      else

        do k_b = k_min, k_max
          num1 = func_BoundGridNum(2, amax, a_ArrX, b_ArrX(k_b))
          b_ArrY(k_b) = func_LinearInterPol(          &
            &           a_ArrX(num1-1), a_ArrX(num1), &
            &           a_ArrY(num1-1), a_ArrY(num1), &
            &                           b_ArrX(k_b)  )
        end do

      end if

    else
      print *, 'Flag is invalid value (Flag = 1 or 2), Flag:', Flag
      stop
    end if


  end subroutine ArrInterPol

  function func_BoundGridNum( &           ! 2016/11/10
    & Flag, nmax, n_Arr, XX ) result(kt)
    !
    ! Flag=1: ascending sequence
    ! Flag=2: descending sequence
    !
    ! Flag=1: g(i-1) < XX <= g(i), return i
    ! Flag=2: g(i-1) > XX >= g(i), return i
    !
    integer :: Flag
    integer :: nmax
    real(8) :: n_Arr(1:nmax)
    real(8) :: XX
    integer :: kt
    integer :: ks


    if(nmax < 1) then
      print *, 'dimension of array is smaller than 1:', nmax
      print *, 'stop'
      stop
    end if


    if(Flag == 1) then

      if( XX < n_Arr(1) ) then
        kt = 1!2  ! 2016/11/08
      elseif( XX >= n_Arr(nmax) ) then
        kt = nmax+1 ! nmax ! 2016/11/08
      !elseif( XX == n_Arr(nmax) ) then
      !  kt = nmax
      else
        do ks = 2, nmax
          if( (XX >= n_Arr(ks-1)) .and. (XX < n_Arr(ks))) then
            kt = ks
            exit
          end if
        end do
      end if

    elseif(Flag == 2) then

      if( XX > n_Arr(1) ) then
        kt = 1!2 ! 2016/11/08
      elseif( XX <= n_Arr(nmax) ) then
        kt = nmax+1 !nmax ! 2016/11/08
      else
        do ks = 2, nmax
          if( (XX <= n_Arr(ks-1)) .and. (XX > n_Arr(ks))) then
            kt = ks
            exit
          end if
        end do
      end if

    else
      print *, 'Flag is invalid value(Flag = 1 or 2), Flag:', Flag
      stop
    end if

  end function func_BoundGridNum

end module interpol
