;+
; NAME:
;	Planck
; PURPOSE:
;	Calculation of the blackbody differential brightness curve for 
;       a specified temperature and array of wavelengths.
;       This procedure evaluates the Planck function for the spectral
;       photon flux at the specified wavelengths, emitted by a blackbody
;       at the specified temperature, and returns it, multiplied by the
;       user-specified 'factor'.
; CATEGORY:
;	Calibration utility.
; CALLING SEQUENCE:
;	Planck, Wavelength, Parameters, DiffBrightness, PartDeriv
; INPUTS:
;	Wavelength: Array of wavelengths (in Angstroms!).
;       Parameters: Two-element array:
;               (0): 'Temp': the temperature value (K) (assumed scalar);
;               (1): 'Factor': a multiplier to reduce the blackbody 
;                    brightness to the PLBS screen brightness, for example.
; OPTIONAL INPUTS:
;	None.
; KEYWORD PARAMETERS:
;	None.
; OUTPUTS:
;	DiffBrightness: Array of differential brightnesses (kR/A)
;               corresponding to the elements of Wavelength.
; OPTIONAL OUTPUTS:
;	PartDeriv: If the routine is to be used in fitting a blackbody
;               curve to the measured differential brightness of the
;               Polar Camera PLBS, the fitting routine requires the
;               partial derivatives of the function to be fit with respect
;               to the function parameters.  If a named variable is
;               supplied as the fourth parameter to the procedure, it is
;               returned containing an array of partial derivatives of 
;               DiffBrightness w.r.t. 'Temp' and 'Factor';
;               PartDeriv(i,j) = p.d. of DiffBrightness at i'th wavelength
;               w.r.t. j'th variable
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	None.
; RESTRICTIONS:
;	The wavelength array must be expressed in Angstroms.  The
;       differential brightness is expressed in kiloRayleighs per
;       Angstrom.
; PROCEDURE:
;	Straightforward.  The combinations of fundamental constants are
;       telescoped into the numeric parameters k1 and k2.
; EXAMPLE:
;       wl=FIndGen(5)*1000+5000                 ; 5000, 6000, ..., 9000 A
;	Planck,wl,[2000.,1E-12],DiffBr,PDer
; SEE ALSO:
;	
; MODIFICATION HISTORY:
;   Written 92/10 by DPS
;   Modified 93/04 by DPS to be compatible with IDL User Library
;       procedure CURVEFIT
;   Modified 93/09 by DPS to use modified CURVEFIT, called MYCURFIT, in
;       ~/idl.
;   Documented 960605 by DPS.
;-
PRO planck,wl,params,bb,pder
;
; This procedure evaluates the Planck function for the spectral
; photon flux at the specified wavelengths, emitted by a blackbody
; at the specified temperature, and returns it, multiplied by the
; user-specified 'factor'.
;
; INPUTS:
;   wl     = wavelength array (A)
;   params = 2-element array:
;   params(0) = t      = temperature value (K) (assumed scalar)
;   params(1) = factor = multiplier to reduce blackbody brightness to
;                        PLBS screen brightness, for example
; OUTPUT:
;   bb     = array of blackbody output brightnesses (kR/A)
; OPTIONAL OUTPUT:
;   pder   = array of partial derivatives of bb w.r.t. temperature and
;            'factor'; pder(i,j) = p.d. of bb at i'th wavelength
;            w.r.t. j'th variable
;
; MODIFICATION HISTORY:
;   Written 92/10 by DPS
;   Modified 93/04 by DPS to be compatible with IDL User Library
;       procedure CURVEFIT
;   Modified 93/09 by DPS to use modified CURVEFIT, called MYCURFIT, in
;       ~/idl.
;

IF N_Params() LT 3 THEN BEGIN
    Doc_Library,'planck'
    Return
    ENDIF
    
ON_ERROR,2
wl=DOUBLE(wl)
t=DOUBLE(params(0))
factor=DOUBLE(params(1))
k1=1.883651D30
k2=1.438769D08
expon=k2/(wl*t)
denom=((wl^4)*(EXP(expon)-1))
bb=factor*k1/denom
IF N_PARAMS() LE 3 THEN RETURN
; otherwise calculate matrix of partial derivatives as well
PDER=DBLARR(N_ELEMENTS(wl),2)
PDER(0,0)=bb^2*expon*wl^4*EXP(expon)/k1	; p.d. w.r.t. 't'
PDER(0,1)=bb/factor			; p.d. w.r.t. 'factor'
RETURN
END