module rad_blackbody_mod
      
      implicit none
      
      contains
      
      subroutine rad_blackbody_planck_law_wavelength(blalong,blat,blae)

      implicit none
      
      double precision, intent(in) :: blalong
      double precision, intent(in) :: blat
      double precision, intent(out) :: blae
      
!      double precision :: sigma,pi,c0,h,cbol,rind,c,c1,c2

      ! physical constants
      double precision, parameter :: sigma=5.670374D-8
      double precision, parameter :: pi=atan(1.d0)*4.d0
      double precision, parameter :: c0=2.997925d+08
      double precision, parameter :: h=6.62607d-34
      double precision, parameter :: cbol=1.380649d-23
      double precision, parameter :: rind=1.d0
      double precision, parameter :: c=c0/rind
      double precision, parameter :: c1=h*(c**2)
      double precision, parameter :: c2=h*c/cbol


      blae=2.d0*pi*c1/blalong**5/(exp(c2/blalong/blat)-1.d0)

      end subroutine rad_blackbody_planck_law_wavelength

      subroutine rad_blackbody_planck_law_wavenumber(blalong,blat,blae)

      implicit none

      double precision, intent(in) :: blalong
      double precision, intent(in) :: blat
      double precision, intent(out) :: blae
      
!      double precision :: sigma,pi,c0,h,cbol,rind,c,c1,c2

      ! physical constants
      double precision, parameter :: sigma=5.670374D-8
      double precision, parameter :: pi=atan(1.d0)*4.d0
      double precision, parameter :: c0=2.997925d+08
      double precision, parameter :: h=6.62607d-34
      double precision, parameter :: cbol=1.380649d-23
      double precision, parameter :: rind=1.d0
      double precision, parameter :: c=c0/rind
      double precision, parameter :: c1=h*(c**2)
      double precision, parameter :: c2=h*c/cbol


      blae=2.d0*pi*c1*blalong**3/(exp(c2*blalong/blat)-1.d0)

      end subroutine rad_blackbody_planck_law_wavenumber
      
      end module rad_blackbody_mod