;******************************************************************************
;+
;*NAME:
;
;	BBCMP     30-MAR-83
;  
;*CLASS:
;  
;*CATEGORY:
;  
;*PURPOSE:
;
;	To compute the rms and average deviations of an observed spectrum from
;	Planck curves calculated at user-specified temperatures.
;  
;*CALLING SEQUENCE:
;
;	BBCMP,WAVE,FLUX,TEMP,AVGD,RMS
;  
;*PARAMETERS:
;
;    	WAVE	(REQ) (I) (0 1) (B I L F D)
;		Wavelength vector for the spectrum.
;
;    	FLUX	(REQ) (I) (0 1) (B I L F D)
;		Flux vector for the spectrum.
;
;    	TEMP	(REQ) (I) (0 1) (B I L F D)
;		Scalar or vector of temperature of comparison blackbody 
;		curves.
;
;    	AVGD	(REQ) (O) (0 1) (F D)
;		Scalar or vector of average deviations of observed spectrum .
;               from blackbody curves.
;
;    	RMS	(REQ) (O) (0 1) (F D)
;		Scalar or array of root-mean-square deviation of spectrum .
;               from optimally scaled blackbody of temperature TEMP.
;  
;*INTERACTIVE INPUT:
;
;	none
;  
;*FILES USED:
;
;	none
;  
;*SYSTEM VARIABLES USED:
;
;        !NOPRINT
;  
;*SUBROUTINES CALLED:
;
;	PCHECK
;	PARCHECK
;	PLANCK
;  
;*SIDE EFFECTS:
;  
;*RESTRICTIONS:
;  
;*NOTES:
;
;       The blackbody curves are calculated on the wavelength grid of WAVE.
;	The blackbody fluxes are scaled by
;
;		BBFLUX = BBFLUX*TOTAL(FLUX*BBFLUX)/TOTAL(BBFLUX*BBFLUX)
;
;       before the average and rms deviations are calculated.
;
;*PROCEDURE:
;
;       Input arrays are verified by PCHECK. WAVE and FLUX must be of the same
;	size and type.  If TEMP is scalar, it is forced to be a 1-element
;	floating-point array.  PLANCK returns black body fluxes which are
;	normalized and then scaled as described above.  AVGD and RMS are
;	calculated as follows:
;
;       		F = input flux array
;       		BBF = scaled black body flux array
;       		DIFF = BBF - F
;       then
;       		AVGD = TOTAL(ABS(DIFF))/NELEMENTS(F)-1
;       		RMS = SQRT(TOTAL(DIFF*DIFF)/NELEMENTS(F)-1
;
;       Reference:  Bevington, P.R. 1969, Data Reduction and Error Analysis for
;		    the Physical Sciences, McGraw-Hill (New York), p. 14- 19.
;  
;*INF_1:
;  
;*EXAMPLES:
;
;       To test several temperatures:
;
;		bbcmp,w,f,[20000,40000,60000],avgd,rms
;  
;*MODIFICATION HISTORY:
;
;       Programmer: R.J. Panek 30 March 1983
;        7-30-86 CAG corrected error in rms deviation calculation, added
;		     average deviation calculation, and scaled BBFLUX by its
;		     maximum value to avoid floating point overflow errors.
;       12-17-86 RWT change print format, use !MODE for !VAR5, and remove
;		     FORMAT for VAX conversion
;        4-13-87 RWT VAX mods: add PARCHECK
;        7-24-87 RWT use !NOPRINT
;        4-26-88 RWT add procedure call listing and remove USAGE
;        7-18-90 RWT update prolog & fix error when multiple temperatures are
;		     specified.
;        3-04-91 PJL mofified for unix/sun
;        6-19-91 PJL cleaned up; tested on SUN and VAX; updated prolog
;       11-27-91 PJL corrected typo in prolog
;	13 Sep 93  PJL  allow wave and flux to be of different data types, but
;			they still must have the same number of elements
;
;-
;******************************************************************************
 pro bbcmp,wave,flux,temp,avgd,rms
;     
;  check the data types for the inputs 
;
 npar = n_params(0)
 if npar eq 0 then begin
    print,'BBCMP,WAVE,FLUX,TEMP,AVGD,RMS'
    retall
 endif  ; npar
 parcheck,npar,5,'BBCMP'
 pcheck,wave,1,110,0111
 pcheck,flux,2,110,0111
 pcheck,temp,3,110,0111
 s1 = size(wave) 
 s2 = size(flux)
;
;  check that wave and flux have the same number of data points
;
 if (s1(s1(0)) ne s2(s2(0))) then begin  ; bad input
    print,'Wave and flux must have the same number of elements.'
    retall
 endif  ; bad input
;
;  find out how many temperatures are to be compared
;
 s = size(temp)
 if (s(0) ne 1) then tt = fltarr(1) + temp else tt = float(temp)  ; force vector
 ntmp = n_elements(tt)
 rms = fltarr(ntmp)   ; define RMS array
 avgd = rms           ; define average deviation array
;
;  Calculate the blackbody curve for each temperature. Scale the 
;  curve to the maximum in the allowed wavelength range to avoid
;  floating point overflow errors when the flux weighting is done
;
 for i=0,ntmp-1 do begin
    planck,tt(i),wave,bbflux
    bbflux = bbflux/max(bbflux) ; scale for numerical purposes
;
;  flux weighted scaling
;
    bbflux = bbflux * total(flux*bbflux)/total(bbflux*bbflux)
    diff = (bbflux-flux)
;
;  use definition of rms deviation from Bevington p. 19
;  also add average deviation for further comparisons
;
    avgd(i) = total(abs(diff))/(s2(1)-1)
    rms(i) = sqrt( total( (diff*diff)>1.e-36)/(s2(1)-1))
 endfor  ; i
;
;  use system variable !noprint to control the optional tabular printout
;
 if (!noprint eq 0) then begin  ; printout
    print,'  temp     rms error     avg deviation'
    for i=0,ntmp-1 do print,"$(1x,f8.1,2e14.5)",tt(i),rms(i),avgd(i)
 endif  ; printout
;
 return
 end  ; bbcmp