c     kspectrum (http://www.meso-star.com/en_Products.html) - This file is part of kspectrum
c     Copyright (C) 2008-2015 - Méso-Star - Vincent Eymet
c     
c     This file must be used under the terms of the CeCILL license.
c     This source file is licensed as described in the file COPYING, which
c     you should have received as part of this distribution.  The terms
c     are also available at
c     http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
c     
      double precision function Planck(T,nu1,nu2)
      implicit none
      double precision T,lambda_c,nu1,nu2,nu_c,B
      double precision lambda1,lambda2,dlambda,dnu
c
c     Purpose: to compute the mean value of the Planck intensity (W/m²/sr/cm¯¹)
c     over a given spectral interval
c
c     Inputs:
c       + T: temperature (K)
c       + nu1: wavenumber at the beginning of the spectral interval (cm¯¹)
c       + nu2: wavenumber at the end of the spectral interval (cm¯¹)
c
c     Output: Planck (W/m²/sr/cm¯¹)
c

      dnu=nu2-nu1
c     Debug
      if (dnu.le.0.0D+0) then
         write(*,*) 'Error from routine Planck:'
         write(*,*) 'dnu=',dnu
         stop
      endif
c     Debug
      call inversecm2micrometer(nu1,lambda1)
      call inversecm2micrometer(nu2,lambda2)
      dlambda=lambda1-lambda2
c     Debug
      if (dlambda.le.0.0D+0) then
         write(*,*) 'Error from routine Planck:'
         write(*,*) 'dlambda=',dlambda
         stop
      endif
c     Debug
      nu_c=(nu1+nu2)/2.0D+0 ! cm¯¹
      call inversecm2micrometer(nu_c,lambda_c)

      Planck=B(T,lambda_c) ! W/m²/sr/µm
      Planck=Planck*dlambda/dnu ! W/m²/sr/cm¯¹

      return
      end




      function B(T,lambda)
      implicit none
c     
c     Planck intensity (W/m²/sr/µm) as function of wavelenght (µm) and temperature (K)
c     
c     Inputs:
c       + T: temperature in K
c       + lambda: wavelenght (µm)
c
c     Output: B (W/m²/sr/µm)
c
      double precision B,T,lambda
      double precision c1,c2,in
      include 'parameters.inc'

      c1=2*pi*hPl*c0**2 ! W.m²
      c2=hPl*c0/kBz*1.D+6 ! µm.K
      in=lambda*T
      if (T.eq.0) then
         B=0.0D+0
      else
         B=c1/(in**5*(dexp(c2/in)-1))*T**5*1.D+24/pi ! W/m²/sr/µm
      endif

      return
      end


      subroutine inversecm2micrometer(nu,lambda)
      implicit none
c     Purpose: to convert a wavenumber (in cm¯¹) to a wavelenght (in µm)
      double precision nu,lambda

      lambda=1.0D+4/nu ! µm

      return
      end


      subroutine micrometer2inversecm(lambda,nu)
      implicit none
c     Purpose: to converta wavelenght (in µm)  to a wavenumber (in cm¯¹) 
      double precision nu,lambda

      nu=1.0D+4/lambda ! cm¯¹

      return
      end