c ... a quick hack to test the MAS IR calibration routines

c ... Tree of subprogram calls:
c ... MAIN
c ...     TEMBCK
c ...         INITCF
c ...             GETCOF
c ...                 CENWAV
c ...                 HMWSRF
c ...                     GETLUN
c ...                     INTERP
c ...                 BNDFIT
c ...                     WBRITE
c ...                     WPLANC
c ...                     LINREG
c ...         MASRAD
c ...             PLANCK
c ...         MASBRT
c ...             BRIGHT
c ...     MASCAL
c ...         MASRAD (see above)
c ...         INITCF (see above)

      implicit none
      
      integer band, bb1c( 50 ), bb2c( 50 ), count( 50 ), headc
      real bb1t, bb2t, headt, slope, intercept, rad, brt, masbrt,
     &  temp( 50 ), tembck
      character answer*1
      external masbrt, tembck

c ... the following numbers are from Chris Moellers EMISCAL.f
c ... test input dataset

c ... blackbody 1 counts (for all 50 bands)

      data bb1c / 25 * 0,
     &    2119,  2066,  2359,  2163,  2342,  2042,  2491,  2213,  2140,
     &    2415,  2311,  2243,  2083,  2421,  2070,  2429,  2429,  2436,
     &    2422,  2400,  2436,  2356,  2199,  2297,  6492 /

c ... blackbody 2 counts

      data bb2c / 25 * 0,
     &    2175,  2169,  2619,  2646,  3174,  3312,  4291,  4646,  5322,
     &    6506,  7160,  7646,  8618,  9606,  9712,  9304,  8666, 16578,
     &   24144, 25486, 15082, 23782, 32745, 26221, 19857 /

c ... earth scene counts

      data count / 25 * 0,
     &    2163,  2134,  2511,  2534,  3091,  3149,  3520,  3971,  2454,
     &    2220,  3925,  6263,  6360,  6835,  6159,  4859,  8108, 13568,
     &   23243, 24666, 14549, 21436, 24512,  7883,  3373 /

c ... earth scene temperatures computed in Chris Moellers test program

      data temp / 25 * 0.0,
     &  291.198, 288.018, 285.583, 290.225, 293.159,
     &  292.414, 283.792, 288.217, 257.119, 234.307,
     &  272.331, 288.150, 285.046, 283.416, 280.148,
     &  271.592, 291.064, 285.888, 292.434, 292.600,
     &  292.241, 289.403, 282.485, 258.376, 230.289 /

c ... Data for computing head temperature

      band = 45
      headc = 10612
      bb1t = 243.70
      bb2t = 295.61

c ... Test error handling if required

      write(*,'(''Test error handling (y/n)'')')
      read(*,'(bn,a1)') answer
      if( answer .eq. 'y' .or. answer .eq. 'Y') then
        headc = -1
        bb2c( 31 ) = -1
      endif

c ... Estimate head temperature using band 45 head count,
c ... or band 47 head count if band 45 fails

      headt = tembck( band, bb1c( band ), bb2c( band ), headc, bb1t,
     &  bb2t )
      if( headt .eq. -1.0 ) then
        band = 47
        headc = 10200
        headt = tembck( band, bb1c( band ), bb2c( band ), headc, bb1t,
     &    bb2t )
      endif
        
      write(*,'(''BB1 TEMP (K) = '',f6.2)') bb1t
      write(*,'(''BB2 TEMP (K) = '',f6.2)') bb2t
      write(*,'(''BAND '',i2,'' HEAD COUNT    = '',i5)') band, headc
      write(*,'(''BAND '',i2,'' HEAD TEMP (K) = '',f6.2)') band, headt
      write(*,'(a)') 'BAND   BB1C   BB2C  COUNT    SLOPE       ' //
     &  'INTERCEPT   RADIANCE   TEMP    REFTEMP    DIFF'

c ... Compute earth scene temperature for each band

      do band = 26, 50
            
        call mascal( band, bb1c( band ), bb2c( band ), bb1t, bb2t,
     &    headt, slope, intercept )
        rad = slope * real( count( band ) ) + intercept
        brt = masbrt( band, rad )
        
        write(*,fmt='(1x,i2,1x,3i7,2x,e12.6,2f11.6,3f9.3)')
     &    band, bb1c( band ), bb2c( band ), count( band ),
     &    slope, intercept, rad, brt, temp( band ), brt - temp( band )
      
      end do

      end