function BT =  modis_bright(rad, band,units)

% FUNCTION BT =  MODIS_BRIGHT(RAD, BAND, UNITS)
% DESCRIPTION:
%    Compute brightness temperature for an EOS-AM MODIS infrared band.
%
%    Spectral responses for each IR detector were obtained from MCST:
%    ftp://mcstftp.gsfc.nasa.gov/incoming/MCST/PFM_L1B_LUT_4-30-99
%
%    An average spectral response for each infrared band was
%    computed. The band-averaged spectral response data were used
%    to compute the effective central wavenumbers and temperature
%    correction coefficients included in this module.
%
% USAGE:
%    BT = MODIS_BRIGHT(RAD, BAND, UNITS)
%
% INPUT PARAMETERS:
%    RAD           Planck radiance (units are determined by UNITS)
%    BAND          MODIS IR band number (20-25, 27-36)
%    UNITS         Flag defining radiance units
%                  0 => milliWatts per square meter per steradian per
%                       inverse centimeter
%                  1 => Watts per square meter per steradian per micron
%
% OUTPUT PARAMETERS:
%    MODIS_BRIGHT  Brightness temperature (Kelvin)
%                  Note that a value of -1.0 is returned if
%                  BAND is not in range 20-25, 27-36.
%
% MODIFICATION HISTORY:
% Liam.Gumley@ssec.wisc.edu
% http://cimss.ssec.wisc.edu/~gumley
% ...Id: modis_bright.pro,v 1.3 2000/03/15 20:33:50 gumley Exp $
%
% Translated to Matlab by Shaima.Nasiri@ssec.wisc.edu
%                         2000/06/20
%
% Copyright (C) 1999, 2000 Liam E. Gumley
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
%-


%- Check input parameters
 nargs = nargin;
 if (nargs ~= 3) 
   error(['Usage: BT =  MODIS_BRIGHT(RAD, BAND, UNITS)']);
 end
 if isempty(rad)
   error('Argument RAD is undefined');
 end
 if isempty(band)
   error( 'Argument BAND is undefined');
 end
 if isempty(units)
   error('Argument UNITS is undefined');
 end
 if ((band < 20) | (band > 36) | (band == 26) )
   error('Argument BAND must be in the range [20-25, 27-36]');
 end
 

%- BAND-AVERAGED MODIS SPECTRAL RESPONSE FUNCTIONS FOR EOS-AM
%- TEMPERATURE RANGE FOR FIT WAS  180.00 K TO  320.00 K
%- BANDS
%-   20,  21,  22,  23,
%-   24,  25,  26,  27,
%-   28,  29,  30,  31,
%-   32,  33,  34,  35,
%-   36
%- NOTE THAT BAND 26 VALUES ARE SET TO ZERO

%- Effective central wavenumber (inverse centimenters)
 cwn = [...
   2.641775E+03, 2.505277E+03, 2.518028E+03, 2.465428E+03, ...
   2.235815E+03, 2.200346E+03, 0.0,          1.477967E+03, ...
   1.362737E+03, 1.173190E+03, 1.027715E+03, 9.080884E+02, ...
   8.315399E+02, 7.483394E+02, 7.308963E+02, 7.188681E+02, ...
   7.045367E+02];
 
%- Temperature correction slope (no units)
 tcs = [...
   9.993411E-01, 9.998646E-01, 9.998584E-01, 9.998682E-01, ...
   9.998819E-01, 9.998845E-01, 0.0,          9.994877E-01, ...
   9.994918E-01, 9.995495E-01, 9.997398E-01, 9.995608E-01, ...
   9.997256E-01, 9.999160E-01, 9.999167E-01, 9.999191E-01, ...
   9.999281E-01];
 
%- Temperature correction intercept (Kelvin)
 tci = [...
   4.770532E-01, 9.262664E-02, 9.757996E-02, 8.929242E-02, ...
   7.310901E-02, 7.060415E-02, 0.0,          2.204921E-01, ...
   2.046087E-01, 1.599191E-01, 8.253401E-02, 1.302699E-01, ...
   7.181833E-02, 1.972608E-02, 1.913568E-02, 1.817817E-02, ...
  1.583042E-02];

%- Compute brightness temperature
 if (units == 1) 
 %- Radiance units are
 %- Watts per square meter per steradian per micron
  BT = ( bright_m(1.0e+4 / cwn(band - 19), rad) - tci(band - 19) ) ...
       / tcs(band - 19);

 else

 %- Radiance units are
 %- milliWatts per square meter per steradian per wavenumber
  BT = ( brite_m(cwn(band - 19), rad) - tci(band - 19) ) ...
       / tcs(band - 19);

 end

 return

%-------------------------------------------------------------------- 
 
function brite = brite_m(v, r)

% FUNCTION BRITE = BRITE_M(V, R)
%
% DESCRIPTION:
%    Compute brightness temperature given monochromatic Planck radiance.
%
% USAGE:
%    RESULT = BRITE_M(V, R)
%
% INPUT PARAMETERS:
%    V         Wavenumber (inverse centimeters)
%    R         Planck radiance (milliWatts per square meter per
%              steradian per inverse centimeter)
%
% OUTPUT PARAMETERS:
%    BRIGHT_M    Brightness temperature (Kelvin)
%
% MODIFICATION HISTORY:
% Liam.Gumley@ssec.wisc.edu
% http://cimss.ssec.wisc.edu/~gumley
% $Id: modis_bright.m,v 1.1 2001/07/20 16:39:57 kevinb Exp $
% Translated to Matlab by Shaima Nasiri
%                        2000/06/20
%
% Copyright (C) 1999, 2000 Liam E. Gumley
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
%-

%- Constants are from "The Fundamental Physical Constants",
%- Cohen, E. R. and B. N. Taylor, Physics Today, August 1993.

%- Planck constant (Joule second)
 h = 6.6260755e-34;

%- Speed of light in vacuum (meters per second)
 c = 2.9979246e+8;

%- Boltzmann constant (Joules per Kelvin)      
 k = 1.380658e-23;

%- Derived constants      
 c1 = 2.0 * h * c * c;
 c2 = (h * c) / k;

%- Convert wavenumber to inverse meters
 vs = 1.0e+2 * v;
      
%% NOTE BY SWW 13 OCT 00: alog is not a valid matlab function
%% sww 13 dec 00: changing alog to matlab's log (natural logarithm)
%- Compute brightness temperature
% brite =  c2 * vs ./ alog(c1 * vs^3 ./ (1.0e-5 * r) + 1.0);
 brite =  c2 * vs ./ log(c1 * vs^3 ./ (1.0e-5 * r) + 1.0);

 return
 
