!2013/02/14
!2013/07/10

!flux の計算プログラム

module mod_flux

  implicit none

  real(8) :: PI = 3.1415926535897932_8
                              ! $ \pi $ .
                              ! Circular constant
  real(8) :: LightSpeed = 2.99792458e+8_8
                              ! $ \c $ [m s-1].
                              ! Light Speed
  real(8) :: PlanckConst = 6.62606957e-34_8 !CODATA 2010
                              ! $ \h $ [J s].
  real(8) :: BoltzmannConst = 1.3806488e-23_8 !CODATA 2010
                              ! $ \k_b $ [J K-1]
  real(8) :: SecRadConst = 6.62606957e-34_8 * 2.99792458e+8_8 / 1.3806488e-23_8 !CODATA 2010
                              ! $ \c_2 $ [m K-1]
                              ! Second Radiation Constant (hc/k)
  real(8) :: StB = 5.670373e-8_8 !CODATA 2010
                              ! $ \sigma_{SB} $ [W m-2 K-4].
                              ! ステファンボルツマン定数.
                              ! Stefan-Boltzmann constant

contains

  subroutine calc_Flux(                           &
    & eps, WaveNumMin, WaveNumMax,                & !(in )
    & z_DelOptDep, r_TempIn, r_TempOut, r_FluxIn, & !(in )
    & r_FluxOut                                   & !(out)
    &                                             )
    !
    ! フラックスの計算
    ! Calculate flux
    !

    ! 宣言文 ; Declaration statements
    !
    real(8), intent(in ) :: eps

    real(8), intent(in ) :: WaveNumMin

    real(8), intent(in ) :: WaveNumMax

    real(8), intent(in ) :: z_DelOptDep

    real(8), intent(in ) :: r_TempIn

    real(8), intent(in ) :: r_TempOut

    real(8), intent(in ) :: r_FluxIn

    real(8), intent(out) :: r_FluxOut

    integer :: DivNum
                         ! 分割数
    real(8) :: OptDepH

    real(8) :: err

    real(8) :: pre_r_FluxOut

    integer :: k, kk

    real(8) :: DifFactor = 5.0_8/3.0_8

    !
    r_FluxOut = r_FluxIn
    DivNum = int(z_DelOptDep) + 1
    OptDepH = z_DelOptDep/ DivNum

    do k = 1, DivNum
      r_FluxOut = PI * integral_planck(WaveNumMin, WaveNumMax,                  & !!change
        & r_TempIn - (r_TempIn - r_TempOut) * (2 * k - 1)/(2 * DivNum) ) *      & !!change
        & ( 1.0_8 + exp(-DifFactor * OptDepH) ) * 0.5_8 * DifFactor * OptDepH + &
        & r_FluxOut * exp(-DifFactor * OptDepH)
    end do

    pre_r_FluxOut = r_FluxOut
    ! print *, "first r_FluxOut=", pre_r_FluxOut

    do kk = 1, 1000
      r_FluxOut = r_FluxIn
      DivNum = DivNum * 2
      OptDepH = z_DelOptDep/ DivNum

      do k = 1, DivNum
        r_FluxOut = PI * integral_planck(WaveNumMin, WaveNumMax,                  & !!change
          & r_TempIn - (r_TempIn - r_TempOut) * (2 * k - 1)/(2 * DivNum) ) *      & !!change
          & ( 1.0_8 + exp(-DifFactor * OptDepH) ) * 0.5_8 * DifFactor * OptDepH + &
          & r_FluxOut * exp(-DifFactor * OptDepH)
      end do ! k in kk
      ! print *, "second r_FluxOut=", r_FluxOut
      err = abs( (r_FluxOut - pre_r_FluxOut)/r_FluxOut )
      ! print *, "err=", err
      if( ( err < eps ) .or.                                    &
        & ( log(abs(r_FluxOut)) < minexponent(abs(r_FluxOut)) ) &
        &                                                     ) then

        if( log(abs(r_FluxOut)) < minexponent(abs(r_FluxOut)) ) then
          r_FluxOut = 0.0_8
        end if
        exit
      end if
      pre_r_FluxOut = r_FluxOut
    end do ! kk

  end subroutine calc_Flux

!  subroutine calc_Flux_gray(                           &
!    & eps, z_DelOptDep, r_TempIn, r_TempOut, r_FluxIn, & !(in )
!    & r_FluxOut                                        & !(out)
!    &                                                  )
!    !
!    ! フラックスの計算
!    ! Calculate flux
!    !
!
!    ! 宣言文 ; Declaration statements
!    !
!    real(8), intent(in ) :: eps
!
!    real(8), intent(in ) :: z_DelOptDep
!
!    real(8), intent(in ) :: r_TempIn
!
!    real(8), intent(in ) :: r_TempOut
!
!    real(8), intent(in ) :: r_FluxIn
!
!    real(8), intent(out) :: r_FluxOut
!
!    integer :: DivNum
!                         ! 分割数
!    real(8) :: OptDepH
!
!    real(8) :: err
!
!    real(8) :: pre_r_FluxOut
!
!    integer :: k, kk
!
!    !
!    r_FluxOut = r_FluxIn
!    DivNum = int(z_DelOptDep) + 1
!    OptDepH = z_DelOptDep/ DivNum
!
!    do k = 1, DivNum
!      r_FluxOut = 1.5_8 * StB *                                                 &
!        & ( r_TempIn - (r_TempIn - r_TempOut) * (2 * k - 1)/(2 * DivNum) )**4 * &
!        & ( 1.0_8 + exp(-1.5 * OptDepH) ) * 0.5_8 * OptDepH +                   &
!        & r_FluxOut * exp(-1.5 * OptDepH)
!    end do
!
!    pre_r_FluxOut = r_FluxOut
!    print *, "first r_FluxOut=", pre_r_FluxOut
!
!    do kk = 1, 1000
!      r_FluxOut = r_FluxIn
!      DivNum = DivNum * 2
!      OptDepH = z_DelOptDep/ DivNum
!
!      do k = 1, DivNum
!        r_FluxOut = 1.5_8 * StB *                                                 &
!          & ( r_TempIn - (r_TempIn - r_TempOut) * (2 * k - 1)/(2 * DivNum) )**4 * &
!          & ( 1.0_8 + exp(-1.5 * OptDepH) ) * 0.5_8 * OptDepH +                   &
!          & r_FluxOut * exp(-1.5 * OptDepH)
!      end do ! k in kk
!      print *, "second r_FluxOut=", r_FluxOut
!      err = abs( (r_FluxOut - pre_r_FluxOut)/r_FluxOut )
!      print *, "err=", err
!      if( ( err < eps ) .or.                                    &
!        & ( log(abs(r_FluxOut)) < minexponent(abs(r_FluxOut)) ) &
!        &                                                     ) then
!
!        if( log(abs(r_FluxOut)) < minexponent(abs(r_FluxOut)) ) then
!          r_FluxOut = 0.0_8
!        end if
!        exit
!      end if
!      pre_r_FluxOut = r_FluxOut
!    end do ! kk
!
!  end subroutine calc_Flux_gray

  function integral_planck(WaveNumMin, WaveNumMax, Temp) result(a)
    real(8) :: WaveNumMin, WaveNumMax, Temp
    real(8) :: a
    real(8) :: dWaveNum
    integer :: k_wavenum, k_fpw

    dWaveNum = WaveNumMax - WaveNumMin

    if(dWaveNum <= 1.0_8) then
      a = ( fBk(WaveNumMin, Temp) + fBk(WaveNumMax, Temp) ) * 0.5_8 * &
        & dWaveNum
    else
      k_wavenum = int(dWaveNum) + 1
      dWaveNum = dWaveNum/k_wavenum
      a = ( fBk(WaveNumMin, Temp) + fBk(WaveNumMax, Temp) ) * 0.5_8
      do k_fpw = 1, k_wavenum -1
        a = a + fBk(WaveNumMin + k_fpw * dWaveNum, Temp)
      end do
      a = a * dWaveNum
    end if

  end function integral_planck

!  function fplanck(WaveNum, Temp) result(a)
!    real(8) :: WaveNum, Temp
!    real(8) :: a
!    real(8) :: x
!
!    x = SecRadConst * WaveNum / Temp
!    a = 15.0_8 /(PI * PI * PI * PI) * StB * Temp * Temp * Temp * Temp * &
!      & (x * x * x) / (exp(x) - 1.0_8)
!  end function fplanck

  function fBk(WaveNum, Temp) result(a)
    real(8) :: WaveNum, Temp
    real(8) :: a

    a = 2.0_8 * PlanckConst * LightSpeed * LightSpeed * WaveNum * WaveNum * WaveNum/ &
      & ( exp(SecRadConst * WaveNum/ Temp) - 1.0_8 )

  end function fBk
end module mod_flux
