; IDL script for Radiative Processes Planck function lab


; First define the two Planck functions

FUNCTION planck_lambda, lambda, temp
  ; Planck radiance function with wavelength spectral interval
  ; Inputs: 
  ;   lambda   wavelength in microns
  ;   temp     temperature in Kelvins
  ; Outputs:
  ;   Blambda  returns Planck intensity function value in W/(m^2 ster um)


  RETURN, Blambda
END


FUNCTION planck_nu, nu, temp
  ; Planck radiance function with wavenumber spectral interval
  ; Inputs: 
  ;   nu       wavenumber in inverse centimeters (cm^-1)
  ;   temp     temperature in Kelvins
  ; Outputs:
  ;   Bnu      returns Planck intensity function value in W/(m^2 ster cm^-1)


  RETURN, Bnu
END




;   The following sections are available:
 section='PlanckSpectrum'
; section='IntegratePlanck'



; This script supports 3 output formats: 
;   X windows (out='x'), Postscript ('ps'), Encapsulated Postscript ('eps')
out='x'
psfile=section+"."+out

if (out eq 'ps') then begin
  ;  For Postscript use device fonts (pick Times Roman)
  set_plot, 'ps'
  device, filename=psfile, /color, /portrait, /helvetica, $
   	 xsize=18, ysize=24, xoffset=1.5, yoffset=2.0
  !p.font=0
  !p.charsize=1.2
endif
if (out eq 'eps') then begin
  set_plot, 'ps'
  device, filename=psfile, /color, /portrait, /encapsulated, /helvetica, $
   	xsize=9.0, ysize=12.0, xoffset=1.0, yoffset=1.0
  !p.font=0
  !p.charsize=0.7
endif
if (out eq 'x') then begin
  ;  For X windows use the Hershey stroked fonts
  set_plot, 'x'
  window, xsize=525, ysize=700, title='ATOC/ASTR 5560 Plots'
  !p.font=-1
  !p.charsize=1.2
endif

; Set other global plotting variables
!p.multi = [0,1,2]
!p.thick=1.0  &  !x.thick=1.0  &  !y.thick=1.0
!x.ticklen = 0.013  &  !y.ticklen = 0.010




;  Plot Planck functions vs wavelength and wavenumber
if (section eq 'PlanckSpectrum') then begin
  temps = [300., 280]          ; list of temperatures (K) to plot
  lambdamin=1.0                ; minimum wavelength (micron)
  lambdamax=40.                ; maximum wavelength
  lambdastep=0.05              ; wavelength step size
  Blambdamax=10                ; maximum Planck value (/cm^-1)
  numin=10.                    ; minimum wavenumber (cm^-1)
  numax=2500.                  ; maximum wavenumber 
  nustep=10                    ; wavenumber step size
  Bnumax=0.2                   ; maximum Planck value (/micron)

  ntemp=n_elements(temps)

  ;  First do Planck function vs wavelength
  nlambda = fix( (lambdamax-lambdamin)/lambdastep )
  lambdagrid = lambdamin + lambdastep*findgen(nlambda)

  Blambda = fltarr(n_elements(lambdagrid),ntemp)
  for i = 0, ntemp-1 do $
    Blambda[*,i] = planck_lambda (lambdagrid,temps[i])

  ;  Plot Planck function for wavelength
  plot, lambdagrid, Blambda[*,0], /xstyle, yrange=[0.0,Blambdamax], linestyle=0, $
	xtitle='Wavelength (micron)', ytitle='Planck Radiance (W m!U-2!N sr!U-1!N um!U-1!N)', $
	title='Planck function in wavelength'
  for i = 1, ntemp-1 do $
    oplot, lambdagrid, Blambda[*,i], linestyle=i

  ;  Plot legend
  x0=lambdagrid(nlambda/2)   &  dx=lambdagrid(nlambda-1)-lambdagrid(0)
  x1=x0+0.02*dx  &  x2=x0+0.08*dx      &  x3=x0+0.10*dx
  y0=Blambdamax  &  dy=Blambdamax-0.0  &  ey=0.02*dy
  for i = 0, ntemp-1 do begin
    y1=y0-(0.06*(ntemp-i)-0.01)*dy
    oplot, [x1,x2], [y1,y1], linestyle=i
    xyouts, x3, y1-ey, 'Temperature ='+string(temps[i],format='(F4.0)'), $
	charsize=1.0*!p.charsize
  endfor


  ;  Then do Planck function vs wavenumber
  nnu = fix( (numax-numin)/nustep )
  nugrid = numin + nustep*findgen(nnu)

  Bnu = fltarr(n_elements(nugrid),ntemp)
  for i = 0, ntemp-1 do $
    Bnu[*,i] = planck_nu (nugrid,temps[i])

  ;  Plot Planck function for wavenumber
  plot, nugrid, Bnu[*,0], /xstyle, yrange=[0.0,Bnumax], linestyle=0, $
	xtitle='Wavenumber (cm!U-1!N)', ytitle='Planck Radiance (W m!U-2!N sr!U-1!N cm)', $
	title='Planck function in wavenumber'
  for i = 1, ntemp-1 do $
    oplot, nugrid, Bnu[*,i], linestyle=i

  ;  Plot legend
  x0=nugrid(nnu/2)   &  dx=nugrid(nnu-1)-nugrid(0)
  x1=x0+0.02*dx  &  x2=x0+0.08*dx      &  x3=x0+0.10*dx
  y0=Bnumax  &  dy=Bnumax-0.0  &  ey=0.02*dy
  for i = 0, ntemp-1 do begin
    y1=y0-(0.06*(ntemp-i)-0.01)*dy
    oplot, [x1,x2], [y1,y1], linestyle=i
    xyouts, x3, y1-ey, 'Temperature ='+string(temps[i],format='(F4.0)'), $
	charsize=1.0*!p.charsize
  endfor
endif





; Numerically integrates the Planck function across part of the spectrum
if (section eq 'IntegratePlanck') then begin
  temps = [300., 250.]
  numin=5.                    ; minimum wavenumber (cm^-1)
  numax=100.                 ; maximum wavenumber 
  nustep=5.0                  ; wavenumber step size

  nnu = fix( (numax-numin)/nustep + 0.5)
  nustep = (numax-numin)/nnu
  sum = 0.5*planck_nu (numin,temps)
  sum = sum + 0.5*planck_nu (numax,temps)
  for i = 1, nnu-1 do begin
    nu = numin + i*nustep
    sum = sum + planck_nu (nu,temps)
  endfor
  sum = sum*nustep

  print, 'Planck integral. Wavenumber range: ',numin, numax, nustep
  print, 'Temperatures: ',temps
  print, 'Integrated Planck values: ',sum


  ; Implement the spectral integrated Planck series expansion here:
  c1 = 1.1911E-8
  c2 = 1.4388
  x1 =  c2*numin/temps
  x2 =  c2*numax/temps
  sum1 = 0.0
  sum2 = 0.0
  for n = 1, 3 do begin
    sum1 = sum1 ; + expansion in x1
    sum2 = sum2 ; + expansion in x2 
  endfor
  planck = c1*(temps/c2)^4 * (sum1-sum2)

  print, 'Integrated Planck values: ', planck

endif


if (out ne 'x') then device, /close
end