c-----------------------------------------------------------------------
c
c  Name:
c      SUBROUTINE BNDFIT
c
c  Purpose:
c      Compute monochromaticity temperature correction coefficients
c      given the spectral response for an infrared band.
c
c  Usage:
c      CALL BNDFIT( N, W, R, WC, T1, T2, TCS, TCI )
c 
c  Input:
c      N       Number of points in wavenumber and spectral response
c              arrays
c      W       Array of wavenumbers where spectral response is defined,
c              typically at 0.1 wavenumber resolution
c              (inverse centimeters)
c      R       Spectral response values (normalized to 1.0 or
c              un-normalized are both acceptable)
c      WC      Effective central wavenumber for the band
c              (see function CENWAV)
c      T1      Lower limit of expected earth scene temperature range (K)
c      T2      Upper limit of expected earth scene temperature range (K)
c
c  Output:
c      TCS     Temperature correction slope (no units)
c      TCI     Temperature correction intercept (K)
c
c  Revised:
c      Liam.Gumley@ssec.wisc.edu
c      $Id: bndfit.f,v 1.4 1999/11/11 20:04:43 gumley Exp $
c
c  Notes:
c      The temperature correction coefficients are applied in the
c      Planck function as 
c
c      R = B( WC, TCS*T+TCI )
c
c      where R is the Planck radiance in the infrared band,
c            WC is the effective central wavenumber,
c            T is the blackbody temperature,
c            B( W, TCS*T+TCI ) is the Planck function.
c
c-----------------------------------------------------------------------

      subroutine bndfit( n, w, r, wc, t1, t2, tcs, tci )

      implicit none

c ... arguments

      integer n
      real w( * ), r( * ), wc, t1, t2, tci, tcs

c ... local variables
      
      integer npts, i, j
      parameter( npts = 100 )
      real r1, r2, rad, te( npts ), ti( npts ), a, b, dw, rx, wx,
     &  wbrite, wplanc
      double precision xbar, ybar, xsig, ysig, yint, slope,
     &  corr, evar, serr
      external wbrite, wplanc

c ... establish 100 point equal spacing in radiance between t1 and t2,
c ... and then convert radiances to brightness temperatures

      r1 = wplanc( wc, t1 )
      r2 = wplanc( wc, t2 )
      do i = 1, npts
        rad = ( r2 - r1 ) * real( i - 1 ) / real( npts - 1 ) + r1
        te( i ) = wbrite( wc, rad )
      end do
      
c ... compute true response weighted radiance for each element of te,
c ... and compute corresponding brightness temperature

      do j = 1, npts
        a = 0.0
        b = 0.0
        do i = 1, n - 1
            dw = w( i + 1 ) - w( i )
            wx = w( i ) + 0.5 * dw
            rx = 0.5 * ( r( i ) + r( i + 1 ) )
            a = a + dw * rx * wplanc( wx, te( j ) )
            b = b + dw * rx
         enddo
         rad = a / b
         ti( j ) = wbrite( wc, rad )
      enddo

c ... compute linear fit for ti vs. te,
c ... and return slope and intercept

      call linreg( npts, te, ti, 
     &  xbar, ybar, xsig, ysig, yint, slope, corr, evar, serr )
      tcs = sngl( slope )
      tci = sngl( yint )
      
      end