!= planck 関数計算用モジュール
!
!= planck function calculation module
!
! Authors::   Masanori Onishi
! Version::   $Id: main_dcrtm.f90,v 0.00 onishi$
! Tag Name::  $Name: dcrtm-20150331 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2015. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]

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

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

  ! 物理定数設定2
  ! physical constants setteing 2
  !
  use planetconst, only: LightSpeed,  &
    &                    PlanckConst, &
    &                    SecRadConst


  ! 宣言文 ; Declaration statements
  !
  implicit none

  private

  public :: integral_planck
  public :: fBk
  public :: integral_DBkDT
  public :: fDBkDT

contains


  function integral_planck(WaveNumMin, WaveNumMax, Temp) result(a)
    !
    ! プランク関数の波数積分
    ! wavenuber integral of planck function
    !
    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 fBk(WaveNum, Temp) result(a)
    !
    ! プランク関数 
    ! planck function
    !
    real(8) :: WaveNum, Temp
    real(8) :: a

    if( WaveNum == 0.0_8 ) then ! 2014/06/14
      a = 0.0_8
    else
      a = 2.0_8 * PlanckConst * LightSpeed * LightSpeed * WaveNum * WaveNum * WaveNum/ &
        & ( exp(SecRadConst * WaveNum/ Temp) - 1.0_8 )
    end if

  end function fBk

  function integral_DBkDT(WaveNumMin, WaveNumMax, Temp) result(a)
    ! 
    ! プランク関数の温度微分の積分
    ! integral of dBk/dT
    !
    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 = ( fDBkDT(WaveNumMin, Temp) + fDBkDT(WaveNumMax, Temp) ) * 0.5_8 * &
        & dWaveNum
    else
      k_wavenum = int(dWaveNum) + 1
      dWaveNum = dWaveNum/k_wavenum
      a = ( fDBkDT(WaveNumMin, Temp) + fDBkDT(WaveNumMax, Temp) ) * 0.5_8
      do k_fpw = 1, k_wavenum -1
        a = a + fDBkDT(WaveNumMin + k_fpw * dWaveNum, Temp)
      end do
      a = a * dWaveNum
    end if

  end function integral_DBkDT


  function fDBkDT(WaveNum, Temp) result(a)
    ! 
    ! プランク関数の温度微分
    ! dBk/dT
    !
    real(8) :: WaveNum, Temp
    real(8) :: a

    if( WaveNum == 0.0_8 ) then
      a = 0.0_8 !! really? 2014/09/30 need to check
    else
      a = 2.0_8 * PlanckConst * &
        & LightSpeed * LightSpeed * SecRadConst * &
        & WaveNum * WaveNum * WaveNum * WaveNum * &
        & exp(SecRadConst * WaveNum/ Temp)/ &
        & (Temp * Temp * &
        & (exp(SecRadConst * WaveNum/ Temp) - 1.0_8) * &
        & (exp(SecRadConst * WaveNum/ Temp) - 1.0_8) )
    end if

  end function fDBkDT

end module planckf
