# ------------------------------------------------------------
# SRFI.sml
# ------------------------------------------------------------
# SET WARNING LEVEL: Refer to B5 & B6.
  $warnings 3
# ------------------------------------------------------------
# DEFINE PROCEDURE writeTitle: Refer to B12 & B13.
# PURPOSE: WRITES TITLE & AUTHOR INFO TO CONSOLE WINDOW.
  proc writeTitle() begin
     printf("SRFI.sml:\n");
     printf("     SML VERSION: November 15, 2005\n");
     printf(" TNTmips VERSION: Released 7.1\n"); 
     printf("         PURPOSE: CONVERTS MS DNs ");
     printf("TO SRFI VALUES.\n");
     printf("     MAY PRODUCE: PVI & PBI RASTERS\n");
     printf("         DETAILS: FAQs_by_Jack A & B\n");
     printf("          AUTHOR: Dr. Jack F. Paris\n");
     printf("    CONTACT INFO: jparis37@msn.com ");
     printf(" 303-775-1195\n");
     printf("     ALLOWED USE: ONLY NON-COMMERCIAL\n\n");
  end
# ------------------------------------------------------------
# DEFINE VARIABLES FOR GENERAL CHARACTER STRINGS: Refer to B7.
  string t$,p$,p1$,p2$,p3$,p4$,p5$,p6$,p7$,p8$,p9$;
  string p10$,p11$,p12$,p13$,p14$,p15$,p16$,p17$,p18$,p19$;
  string p20$,p21$;
# ------------------------------------------------------------
# CLEAR CONSOLE WINDOW & ASK USER TO RE-POSITION IT:
# Refer to B8 & B9.
  clear();
  p1$ = "CONSOLE-WINDOW ADJUSTMENT:\n";
  p2$ = "* REPOSITION the CONSOLE WINDOW.\n";
  p3$ = "* Then, CLICK OK.";
  p$  = p1$ + p2$ + p3$;
  PopupMessage(p$);
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO DATETIME & CONTRAST CLASSES:
# Refer to B10 - B11.
  class DATETIME cdate$;
  class CONTRAST smlContrast;
  string clut1$,clut2$,clut3$;
# ------------------------------------------------------------
# ASSIGN VALUES TO VARIABLES RELATED TO smlContrast:
  clut1$  = "exp"; clut2$  = "Exponential";
  clut3$  = "Exponential contrast table";
  smlContrast.Power            =   0.8;
  smlContrast.InputLowerLimit  =   1;
  smlContrast.OutputLowerLimit =   1;
  smlContrast.OutputUpperLimit = 255;
# ------------------------------------------------------------
# WRITE TITLE & AUTHOR INFORMATION: Refer to B13.
  writeTitle();
# ------------------------------------------------------------
# DECLARE USER-INPUT VARIABLES: Refer to B14 & B15.
  string  site$,gc$;
  numeric cdate,pdate,dd,d1,d2;
  numeric sunelevang,imager,vnir;
  numeric atcor,delcf,msfac,icRL;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO DATA TYPES OF RASTERS:
# Refer to B14 & B15.
  string dntype$,srfitype$,ptype$;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO USER-INPUT VARIABLES:
# Refer to B14 & B15.
  numeric he,doy,sc,dnlow,dnhigh;
  string imager$,atcor$,gi$;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO SRFapc-to-SRFsfc CONVERSION:
# Refer to B14 & B15.
  numeric cCB,cBL,cGL,cYL,cRL,cRE,cNA,cNB;
  numeric cMA,cMB,cMC,cMD,cME,cMF,cMG,cRLm1,pc;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO BOOLEAN ENABLERS: 
# Refer to B14 & B15.
  numeric pF,pCB,pBL,pYL,pRE,pNA,pNB;
  numeric pMA,pMB,pMC,pMD,pME,pMF,pMG,nc;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO WAVELENGTH-BASED MODELS:
# Refer to B14 & B15.
  numeric p,p2,wLenCB,wLenBL,wLenGL,wLenYL,wLenRL,wLenRE;
  numeric wLenNA,wLenNB,wLenMA,wLenMB,wLenMC,wLenMD,wLenME;
  numeric wLenMF,wLenMG;
  numeric wLenV,logwLenRLdwLenV,wtBL;
  numeric srfpathV,logsrfpathV;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO DN-to-SR CONVERSION:
# Refer to B14 & B15.
  numeric skCB,skBL,skGL,skYL,skRL,skRE,skNA,skNB;
  numeric skMA,skMB,skMC,skMD,skME,skMF,skMG;
  numeric dnbCB,dnbBL,dnbGL,dnbYL,dnbRL,dnbRE,dnbNA,dnbNB;
  numeric dnbMA,dnbMB,dnbMC,dnbMD,dnbME,dnbMF,dnbMG;
# ------------------------------------------------------------
# DECLARE DIRECT-SOLAR-SPECTRAL-IRRADIANCE VARIABLES:
# Refer to B14 & B15.
  numeric dssiCB,dssiBL,dssiGL,dssiYL,dssiRL,dssiRE;
  numeric dssiNA,dssiNB,dssiMA,dssiMB,dssiMC,dssiMD,dssiME;
  numeric dssiMF,dssiMG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO PVI & PBI: Refer to B14 & B15.
  numeric pfac,angrad,sinang,cosang;
  numeric pbioff,pvioff,srfiNRint,bslineslope,maxpvipbi;
  numeric srfiNR,tsrfiNR,pbif,pbi,pvif,pvi;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO ARRAY PARAMETERS: 
# Refer to B14 & B15.
  numeric maxi,maxip1;
# ------------------------------------------------------------
# DECLARE LIST OF POSSIBLE INPUT RASTERS: Refer to B14 & B15.
  raster CB,BL,GL,YL,RL,RE,NA,NB,MA,MB,MC,MD,ME,MF,MG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO LOOPS: Refer to B14 & B15.
  numeric i,nlins,ncols;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO INPUT RASTERS: 
# Refer to B14 & B15.
  numeric dnCB,dnBL,dnGL,dnYL,dnRL,dnRE,dnNA,dnNB;
  numeric dnMA,dnMB,dnMC,dnMD,dnME,dnMF,dnMG;
  numeric dnminCB,dnminBL,dnminGL,dnminYL,dnminRL,dnminRE;
  numeric dnminNA,dnminNB,dnminMA,dnminMB,dnminMC,dnminMD;
  numeric dnminME,dnminMF,dnminMG;
  numeric dnmaxCB,dnmaxBL,dnmaxGL,dnmaxYL,dnmaxRL,dnmaxRE;
  numeric dnmaxNA,dnmaxNB,dnmaxMA,dnmaxMB,dnmaxMC,dnmaxMD;
  numeric dnmaxME,dnmaxMF,dnmaxMG;
# ------------------------------------------------------------
# DECLARE VARIABLES FOR HISTOGRAM-EDGE FINDING ALGORITHM:
# Refer to B14 & B15.
  numeric dnheCB,dnheBL,dnheGL,dnheYL,dnheRL,dnheRE;
  numeric dnheNA,dnheNB,dnheMA,dnheMB,dnheMC,dnheMD;
  numeric dnheME,dnheMF,dnheMG;
  numeric logsrfpathRL,difoflogs;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO SItoa & SRFtoa: 
# Refer to B14 & B15.
  numeric pi100,esd,sitoafac;
  numeric sitoaCB,sitoaBL,sitoaGL,sitoaYL,sitoaRL,sitoaRE;
  numeric sitoaNA,sitoaNB,sitoaMA,sitoaMB,sitoaMC,sitoaMD;
  numeric sitoaME,sitoaMF,sitoaMG;
  numeric aCB,aBL,aGL,aYL,aRL,aRE,aNA,aNB;
  numeric aMA,aMB,aMC,aMD,aME,aMF,aMG;
  numeric mCB,mBL,mGL,mYL,mRL,mRE,mNA,mNB;
  numeric mMA,mMB,mMC,mMD,mME,mMF,mMG;
  numeric srftoaCB,srftoaBL,srftoaGL,srftoaYL,srftoaRL;
  numeric srftoaRE,srftoaNA,srftoaNB,srftoaMA,srftoaMB;
  numeric srftoaMC,srftoaMD,srftoaME,srftoaMF,srftoaMG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO DNpath: Refer to B14 & B15.
  numeric dnpathCB1,dnpathBL1,dnpathGL1,dnpathYL1,dnpathRE1;
  numeric dnpathNA1,dnpathNB1,dnpathMA1,dnpathMB1,dnpathMC1;
  numeric dnpathMD1,dnpathME1,dnpathMF1,dnpathMG1;
  numeric dnpathCB,dnpathBL,dnpathGL,dnpathYL,dnpathRL;
  numeric dnpathRE,dnpathNA,dnpathNB,dnpathMA,dnpathMB;
  numeric dnpathMC,dnpathMD,dnpathME,dnpathMF,dnpathMG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO SRFpath: Refer to B14 & B15.
  numeric srfpathCB1,srfpathBL1,srfpathGL1,srfpathYL1;
  numeric srfpathRE1,srfpathNA1,srfpathNB1,srfpathMA1;
  numeric srfpathMB1,srfpathMC1,srfpathMD1,srfpathME1;
  numeric srfpathMF1,srfpathMG1;
  numeric srfpathCB2,srfpathBL2,srfpathGL2,srfpathYL2;
  numeric srfpathRE2,srfpathNA2,srfpathNB2,srfpathMA2;
  numeric srfpathMB2,srfpathMC2,srfpathMD2,srfpathME2;
  numeric srfpathMF2,srfpathMG2;
  numeric srfpathCB,srfpathBL,srfpathGL,srfpathYL,srfpathRL;
  numeric srfpathRE,srfpathNA,srfpathNB,srfpathMA,srfpathMB;
  numeric srfpathMC,srfpathMD,srfpathME,srfpathMF,srfpathMG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO SRFapc: Refer to B14 & B15.
  numeric srfapcCB,srfapcBL,srfapcGL,srfapcYL,srfapcRL;
  numeric srfapcRE,srfapcNA,srfapcNB,srfapcMA,srfapcMB;
  numeric srfapcMC,srfapcMD,srfapcME,srfapcMF,srfapcMG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO SRFsfc & SRFI:
# Refer to B14 & B15.
  numeric srfsfcCB,srfsfcBL,srfsfcGL,srfsfcYL,srfsfcRL;
  numeric srfsfcRE,srfsfcNA,srfsfcNB,srfsfcMA,srfsfcMB;
  numeric srfsfcMC,srfsfcMD,srfsfcME,srfsfcMF,srfsfcMG;
  numeric srfiCB,srfiBL,srfiGL,srfiYL,srfiRL,srfiRE,srfiNA;
  numeric srfiNB,srfiMA,srfiMB,srfiMC,srfiMD,srfiME,srfiMF;
  numeric srfiMG;
# ------------------------------------------------------------
# DECLARE LIST OF POSSIBLE OUTPUT RASTERS: Refer to B14 & B15.
  raster SRFICB,SRFIBL,SRFIGL,SRFIYL,SRFIRL,SRFIRE,SRFINA;
  raster SRFINB,SRFIMA,SRFIMB,SRFIMC,SRFIMD,SRFIME,SRFIMF;
  raster SRFIMG,PBI,PVI;
# ------------------------------------------------------------
# GET NAME OF SITE FROM THE USER: Refer to B16.
  t$  = "SITE NAME";
  p1$ = "SITE-NAME:\n";
  p2$ = "* ENTER SITE NAME.\n";
  p3$ = "* Then, CLICK OK.\n";
  p4$ = "NAME ENTERED:";
  p$  = p1$ + p2$ + p3$ + p4$;
  site$ = PopupString(p$,"NAME",t$);
  printf("       SITE NAME: %s\n",site$);
# ------------------------------------------------------------
# GET COLLECTION DATE FROM THE USER: Refer to B17.
  p1$ = "COLLECTION-DATE:\n";
  p2$ = "  FORMAT: YYYYMMDD\n";
  p3$ = "* ACCEPT the Default DATE,\n";
  p4$ = "* Or, ENTER the Correct DATE.\n";
  p5$ = "* Then, Click OK.\n\n";
  p6$ = "DATE ENTERED:";
  p$  = p1$ + p2$ + p3$ + p4$ + p5$ + p6$;
  dd = 20010509; d1 = 19710723; d2 = 29991231;
  cdate = PopupNum(p$,dd,d1,d2,0);
  dd = cdate; d1 = cdate;
  cdate$.SetDateYYYYMMDD(cdate);
  doy = cdate$.GetDayOfYear();
  printf(" COLLECTION DATE: %8d\n",cdate);
  printf(" DAY OF THE YEAR: %d\n",doy);
# ------------------------------------------------------------
# GET SUN ELEVATION ANGLE FROM THE USER: Refer to B18.
  p1$ = "SUN-ELEVATION-ANGLE:\n";
  p2$ = "   UNITS: deg above the Horizon\n";
  p3$ = "  FORMAT: NN.NN\n";
  p4$ = "   RANGE: 1.00  to 90.00 deg\n";
  p5$ = "* ACCEPT the Default ANGLE,\n";
  p6$ = "* Or, ENTER the Correct ANGLE.\n";
  p7$ = "* Then, CLICK OK.\n\n";
  p8$ = "ANGLE ENTERED:";
  p$  = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
  sunelevang = PopupNum(p$,66.99,1,90,2);
  printf(" SUN ELEV. ANGLE: %5.2f deg\n",sunelevang);
# ------------------------------------------------------------
# GET IMAGER NUMBER FROM THE USER: Refer to B19. !!!!!!!!!!!!!!!!!! to mozna skrocic !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  p1$  = "IMAGER-NUMBER:\n";
  p2$  = "  IMAGER\n";
  p3$  = "  NUMBER:  NAME OF IMAGER      STATUS\n";
  p4$  = "       1:  World View II     Available\n";
  p5$  = "       2:  Ikonos MS        Available\n";
  p6$  = "       3:  OrbView-3 MS           Not\n";
  p7$  = "       4:  Landsat-7 ETM+   Available\n";
  p8$  = "       5:  Landsat-5 TM     Available\n";
  p9$  = "       6:  Landsat-5 MSS          Not\n";
  p10$ = "       7:  Landsat-4 TM     Available\n";
  p11$ = "       8:  Landsat-4 MSS          Not\n";
  p12$ = "       9:  Landsat-3 MSS          Not\n";
  p13$ = "      10:  Landsat-2 MSS          Not\n";
  p14$ = "      11:  Landsat-1 MSS          Not\n";
  p15$ = "      12:  Terra ASTER      Available\n";
  p16$ = "      13:  Terra MODIS            Not\n";
  p17$ = "      14:  Aqua MODIS             Not\n";
  p18$ = "* ACCEPT the Default IMAGER NUMBER,\n";
  p19$ = "* Or, SELECT the Correct IMAGER NUMBER.\n";
  p20$ = "* Then, CLICK OK.\n\n";
  p21$ = "NUMBER ENTERED:";
  p$   = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$ + p9$;
  p$   = p$ + p10$ + p11$ + p12$ + p13$ + p14$ + p15$ + p16$;
  p$   = p$ + p17$ + p18$ + p19$ + p20$ + p21$;
  imager = PopupNum(p$,12,14,1,0);  
# ------------------------------------------------------------
# GET ATMOSPHERIC CORRECTION (atcor) LEVEL FROM THE USER:
# Refer to B20.
  p1$ = "ATMOSPHERIC-CORRECTION LEVEL:\n";
  p2$ = "     LEVEL: TYPE OF ATMOSPHERIC CORRECTION\n";
  p3$ = "         1: No Atmospheric Correction\n";
  p4$ = "         2: Atmospheric Reflectance Only\n";
  p5$ = "         3: All Atmospheric Effects\n";
  p6$ = "* ACCEPT the Default LEVEL,\n";
  p7$ = "* Or, SELECT a Differernt LEVEL.\n";
  p8$ = "* Then, CLICK OK.\n\n";
  p9$ = "LEVEL ENTERED:";
  p$  = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$ + p9$;
  atcor = PopupNum(p$,3,1,3,0);
# ------------------------------------------------------------
# WRITE CORRECTION LEVEL TO THE CONSOLE WINDOW:
# Refer to B20.
  printf("CORRECTION LEVEL: ");
  p4$ = "No Atmospheric Corrections";
  p5$ = "Corrected for Atmospheric Reflectance";
  p6$ = "Corrected for All Atmospheric Effects";
  if (atcor==1) then printf("%s\n",p4$);
  if (atcor==2) then printf("%s\n",p5$);
  if (atcor==3) then printf("%s\n",p6$);
# ------------------------------------------------------------
# IF atcor > 1, GET delcf PARAMETER FROM USER:
# Refer to B22.
  if (atcor > 1) then begin
     p1$ = "HISTOGRAM-EDGE FINDING PARAMETER:\n";
     p2$ = "  FORMAT: N.NN (Units: Percentage Points)\n";
     p3$ = "   RANGE: 0.01 to 1.00 Percentage Points\n";
     p4$ = "* Either, ACCEPT the Default VALUE,\n";
     p5$ = "* Or, ENTER a Better VALUE.\n";
     p6$ = "* Then, CLICK OK.\n\n";
     p7$ = "VALUE ENTERED:";
     p$  = p1$ + p2$ +p3$ + p4$ + p5$ + p6$ + p7$;
     delcf = PopupNum(p$,0.05,0.01,1.00,2);
     printf("           delcf: %4.2f Percentage",delcf);
     printf(" Points\n");
  end
  he = 0.01 * delcf;
# ------------------------------------------------------------
# IF atcor = 3, GET the msfac VALUE: Refer B23.
  if (atcor == 3) then begin
     p1$ = "msfac-VALUE ENTRY:\n";
     p2$ = "   msfac: Scaling Factor for All MS Bands.\n";
     p3$ = "  FORMAT: N.NNNN\n";
     p4$ = "   RANGE: 0.1000 to 2.0000\n";
     p5$ = "* ACCEPT the Default VALUE,\n";
     p6$ = "* Or, ENTER a Better VALUE.\n";
     p7$ = "  NOTE: Don't Change msfac Unless Justified.\n";
     p8$ = "* Then, CLICK OK.\n\n";
     p9$ = "VALUE ENTERED:";
     p$  = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
     p$ = p$ + p9$;
     msfac = PopupNum(p$,1,0.1,2,4);
     printf("           msfac: %6.4f\n",msfac);
# ------------------------------------------------------------
# IF atcor = 3, GET icRL VALUE FROM THE USER:
# Refer to B24.
     p1$ = "icRL-VALUE ENTRY:\n";
     p2$ = "    icRL: Input c-Factor for the RL Band.\n";
     p3$ = "  FORMAT: N.NNNN\n";
     p4$ = "   RANGE: 1.0000 to 3.0000\n";
     p5$ = "* ACCEPT the Default Value,\n";
     p6$ = "* Or, ENTER a Better Value.\n";
     p7$ = "* Then, CLICK OK.\n\n";
     p8$ = "VALUE ENTERED:";
     p$  = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
     icRL = PopupNum(p$,1.34,1,2,4);
     printf("            icRL: %6.4f\n",icRL);
     cRL = icRL;
  end
# ------------------------------------------------------------
# CALCULATE cRLm1 & ASSIGN THE pc VALUE:
# Refer to B26.
  cRLm1 = cRL - 1; pc = 2.27;
# ------------------------------------------------------------
# GENERATE IMAGER-SPECIFIC PARAMETERS: Refer to B25. !!!!!!!!!!!! Modyfications to World View 2 spectral channels !!!!!!!!!!!!!!!!!!!!
# besed on the information published on www.digitalglobe.com
  if (imager == 1) then begin
     imager$ = "QuickBird MS";
     pCB=1;pBL=1; pGL=1; pYL=1; pRL=1; pRE=1; pNA=1; pNB=1;
	 maxi=2100; dnlow=1; dnhigh=2047;
     wLenCB=0.427; wLenBL=0.478; wLenGL=0.546; wLenYL=0.608 ; wLenRL=0.659; wLenRE=0.724; wLenNA=0.831; wLenNB=0.908;
     #skBL=0.2359; skGL=0.1453; skRL=0.1785; skNA=0.1353;
     skCB=0.18591308; skBL=0.2101375; skGL=0.13875816; skYL=0.145745375; skRL=0.18393717; skRE=0.1297034; skNA=0.0979504; skNB=0.05023463;
 
	 dnbCB=0; dnbBL=0; dnbGL=0; dnbYL=0; dnbRL=0; dnbRE=0; dnbNA=0; dnbNB=0;
	 dssiCB= 1710; dssiBL=1925; dssiGL=1843; dssiYL= 1740; dssiRL=1575; dssiRE=1325; dssiNA=1056; dssiNB=870;
     #dssiBL=1925; dssiGL=1843; dssiRL=1575; dssiNA=1113;
  end
# dssiCB= 1710; dssiBL=1925; dssiGL=1843; dssiYL= 1740; dssiRL=1575; dssiRE=1325; dssiNA=1056; dssiNB=870;
#skCB=0.18591308; skBL=0.2101375; skGL=0.13875816; skYL=0.145745375; skRL=0.18393717; skRE=0.1297034; skNA=0.0979504; skNB=0.05023463;
# ------------------------------------------------------------
# CALCULATE C-FACTORS & FOR ENABLED MS BANDS: Refer to B26.
  if (pCB) then cCB = 1 + cRLm1 * (wLenRL/wLenCB)^pc;
  if (pBL) then cBL = 1 + cRLm1 * (wLenRL/wLenBL)^pc;
  cGL = 1 + cRLm1 * (wLenRL/wLenGL)^pc;
  if (pYL) then cYL = 1 + cRLm1 * (wLenRL/wLenYL)^pc;
  pc = 1.1;
  if (pRE) then cRE = 1 + cRLm1 * (wLenRL/wLenRE)^pc;
  if (pNA) then cNA = 1 + cRLm1 * (wLenRL/wLenNA)^pc;
  if (pNB) then cNB = 1 + cRLm1 * (wLenRL/wLenNB)^pc;
  if (pMA) then cMA = 1 + cRLm1 * (wLenRL/wLenMA)^pc;
  if (pMB) then cMB = 1 + cRLm1 * (wLenRL/wLenMB)^pc;
  if (pMC) then cMC = 1 + cRLm1 * (wLenRL/wLenMC)^pc;
  if (pMD) then cMD = 1 + cRLm1 * (wLenRL/wLenMD)^pc;
  if (pME) then begin
     cME = 1 + cRLm1 * (wLenRL/wLenME)^pc;
     cME = cME * 1.1;
  end
  if (pMF) then begin
     cMF = 1 + cRLm1 * (wLenRL/wLenMF)^pc;
     cMF = cMF * 1.5;
  end
  if (pMG) then begin
     cMG = 1 + cRLm1 * (wLenRL/wLenMG)^pc;
     cMG = cMG * 2.4;
  end
# ------------------------------------------------------------
# PARAMETERS RELATED TO WVI & WBI CALCULATIONS: Refer to B36.
  srfiNRint = 254; bslineslope = 1.086;
  pbioff = 100; pvioff = 1000; 
  pfac = 0.2723659; maxpvipbi = 3000;
# ------------------------------------------------------------
# IF atcor < 3, SET C-FACTORS EQUAL TO 1: Refer to B26.
  if (atcor < 3) then begin
     cCB = 1; cBL = 1; cGL = 1; cYL = 1; cRL = 1; cRE = 1;
     cNA = 1; cNB = 1; cMA = 1; cMB = 1; cMC = 1; cMD = 1;
     cME = 1; cMF = 1; cMG = 1;
  end
# ------------------------------------------------------------
# SET UP ARRAYS: Refer to B22, B28 & B29.
  maxip1 = maxi + 1;
  array numeric vCB[maxip1]; array numeric vBL[maxip1];
  array numeric vGL[maxip1]; array numeric vYL[maxip1];
  array numeric vRL[maxip1]; array numeric vRE[maxip1];
  array numeric vNA[maxip1]; array numeric vNB[maxip1];
  array numeric vMA[maxip1]; array numeric vMB[maxip1];
  array numeric vMC[maxip1]; array numeric vMD[maxip1];
  array numeric vME[maxip1]; array numeric vMF[maxip1];
  array numeric vMG[maxip1]; array numeric   v[maxip1];
# ------------------------------------------------------------
# DEFINE proc checkHisto(Band)
# PURPOSE: Check that Input Raster has a Full (Unsampled)
#          Histogram.  If Not, Then Recompute Histogram.
  proc checkHisto (raster Band) begin
     local numeric sampleInt;
     sampleInt = HistogramGetSampleInterval(Band);
     if (sampleInt > 1) begin
        DeleteHistogram(Band); CreateHistogram(Band,0);
     end
     else if (sampleInt == -1) then	begin
        CreateHistogram(Band,0);
     end
  end
# ------------------------------------------------------------
# DEFFINE func dnHE(heloc): Refer to B22 * B28.
# PURPOSE: Find Dark Edge of Raster Histogram Based on the
#          DN Value at Which the Cumulative Frequency
#          Distributoin Rises Rapidly (Faster than HEP).
  func dnHE(heloc) begin
     local numeric j, jm1, jp1, sumfq, sumcfq, del, dnhe;
     sumfq = 0;
     for j=1 to maxi begin
        sumfq = sumfq + v[j];
     end
     v[1] = v[1] * 100 / sumfq;
#    v[1] equals percent frequency for DN = 1
     for j=2 to maxi begin
        jm1 = j - 1; v[j] = v[jm1] + v[j] * 100 / sumfq;
#       v[j] equals percent frequency for DN = j
     end
     dnhe = 0;
     for j=1 to maxi begin
        jp1 = j + 1; del = v[jp1] - v[j];
        if (dnhe == 0) then begin
           if (del >  heloc) then dnhe = jp1;
        end
     end
     return(dnhe);
  end
# ------------------------------------------------------------
# IF atcor = 3, RESCALE C-FACTORS BASED ON msfac:
# Refer to B23.
  if (atcor == 3) then begin
    if (pCB) then cCB = cCB * msfac;
    if (pBL) then cBL = cBL * msfac;
    cGL = cGL * msfac; cRL = cRL * msfac;
    if (pYL) then cYL = cYL * msfac;
    if (pRE) then cRE = cRE * msfac;
    if (pNA) then cNA = cNA * msfac;
    if (pNB) then cNB = cNB * msfac;
    if (pMA) then cMA = cMA * msfac;
    if (pMB) then cMB = cMB * msfac;
    if (pMC) then cMC = cMC * msfac;
    if (pMD) then cMD = cMD * msfac;
    if (pME) then cME = cME * msfac;
    if (pMF) then cMF = cMF * msfac;
    if (pMG) then cMG = cMG * msfac;
  end
# ------------------------------------------------------------
# OPEN APPROPRIATE INPUT RASTERS & GET MIN & MAX VALUES: 
# Refer to B27.
  printf("\nOPEN INPUT RASTERS:\n");
  if (pCB) then begin
     printf(" CB");
     GetInputRaster(CB); nlins = NumLins(CB);
     ncols = NumCols(CB); dntype$ = RastType(CB);
     checkHisto(CB);
     dnminCB = GlobalMin(CB); dnmaxCB = GlobalMax(CB);
  end
  if (pBL) then begin
     if (pCB) then begin
        if (pBL) then begin
           printf(" BL");
           GetInputRaster(BL,nlins,ncols,dntype$);
        end
     end
     else begin
        printf(" BL");
        GetInputRaster(BL); nlins = NumLins(BL);
        ncols = NumCols(BL); dntype$ = RastType(BL);
     end
     checkHisto(BL);
     dnminBL = GlobalMin(BL); dnmaxBL = GlobalMax(BL);
  end
  printf(" GL");
  if (pCB or pBL) then begin
     GetInputRaster(GL,nlins,ncols,dntype$);
  end
  else begin
     GetInputRaster(GL);
     nlins = NumLins(GL); ncols = NumCols(GL);
     dntype$ = RastType(GL);
  end
  checkHisto(GL);
  dnminGL = GlobalMin(GL); dnmaxGL = GlobalMax(GL);
  if (pYL) then begin
     printf(" YL");
     GetInputRaster(YL,nlins,ncols,dntype$);
     checkHisto(YL);
     dnminYL = GlobalMin(YL); dnmaxYL = GlobalMax(YL);
  end
  printf(" RL");
  GetInputRaster(RL,nlins,ncols,dntype$);
  checkHisto(RL);
  dnminRL = GlobalMin(RL); dnmaxRL = GlobalMax(RL);
  if (pRE) then begin
     printf(" RE");
     GetInputRaster(RE,nlins,ncols,dntype$);
     checkHisto(RE);
     dnminRE = GlobalMin(RE); dnmaxRE = GlobalMax(RE);
  end
  if (pNA) then begin
     printf(" NA");
     GetInputRaster(NA,nlins,ncols,dntype$);
     checkHisto(NA);
     dnminNA = GlobalMin(NA); dnmaxNA = GlobalMax(NA);
  end
  if (pNB) then begin
     printf(" NB");
     GetInputRaster(NB,nlins,ncols,dntype$);
     checkHisto(NB);
     dnminNB = GlobalMin(NB); dnmaxNB = GlobalMax(NB);
  end
  if (pMA) then begin
     printf(" MA");
     GetInputRaster(MA,nlins,ncols,dntype$);
     checkHisto(MA);
     dnminMA = GlobalMin(MA); dnmaxMA = GlobalMax(MA);
  end
  if (pMB) then begin
     printf(" MB");
     GetInputRaster(MB,nlins,ncols,dntype$);
     checkHisto(MB);
     dnminMB = GlobalMin(MB); dnmaxMB = GlobalMax(MB);
  end
  if (pMC) then begin
     printf(" MC");
     GetInputRaster(MC,nlins,ncols,dntype$);
     checkHisto(MC);
     dnminMC = GlobalMin(MC); dnmaxMC = GlobalMax(MC);
  end
  if (pMD) then begin
     printf(" MD");
     GetInputRaster(MD,nlins,ncols,dntype$);
     checkHisto(MD);
     dnminMD = GlobalMin(MD); dnmaxMD = GlobalMax(MD);
  end
  if (pME) then begin
     printf(" ME");
     GetInputRaster(ME,nlins,ncols,dntype$);
     checkHisto(ME);
     dnminME = GlobalMin(ME); dnmaxME = GlobalMax(ME);
  end
  if (pMF) then begin
     printf(" MF");
     GetInputRaster(MF,nlins,ncols,dntype$);
     checkHisto(MF);
     dnminMF = GlobalMin(MF); dnmaxMF = GlobalMax(MF);
  end
  if (pMG) then begin
     printf(" MG");
     GetInputRaster(MG,nlins,ncols,dntype$);
     checkHisto(MG);
     dnminMG = GlobalMin(MG); dnmaxMG = GlobalMax(MG);
  end
  printf("\n");
# ------------------------------------------------------------
# IF atcor > 1, FIND DARK-EDGES OF DN HISTOGRAMS:
# Refer to B28.
  if (atcor > 1) then begin
     if (pCB) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(CB,v);
        dnheCB = dnHE(he);
     end
     if (pBL) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(BL,v);
        dnheBL = dnHE(he);
        if (dnheBL < dnbBL) then dnheBL = ceil(dnbBL);
     end
     ResizeArrayClear(v,maxip1); ReadHistogram(GL,v);
     dnheGL = dnHE(he);
     if (dnheGL < dnbGL) then dnheGL = ceil(dnbGL);
     if (pYL) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(YL,v);
        dnheYL = dnHE(he);
        if (dnheYL < dnbYL) then dnheYL = ceil(dnbYL);
     end
     ResizeArrayClear(v,maxip1); ReadHistogram(RL,v);
     dnheRL = dnHE(he);
     if (dnheRL < dnbRL) then dnheRL = ceil(dnbRL);
     if (pRE) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(RE,v);
        dnheRE = dnHE(he);
        if (dnheRE < dnbRE) then dnheRE = ceil(dnbRE);
     end
     if (pNA) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(NA,v);
        dnheNA = dnHE(he);
        if (dnheNA < dnbNA) then dnheNA = ceil(dnbNA);
     end
     if (pNB) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(NB,v);
        dnheNB = dnHE(he);
        if (dnheNB < dnbNB) then dnheNB = ceil(dnbNB);
     end
     if (pMA) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(MA,v);
        dnheMA = dnHE(he);
        if (dnheMA < dnbMA) then dnheMA = ceil(dnbMA);
     end
     if (pMB) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(MB,v);
        dnheMB = dnHE(he);
        if (dnheMB < dnbMB) then dnheMB = ceil(dnbMB);
     end
     if (pMC) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(MC,v);
        dnheMC = dnHE(he);
        if (dnheMC < dnbMC) then dnheMC = ceil(dnbMC);
     end
     if (pMD) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(MD,v);
        dnheMD = dnHE(he);
        if (dnheMD < dnbMD) then dnheMD = ceil(dnbMD);
     end
     if (pME) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(ME,v);
        dnheME = dnHE(he);
        if (dnheME < dnbME) then dnheME = ceil(dnbME);
     end
     if (pMF) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(MF,v);
        dnheMF = dnHE(he);
        if (dnheMF < dnbMF) then dnheMF = ceil(dnbMF);
     end
     if (pMG) then begin
        ResizeArrayClear(v,maxip1); ReadHistogram(MG,v);
        dnheMG = dnHE(he);
        if (dnheMG < dnbMG) then dnheMG = ceil(dnbMG);
     end
  end
# ------------------------------------------------------------
# CALCULATE SItoa & DN-to-SRF CONVERSION-COEFFICIENTS: aBAND
# Refer to B29.
  esd = 1 - 0.01672 * cosd(0.9856 * (doy - 4));
  sitoafac = sind(sunelevang) / esd^2;
  if (pCB) then sitoaCB = dssiCB * sitoafac;
  if (pBL) then sitoaBL = dssiBL * sitoafac;
  sitoaGL = dssiGL * sitoafac; sitoaRL = dssiRL * sitoafac;
  if (pYL) then sitoaYL = dssiYL * sitoafac;
  if (pRE) then sitoaRE = dssiRE * sitoafac;
  if (pNA) then sitoaNA = dssiNA * sitoafac;
  if (pNB) then sitoaNB = dssiNB * sitoafac;
  if (pMA) then sitoaMA = dssiMA * sitoafac;
  if (pMB) then sitoaMB = dssiMB * sitoafac;
  if (pMC) then sitoaMC = dssiMC * sitoafac;
  if (pMD) then sitoaMD = dssiMD * sitoafac;
  if (pME) then sitoaME = dssiME * sitoafac;
  if (pMF) then sitoaMF = dssiMF * sitoafac;
  if (pMG) then sitoaMG = dssiMG * sitoafac;
  pi100 = pi * 100;
  if (pCB) then aCB = pi100 * skCB / sitoaCB;
  if (pBL) then aBL = pi100 * skBL / sitoaBL;
  aGL = pi100 * skGL / sitoaGL; aRL = pi100 * skRL / sitoaRL;
  if (pYL) then aYL = pi100 * skYL / sitoaYL;
  if (pRE) then aRE = pi100 * skRE / sitoaRE;
  if (pNA) then aNA = pi100 * skNA / sitoaNA;
  if (pNB) then aNB = pi100 * skNB / sitoaNB;
  if (pMA) then aMA = pi100 * skMA / sitoaMA;
  if (pMB) then aMB = pi100 * skMB / sitoaMB;
  if (pMC) then aMC = pi100 * skMC / sitoaMC;
  if (pMD) then aMD = pi100 * skMD / sitoaMD;
  if (pME) then aME = pi100 * skME / sitoaME;
  if (pMF) then aMF = pi100 * skMF / sitoaMF;
  if (pMG) then aMG = pi100 * skMG / sitoaMG;
# ------------------------------------------------------------
# IF atcor > 1, ESTIMATE SRFpath & DNpath: Refer to B29.
  if (atcor > 1) then begin
     dnpathGL1 = dnheGL; 
     srfpathGL1 = (dnpathGL1 - dnbGL) * aGL;
     dnpathRL = dnheRL;
     srfpathRL  = (dnpathRL - dnbRL)  * aRL;
     logsrfpathRL = log(srfpathRL);
     wLenV = wLenGL; srfpathV = srfpathGL1;
     if (pBL) then begin
        wtBL = 0.55114;
        dnpathBL1 = dnheBL; 
		dnpathCB1=dnheCB;  								# this line was added to script
        srfpathBL1 = (dnpathBL1 - dnbBL) * aBL; 
		srfpathCB1=(dnpathCB1-dnbCB)*aCB;				# this line was added to script
        wLenV = wtBL * wLenBL + (1 - wtBL) * wLenGL;
        srfpathV = (srfpathBL1 + srfpathGL1) / 2;
     end
     logwLenRLdwLenV = log(wLenRL / wLenV);
     logsrfpathV = log(srfpathV);
     difoflogs = logsrfpathV - logsrfpathRL;
     p = difoflogs / logwLenRLdwLenV;
     if (pCB) then begin
        srfpathCB2 =  srfpathRL * (wLenRL / wLenCB)^p;print( srfpathCB2);
        srfpathCB = SetMin(srfpathCB1, srfpathCB2);print( srfpathCB1);print( srfpathCB);print( dnpathCB);
        dnpathCB = round(srfpathCB / aCB) + dnbCB;
     end
     if (pBL) then begin
        srfpathBL2 =  srfpathRL * (wLenRL / wLenBL)^p;
        srfpathBL = SetMin(srfpathBL1, srfpathBL2);
        dnpathBL = round(srfpathBL / aBL) + dnbBL;
     end
     srfpathGL2 =  srfpathRL * (wLenRL / wLenGL)^p;
     srfpathGL = SetMin(srfpathGL1, srfpathGL2);
     dnpathGL = round(srfpathGL / aGL) + dnbGL;
     if (pYL) then begin
        srfpathYL2 =  srfpathRL * (wLenRL / wLenYL)^p;
        dnpathYL1 = dnheYL;
        srfpathYL1 = (dnpathYL1 - dnbYL) * aYL;
        srfpathYL = SetMin(srfpathYL1, srfpathYL2);
        dnpathYL = round(srfpathYL / aYL) + dnbYL;
     end
     p2 = 0.2 * p;
     if (pRE) then begin
        srfpathRE2 =  srfpathRL * (wLenRL / wLenRE)^p2;
        dnpathRE1 = dnheRE;
        srfpathRE1 = (dnpathRE1 - dnbRE) * aRE;
        srfpathRE = SetMin(srfpathRE1, srfpathRE2);
        dnpathRE = round(srfpathRE / aRE) + dnbRE;
     end
     if (pNA) then begin
        srfpathNA2 =  srfpathRL * (wLenRL / wLenNA)^p2;
        dnpathNA1 = dnheNA;
        srfpathNA1 = (dnpathNA1 - dnbNA) * aNA;
        srfpathNA = SetMin(srfpathNA1, srfpathNA2);
        dnpathNA = round(srfpathNA / aNA) + dnbNA;
     end
     if (pNB) then begin
        srfpathNB2 =  srfpathRL * (wLenRL / wLenNB)^p2;
        dnpathNB1 = dnheNB;
        srfpathNB1 = (dnpathNB1 - dnbNB) * aNB;
        srfpathNB = SetMin(srfpathNB1, srfpathNB2);
        dnpathNB = round(srfpathNB / aNB) + dnbNB;
     end
     if (pMA) then begin
        srfpathMA2 =  srfpathRL * (wLenRL / wLenMA)^p2;
        dnpathMA1 = dnheMA;
        srfpathMA1 = (dnpathMA1 - dnbMA) * aMA;
        srfpathMA = SetMin(srfpathMA1, srfpathMA2);
        dnpathMA = round(srfpathMA / aMA) + dnbMA;
     end
     if (pMB) then begin
        srfpathMB2 =  srfpathRL * (wLenRL / wLenMB)^p2;
        dnpathMB1 = dnheMB;
        srfpathMB1 = (dnpathMB1 - dnbMB) * aMB;
        srfpathMB = SetMin(srfpathMB1, srfpathMB2);
        dnpathMB = round(srfpathMB / aMB) + dnbMB;
     end
     if (pMC) then begin
        srfpathMC2 =  srfpathRL * (wLenRL / wLenMC)^p2;
        dnpathMC1 = dnheMC;
        srfpathMC1 = (dnpathMC1 - dnbMC) * aMC;
        srfpathMC = SetMin(srfpathMC1, srfpathMC2);
        dnpathMC = round(srfpathMC / aMC) + dnbMC;
     end
     if (pMD) then begin
        srfpathMD2 =  srfpathRL * (wLenRL / wLenMD)^p2;
        dnpathMD1 = dnheMD;
        srfpathMD1 = (dnpathMD1 - dnbMD) * aMD;
        srfpathMD = SetMin(srfpathMD1, srfpathMD2);
        dnpathMD = round(srfpathMD / aMD) + dnbMD;
     end
     if (pME) then begin
        srfpathME2 =  srfpathRL * (wLenRL / wLenME)^p2;
        dnpathME1 = dnheME;
        srfpathME1 = (dnpathME1 - dnbME) * aME;
        srfpathME = SetMin(srfpathME1, srfpathME2);
        dnpathME = round(srfpathME / aME) + dnbME;
     end
     if (pMF) then begin
        srfpathMF2 =  srfpathRL * (wLenRL / wLenMF)^p2;
        dnpathMF1 = dnheMF;
        srfpathMF1 = (dnpathMF1 - dnbMF) * aMF;
        srfpathMF = SetMin(srfpathMF1, srfpathMF2);
        dnpathMF = round(srfpathMF / aMF) + dnbMF;
     end
     if (pMG) then begin
        srfpathMG2 =  srfpathRL * (wLenRL / wLenMG)^p2;
        dnpathMG1 = dnheMG;
        srfpathMG1 = (dnpathMG1 - dnbMG) * aMG;
        srfpathMG = SetMin(srfpathMG1, srfpathMG2);
        dnpathMG = round(srfpathMG / aMG) + dnbMG;
     end
  end
# ------------------------------------------------------------
# SET UP LOOK-UP TABLES FOR SRFI VALUES (FUNCTION OF DNs):
# Refer to B30.
  if (atcor == 1) then begin
     srfpathCB = 0; srfpathBL = 0; srfpathGL = 0;
     srfpathYL = 0; srfpathRL = 0; srfpathRE = 0;
     srfpathNA = 0; srfpathNB = 0; srfpathMA = 0;
     srfpathMB = 0; srfpathMC = 0; srfpathMD = 0;
     srfpathME = 0; srfpathMF = 0; srfpathMG = 0;
  end
  if (cCB) then ResizeArrayClear(vCB,maxip1);
  if (pBL) then ResizeArrayClear(vBL,maxip1);
  ResizeArrayClear(vGL,maxip1); ResizeArrayClear(vRL,maxip1);
  if (pYL) then ResizeArrayClear(vYL,maxip1);
  if (pRE) then ResizeArrayClear(vRE,maxip1);
  if (pNA) then ResizeArrayClear(vNA,maxip1);
  if (pNB) then ResizeArrayClear(vNB,maxip1);
  if (pMA) then ResizeArrayClear(vMA,maxip1);
  if (pMB) then ResizeArrayClear(vMB,maxip1);
  if (pMC) then ResizeArrayClear(vMC,maxip1);
  if (pMD) then ResizeArrayClear(vMD,maxip1);
  if (pME) then ResizeArrayClear(vME,maxip1);
  if (pMF) then ResizeArrayClear(vMF,maxip1);
  if (pMG) then ResizeArrayClear(vMG,maxip1);
  for i=1 to maxi begin
     if (pCB) then begin
        srftoaCB = (i - dnbCB) * aCB;
        srfapcCB = srftoaCB - srfpathCB;
        srfsfcCB = srfapcCB * cCB;
        srfiCB = round(srfsfcCB * 100);
        if (srfiCB < 1) then srfiCB = 1;
        vCB[i] = srfiCB;
     end
     if (pBL) then begin
        srftoaBL = (i - dnbBL) * aBL;
        srfapcBL = srftoaBL - srfpathBL;
        srfsfcBL = srfapcBL * cBL;
        srfiBL = round(srfsfcBL * 100);
        if (srfiBL < 1) then srfiBL = 1;
        vBL[i] = srfiBL;
     end
     srftoaGL = (i - dnbGL) * aGL; 
     srfapcGL = srftoaGL - srfpathGL;
     srfsfcGL = srfapcGL * cGL;
     srfiGL = round(srfsfcGL * 100);
     if (srfiGL < 1) then srfiGL = 1;
     vGL[i] = srfiGL;
     if (pYL) then begin
        srftoaYL = (i - dnbYL) * aYL;
        srfapcYL = srftoaYL - srfpathYL;
        srfsfcYL = srfapcYL * cYL;
        srfiYL = round(srfsfcYL * 100);
        if (srfiYL < 1) then srfiYL = 1;
        vYL[i] = srfiYL;
     end
     srftoaRL = (i - dnbRL) * aRL;
     srfapcRL = srftoaRL - srfpathRL;
     srfsfcRL = srfapcRL * cRL;
     srfiRL = round(srfsfcRL * 100);
     if (srfiRL < 1) then srfiRL = 1;
     vRL[i] = srfiRL;
     if (pRE) then begin
        srftoaRE = (i - dnbRE) * aRE;
        srfapcRE = srftoaRE - srfpathRE;
        srfsfcRE = srfapcRE * cRE;
        srfiRE = round(srfsfcRE * 100);
        if (srfiRE < 1) then srfiRE = 1;
        vRE[i] = srfiRE;
     end
     if (pNA) then begin
        srftoaNA = (i - dnbNA) * aNA;
        srfapcNA = srftoaNA - srfpathNA;
        srfsfcNA = srfapcNA * cNA;
        srfiNA = round(srfsfcNA * 100);
        if (srfiNA < 1) then srfiNA = 1;
        vNA[i] = srfiNA;
     end
     if (pNB) then begin
        srftoaNB = (i - dnbNB) * aNB;
        srfapcNB = srftoaNB - srfpathNB;
        srfsfcNB = srfapcNB * cNB;
        srfiNB = round(srfsfcNB * 100);
        if (srfiNB < 1) then srfiNB = 1;
        vNB[i] = srfiNB;
     end
     if (pMA) then begin
        srftoaMA = (i - dnbMA) * aMA;
        srfapcMA = srftoaMA - srfpathMA;
        srfsfcMA = srfapcMA * cMA;
        srfiMA = round(srfsfcMA * 100);
        if (srfiMA < 1) then srfiMA = 1;
        vMA[i] = srfiMA;
     end
     if (pMB) then begin
        srftoaMB = (i - dnbMB) * aMB;
        srfapcMB = srftoaMB - srfpathMB;
        srfsfcMB = srfapcMB * cMB;
        srfiMB = round(srfsfcMB * 100);
        if (srfiMB < 1) then srfiMB = 1;
        vMB[i] = srfiMB;
     end
     if (pMC) then begin
        srftoaMC = (i - dnbMC) * aMC;
        srfapcMC = srftoaMC - srfpathMC;
        srfsfcMC = srfapcMC * cMC;
        srfiMC = round(srfsfcMC * 100);
        if (srfiMC < 1) then srfiMC = 1;
        vMC[i] = srfiMC;
     end
     if (pMD) then begin
        srftoaMD = (i - dnbMD) * aMD;
        srfapcMD = srftoaMD - srfpathMD;
        srfsfcMD = srfapcMD * cMD;
        srfiMD = round(srfsfcMD * 100);
        if (srfiMD < 1) then srfiMD = 1;
        vMD[i] = srfiMD;
     end
     if (pME) then begin
        srftoaME = (i - dnbME) * aME;
        srfapcME = srftoaME - srfpathME;
        srfsfcME = srfapcME * cME;
        srfiME = round(srfsfcME * 100);
        if (srfiME < 1) then srfiME = 1;
        vME[i] = srfiME;
     end
     if (pMF) then begin
        srftoaMF = (i - dnbMF) * aMF;
        srfapcMF = srftoaMF - srfpathMF;
        srfsfcMF = srfapcMF * cMF;
        srfiMF = round(srfsfcMF * 100);
        if (srfiMF < 1) then srfiMF = 1;
        vMF[i] = srfiMF;
     end
     if (pMG) then begin
        srftoaMG = (i - dnbMG) * aMG;
        srfapcMG = srftoaMG - srfpathMG;
        srfsfcMG = srfapcMG * cMG;
        srfiMG = round(srfsfcMG * 100);
        if (srfiMG < 1) then srfiMG = 1;
        vMG[i] = srfiMG;
     end
  end
# ------------------------------------------------------------
# IF atcor > 1, THEN GENERATE A HISTOGRAM ANALYSES REPORT:
# Refer to B31.
  if (atcor > 1) then begin
     printf("\nHISTOGRAM ANALYSES REPORT:\n");
     printf("BAND DNmin DNhe DNmax");
     printf(" DNpath SRFpath1 SRFpath\n");
     if (pCB) then begin
        printf("  CB   %3d  %3d",dnminCB,dnheCB);
        printf("  %4d    %3d",dnmaxCB, dnpathCB);
        printf(" %8.3f %7.3f\n",srfpathCB1,srfpathCB);
     end
     if (pBL) then begin
        printf("  BL   %3d  %3d",dnminBL,dnheBL);
        printf("  %4d    %3d",dnmaxBL, dnpathBL);
        printf(" %8.3f %7.3f\n",srfpathBL1,srfpathBL);
     end
     printf("  GL   %3d  %3d",dnminGL,dnheGL);
     printf("  %4d    %3d",dnmaxGL, dnpathGL);
     printf(" %8.3f %7.3f\n",srfpathGL1,srfpathGL);
     if (pYL) then begin
        printf("  YL   %3d  %3d",dnminYL,dnheYL);
        printf("  %4d    %3d",dnmaxYL, dnpathYL);
        printf(" %8.3f %7.3f\n",srfpathYL1,srfpathYL);
     end
     printf("  RL   %3d  %3d",dnminRL,dnheRL);
     printf("  %4d    %3d",dnmaxRL, dnpathRL);
     printf(" %8.3f %7.3f\n",srfpathRL,srfpathRL);
     if (pRE) then begin
        printf("  RE   %3d  %3d",dnminRE,dnheRE);
        printf("  %4d    %3d",dnmaxRE, dnpathRE);
        printf(" %8.3f %7.3f\n",srfpathRE1,srfpathRE);
     end
     if (pNA) then begin
        printf("  NA   %3d  %3d",dnminNA,dnheNA);
        printf("  %4d    %3d",dnmaxNA, dnpathNA);
        printf(" %8.3f %7.3f\n",srfpathNA1,srfpathNA);
     end
     if (pNB) then begin
        printf("  NB   %3d  %3d",dnminNB,dnheNB);
        printf("  %4d    %3d",dnmaxNB, dnpathNB);
        printf(" %8.3f %7.3f\n",srfpathNB1,srfpathNB);
     end
     if (pMA) then begin
        printf("  MA   %3d  %3d",dnminMA,dnheMA);
        printf("  %4d    %3d",dnmaxMA, dnpathMA);
        printf(" %8.3f %7.3f\n",srfpathMA1,srfpathMA);
     end
     if (pMB) then begin
        printf("  MB   %3d  %3d",dnminMB,dnheMB);
        printf("  %4d    %3d",dnmaxMB, dnpathMB);
        printf(" %8.3f %7.3f\n",srfpathMB1,srfpathMB);
     end
     if (pMC) then begin
        printf("  MC   %3d  %3d",dnminMC,dnheMC);
        printf("  %4d    %3d",dnmaxMC, dnpathMC);
        printf(" %8.3f %7.3f\n",srfpathMC1,srfpathMC);
     end
     if (pMD) then begin
        printf("  MD   %3d  %3d",dnminMD,dnheMD);
        printf("  %4d    %3d",dnmaxMD, dnpathMD);
        printf(" %8.3f %7.3f\n",srfpathMD1,srfpathMD);
     end
     if (pME) then begin
        printf("  ME   %3d  %3d",dnminME,dnheME);
        printf("  %4d    %3d",dnmaxME, dnpathME);
        printf(" %8.3f %7.3f\n",srfpathME1,srfpathME);
     end
     if (pMF) then begin
        printf("  MF   %3d  %3d",dnminMF,dnheMF);
        printf("  %4d    %3d",dnmaxMF, dnpathMF);
        printf(" %8.3f %7.3f\n",srfpathMF1,srfpathMF);
     end
     if (pMG) then begin
        printf("  MG   %3d  %3d",dnminMG,dnheMG);
        printf("  %4d    %3d",dnmaxMG, dnpathMG);
        printf(" %8.3f %7.3f\n",srfpathMG1,srfpathMG);
     end
     printf("SRFIpath Model Parameter: %6.4f\n",p);
#    ---------------------------------------------------------
#    PRINT HISTOGRAM ANALYSIS REPORT TO CONSOLE WINDOW:
     printf("\n");
     p1$   = "HISTOGRAM ANALYSES REPORT ASSESSMENT\n";
     p2$   = "* EXAMINE the HISTOGRAM ANALYSES REPORT.\n";
     p3$   = "  - DNheXX Should be Somewhat Higher";
     p4$   = " than the DNminXX.\n";
     p5$   = "  - SRFpath1XX Should be Like SRFpathXX\n";
     p6$   = "       EXCEPT for NA and/or NB Band.\n";
     p7$   = "  - SRFpath Should Decrease Band by Band.\n";
     p8$   = "  - SRFIpath Model Parameter Value Should\n";
     p9$   = "       be in the RANGE: 1.5000 to 5.0000.\n\n";
     p10$  = "* If the REPORT Looks Reasonable, Continue.\n";
     p11$  = "* If the REPORT Looks Unreasonable, STOP.\n\n";
     p12$  = "  - TO STOP (in the NEXT Popup WINDOW),\n";
     p13$  = "       * CLICK CANCEL.\n";
     p14$  = "       * Then, Re-Run SRFI.sml with a\n";
     p15$  = "                   Better delcf Value.\n\n";
     p16$  = "   - To CLOSE This Window, CLICK OK.";
     p$    = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
     p$    = p$ + p9$ + p10$ + p11$ + p12$ + p13$ + p14$;
     p$    = p$ + p15$ + p16$;
     PopupMessage(p$);
  end
# ------------------------------------------------------------
# ASSIGN NAMES & LOCATIONS FOR OUTPUT RASTERS:
# Refer to B33.
  printf("ASSIGN NAMES & LOCATIONS FOR OUTPUT RASTERS:\n");
  srfitype$ = "16-bit unsigned"; ptype$ = "16-bit unsigned";
  if (pCB) then begin
     printf(" SRFICB");
     GetOutputRaster(SRFICB,nlins,ncols,srfitype$);
     SetNull(SRFICB,0);
  end
  if (pBL) then begin
     printf(" SRFIBL");
     GetOutputRaster(SRFIBL,nlins,ncols,srfitype$);
     SetNull(SRFIBL,0);
  end
  printf(" SRFIGL");
  GetOutputRaster(SRFIGL,nlins,ncols,srfitype$);
  SetNull(SRFIGL,0);
  if (pYL) then begin
     printf(" SRFIYL");
     GetOutputRaster(SRFIYL,nlins,ncols,srfitype$);
     SetNull(SRFIYL,0);
  end
  printf(" SRFIRL");
  GetOutputRaster(SRFIRL,nlins,ncols,srfitype$);
  SetNull(SRFIRL,0);
  if (pRE) then begin
     printf(" SRFIRE");
     GetOutputRaster(SRFIRE,nlins,ncols,srfitype$);
     SetNull(SRFIRE,0);
  end
  if (pNA) then begin
     printf(" SRFINA");
     GetOutputRaster(SRFINA,nlins,ncols,srfitype$);
     SetNull(SRFINA,0);
  end
  if (pNB) then begin
     printf(" SRFINB");
     GetOutputRaster(SRFINB,nlins,ncols,srfitype$);
     SetNull(SRFINB,0);
  end
  if (pMA) then begin
     printf(" SRFIMA");
     GetOutputRaster(SRFIMA,nlins,ncols,srfitype$);
     SetNull(SRFIMA,0);
  end
  if (pMB) then begin
     printf(" SRFIMB");
     GetOutputRaster(SRFIMB,nlins,ncols,srfitype$);
     SetNull(SRFIMB,0);
  end
  if (pMC) then begin
     printf(" SRFIMC");
     GetOutputRaster(SRFIMC,nlins,ncols,srfitype$);
     SetNull(SRFIMC,0);
  end
  if (pMD) then begin
     printf(" SRFIMD");
     GetOutputRaster(SRFIMD,nlins,ncols,srfitype$);
     SetNull(SRFIMD,0);
  end
  if (pME) then begin
     printf(" SRFIME");
     GetOutputRaster(SRFIME,nlins,ncols,srfitype$);
     SetNull(SRFIME,0);
  end
  if (pMF) then begin
     printf(" SRFIMF");
     GetOutputRaster(SRFIMF,nlins,ncols,srfitype$);
     SetNull(SRFIMF,0);
  end
  if (pMG) then begin
     printf(" SRFIMG");
     GetOutputRaster(SRFIMG,nlins,ncols,srfitype$);
     SetNull(SRFIMG,0);
  end
  if (atcor == 3) then begin
     printf(" PVI");
     GetOutputRaster(PVI,nlins,ncols,ptype$);
     printf(" PBI");
     GetOutputRaster(PBI,nlins,ncols,ptype$);
     SetNull(PVI,0); SetNull(PBI,0);
  end
  printf("\n\n");
# ------------------------------------------------------------
# FILL SRFI RASTERS WITH VALUES: Refer to B35.
  if (atcor == 1) then begin
     printf("PRODUCING SRFItoa RASTERS:\n");
  end
  if (atcor == 2) then begin
     printf("PRODUCING SRFIapc RASTERS:\n");
  end
  if (atcor == 3) then begin
     printf("PRODUCING SRFIsfc RASTERS:\n");
  end
  if (pCB) then begin
     printf(" CB");
     for each CB begin
        dnCB = CB; nc = IsNull(dnCB);
        if (nc==0) then SRFICB = vCB[dnCB];
     end
     CreateHistogram(SRFICB,0); CreatePyramid(SRFICB,0);
     CopySubobjects(RL,SRFICB,"GEOREF");
     smlContrast.InputUpperLimit = 2820;
     smlContrast.Compute(SRFICB,clut1$,clut2$,clut3$);
     CloseRaster(CB); CloseRaster(SRFICB);
  end
  if (pBL) then begin
     printf(" BL");
     for each BL begin
        dnBL = BL; nc = IsNull(dnBL);
        if (nc==0) then SRFIBL = vBL[dnBL];
     end
     CreateHistogram(SRFIBL,0); CreatePyramid(SRFIBL,0);
     CopySubobjects(RL,SRFIBL,"GEOREF");
     smlContrast.InputUpperLimit = 2820;
     smlContrast.Compute(SRFIBL,clut1$,clut2$,clut3$);
     CloseRaster(BL); CloseRaster(SRFIBL);
  end
  printf(" GL");
  for each GL begin
     dnGL = GL; nc = IsNull(dnGL);
     if (nc==0) then SRFIGL = vGL[dnGL];
  end
  CreateHistogram(SRFIGL,0); CreatePyramid(SRFIGL,0);
  CopySubobjects(RL,SRFIGL,"GEOREF");
  smlContrast.InputUpperLimit = 3230;
  smlContrast.Compute(SRFIGL,clut1$,clut2$,clut3$);
  CloseRaster(GL); CloseRaster(SRFIGL);
  if (pYL) then begin
     printf(" YL");
     for each YL begin
        dnYL = YL; nc = IsNull(dnYL);
        if (nc==0) then SRFIYL = vYL[dnYL];
     end
     CreateHistogram(SRFIYL,0); CreatePyramid(SRFIYL,0);
     CopySubobjects(RL,SRFIYL,"GEOREF");
     smlContrast.InputUpperLimit = 3230;
     smlContrast.Compute(SRFIYL,clut1$,clut2$,clut3$);
     CloseRaster(YL); CloseRaster(SRFIYL);
  end
  printf(" RL");
  for each RL begin
     dnRL = RL; nc = IsNull(dnRL);
     if (nc==0) then SRFIRL = vRL[dnRL];
  end
  CreateHistogram(SRFIRL,0); CreatePyramid(SRFIRL,0);
  CopySubobjects(RL,SRFIRL,"GEOREF");
  smlContrast.InputUpperLimit = 3500;
  smlContrast.Compute(SRFIRL,clut1$,clut2$,clut3$);
  CloseRaster(RL);
  if (pRE) then begin
     printf(" RE");
     for each RE begin
        dnRE = RE; nc = IsNull(dnRE);
        if (nc==0) then SRFIRE = vRE[dnRE];
     end
     CreateHistogram(SRFIRE,0); CreatePyramid(SRFIRE,0);
     CopySubobjects(RL,SRFIRE,"GEOREF");
     smlContrast.InputUpperLimit = 4000;
     smlContrast.Compute(SRFIRE,clut1$,clut2$,clut3$);
     CloseRaster(RE); CloseRaster(SRFIRE);
  end
  if (pNA) then begin
     printf(" NA");
     for each NA begin
        dnNA = NA; nc = IsNull(dnNA);
        if (nc==0) then SRFINA = vNA[dnNA];
     end
     CreateHistogram(SRFINA,0); CreatePyramid(SRFINA,0);
     CopySubobjects(RL,SRFINA,"GEOREF");
     smlContrast.InputUpperLimit = 5000;
     smlContrast.Compute(SRFINA,clut1$,clut2$,clut3$);
     CloseRaster(NA);
  end
  if (pNB) then begin
     printf(" NB");
     for each NB begin
        dnNB = NB; nc = IsNull(dnNB); 					# here was an error: oryginalnie dnNA=NB
        if (nc==0) then SRFINB = vNB[dnNB];
     end
     CreateHistogram(SRFINB,0); CreatePyramid(SRFINB,0);
     CopySubobjects(RL,SRFINB,"GEOREF");
     smlContrast.InputUpperLimit = 5000;
     smlContrast.Compute(SRFINB,clut1$,clut2$,clut3$);
     CloseRaster(NB);
  end
  if (pMA) then begin
     printf(" MA");
     for each MA begin
        dnMA = MA; nc = IsNull(dnMA);
        if (nc==0) then SRFIMA = vMA[dnMA];
     end
     CreateHistogram(SRFIMA,0); CreatePyramid(SRFIMA,0);
     CopySubobjects(RL,SRFIMA,"GEOREF");
     smlContrast.InputUpperLimit = 4200;
     smlContrast.Compute(SRFIMA,clut1$,clut2$,clut3$);
     CloseRaster(MA); CloseRaster(SRFIMA);
  end
  if (pMB) then begin
     printf(" MB");
     for each MB begin
        dnMB = MB; nc = IsNull(dnMB);
        if (nc==0) then SRFIMB = vMB[dnMB];
     end
     CreateHistogram(SRFIMB,0); CreatePyramid(SRFIMB,0);
     CopySubobjects(RL,SRFIMB,"GEOREF");
     smlContrast.InputUpperLimit = 4200;
     smlContrast.Compute(SRFIMB,clut1$,clut2$,clut3$);
     CloseRaster(MB); CloseRaster(SRFIMB);
  end
  if (pMC) then begin
     printf(" MC");
     for each MC begin
        dnMC = MC; nc = IsNull(dnMC);
        if (nc==0) then SRFIMC = vMC[dnMC];
     end
     CreateHistogram(SRFIMC,0); CreatePyramid(SRFIMC,0);
     CopySubobjects(RL,SRFIMC,"GEOREF");
     smlContrast.InputUpperLimit = 3500;
     smlContrast.Compute(SRFIMC,clut1$,clut2$,clut3$);
     CloseRaster(MC); CloseRaster(SRFIMC);
  end
  if (pMD) then begin
     printf(" MD");
     for each MD begin
        dnMD = MD; nc = IsNull(dnMD);
        if (nc==0) then SRFIMD = vMD[dnMD];
     end
     CreateHistogram(SRFIMD,0); CreatePyramid(SRFIMD,0);
     CopySubobjects(RL,SRFIMD,"GEOREF");
     smlContrast.InputUpperLimit = 3500;
     smlContrast.Compute(SRFIMD,clut1$,clut2$,clut3$);
     CloseRaster(MD); CloseRaster(SRFIMD);
  end
  if (pME) then begin
     printf(" ME");
     for each ME begin
        dnME = ME; nc = IsNull(dnME);
        if (nc==0) then SRFIME = vME[dnME];
     end
     CreateHistogram(SRFIME,0); CreatePyramid(SRFIME,0);
     CopySubobjects(RL,SRFIME,"GEOREF");
     smlContrast.InputUpperLimit = 3500;
     smlContrast.Compute(SRFIME,clut1$,clut2$,clut3$);
     CloseRaster(ME); CloseRaster(SRFIME);
  end
  if (pMF) then begin
     printf(" MF");
     for each MF begin
        dnMF = MF; nc = IsNull(dnMF);
        if (nc==0) then SRFIMF = vMF[dnMF];
     end
     CreateHistogram(SRFIMF,0); CreatePyramid(SRFIMF,0);
     CopySubobjects(RL,SRFIMF,"GEOREF");
     smlContrast.InputUpperLimit = 3500;
     smlContrast.Compute(SRFIMF,clut1$,clut2$,clut3$); 
     CloseRaster(MF); CloseRaster(SRFIMF);
  end
  if (pMG) then begin
     printf(" MG");
     for each MG begin
        dnMG = MG; nc = IsNull(dnMG);
        if (nc==0) then SRFIMG = vMG[dnMG];
     end
     CreateHistogram(SRFIMG,0); CreatePyramid(SRFIMG,0);
     CopySubobjects(RL,SRFIMG,"GEOREF");
     smlContrast.InputUpperLimit = 3500;
     smlContrast.Compute(SRFIMG,clut1$,clut2$,clut3$);
     CloseRaster(MG); CloseRaster(SRFIMG);
  end
# ------------------------------------------------------------
# IF atcor = 3, PRODUCE VALUES FOR PVI & PBI RASTERS:
# Refer to B36
  if (atcor == 3) then begin
     printf("\n\nPRODUCING PVI & PBI RASTERS:");
     angrad = -atan(bslineslope);
     sinang = sin(angrad); cosang = cos(angrad);
     for each SRFIRL begin
        srfiRL =  SRFIRL;
        if (srfiRL > 0) then begin
           if (pNA) then srfiNR = SRFINA;
           if (pNB and pNA==0) then srfiNR = SRFINB;
           tsrfiNR = srfiNR - srfiNRint;
           pbif = srfiRL * cosang - tsrfiNR * sinang;
           pbi  = round(pfac * pbif)+ pbioff;
           pvif = srfiRL * sinang + tsrfiNR * cosang;
           pvi  = round(pfac * pvif) + pvioff;
           if (pbi < 1) then pbi = 1;
           if (pvi < 1) then pvi = 1;
           if (pbi > maxpvipbi) then pbi = maxpvipbi;
           if (pvi > maxpvipbi) then pvi = maxpvipbi;
           PBI = pbi; PVI = pvi;
        end
     end
     CreateHistogram(PBI,0); CreateHistogram(PVI,0);
     CreatePyramid(PBI,0); CreatePyramid(PVI,0);
     CopySubobjects(RL,PVI,"GEOREF");
     CopySubobjects(RL,PBI,"GEOREF");
     smlContrast.Power           =    1;
     smlContrast.InputLowerLimit =  800;
     smlContrast.InputUpperLimit = 2000;
     smlContrast.Compute(PVI,clut1$,clut2$,clut3$);
     smlContrast.InputLowerLimit =    1;
     smlContrast.InputUpperLimit = 1500;
     smlContrast.Compute(PBI,clut1$,clut2$,clut3$);
     CloseRaster(PBI); CloseRaster(PVI);
  end
  if (pNA) then CloseRaster(SRFINA);
  if (pNB) then CloseRaster(SRFINB);
  CloseRaster(SRFIRL);
# ------------------------------------------------------------
# CLEAR CONSOLE WINDOW, WRITE TITLE & WRITE PROCESSING-REPORT
# TO CONSOLE WINDOW: Refer to B37.
  clear(); writeTitle();
  printf("SRFI.sml PROCESSING-REPORT:\n\n");
  printf("          SITE NAME: %s\n",site$);
  printf("    COLLECTION DATE: %8d\n",cdate);
  printf("        DAY OF YEAR: %d\n",doy);
  if (pdate>0) then begin
     printf("    PROCESSING DATE: %8d\n",pdate);
  end
  printf("SUN ELEVATION ANGLE: %5.2f deg\n",sunelevang);
  printf("             IMAGER: %s\n",imager$);
  printf("     IMAGE DN SCALE: %d to %d\n",dnlow,dnhigh);
  if (imager == 4) then begin
     printf("          GAIN-CODE: %s\n",gc$);
     p$ = "LPGS (EOS Data Gateway)";
     if (sc == 2) then p$ = "NLAPS (EarthExplorer)";
     printf("        DATA VENDOR: %s\n",p$);
  end
  printf("   CORRECTION LEVEL: ");
  p4$ = "Make No Atmospheric Corrections";
  p5$ = "Correct for Atmospheric Reflectance Only";
  p6$ = "Correct for All Atmospheric Effects";
  if (atcor == 1) then printf("%s\n",p4$);
  if (atcor == 2) then printf("%s\n",p5$);
  if (atcor == 3) then printf("%s\n",p6$);
  printf("              delcf: %4.2f Percentage",delcf);
  printf(" Points\n");
  printf("              msfac: %6.4f\n",msfac);
  printf("               icRL: %6.4f\n\n",icRL);
  if (imager == 4) then begin
     printf("          GAIN-CODE: %s\n",gc$);
     p$ = "LPGS (EOS Data Gateway)";
     if (sc == 2) then p$ = "NLAPS (EarthExplorer)";
     printf("        DATA VENDOR: %s\n",p$);
  end
  printf("   BAND  dnb dnpath srfpath        a");
  printf("     c         m\n");
  if (pCB) then begin
     mCB = aCB * cCB;
     printf("     CB %4d %6d %7.3f",dnbCB,dnpathCB,srfpathCB);
     printf(" %8.6f %5.3f %9.7f\n",aCB,cCB,mCB);
  end
  if (pBL) then begin
     mBL = aBL * cBL;
     printf("     BL %4d %6d %7.3f",dnbBL,dnpathBL,srfpathBL);
     printf(" %8.6f %5.3f %9.7f\n",aBL,cBL,mBL);
  end
  mGL = aGL * cGL;
  if (pYL) then begin
     mYL = aYL * cYL;
     printf("     YL %4d %6d %7.3f",dnbYL,dnpathYL,srfpathYL);
     printf(" %8.6f %5.3f %9.7f\n",aYL,cYL,mYL);
  end
  mRL = aRL * cRL;
  printf("     GL %4d %6d %7.3f",dnbGL,dnpathGL,srfpathGL);
  printf(" %8.6f %5.3f %9.7f\n",aGL,cGL,mGL);
  printf("     RL %4d %6d %7.3f",dnbRL,dnpathRL,srfpathRL);
  printf(" %8.6f %5.3f %9.7f\n",aRL,cRL,mRL);
  if (pRE) then begin
     mRE = aRE * cRE;
     printf("     RE %4d %6d %7.3f",dnbRE,dnpathRE,srfpathRE);
     printf(" %8.6f %5.3f %9.7f\n",aRE,cRE,mRE);
  end
  if (pNA) then begin
     mNA = aNA * cNA;
     printf("     NA %4d %6d %7.3f",dnbNA,dnpathNA,srfpathNA);
     printf(" %8.6f %5.3f %9.7f\n",aNA,cNA,mNA);
  end
  if (pNB) then begin
     mNB = aNB * cNB;
     printf("     NB %4d %6d %7.3f",dnbNB,dnpathNB,srfpathNB);
     printf(" %8.6f %5.3f %9.7f\n",aNB,cNB,mNB);
  end
  if (pMA) then begin
     mMA = aMA * cMA;
     printf("     MA %4d %6d %7.3f",dnbMA,dnpathMA,srfpathMA);
     printf(" %8.6f %5.3f %9.7f\n",aMA,cMA,mMA);
  end
  if (pMB) then begin
     mMB = aMB * cMB;
     printf("     MB %4d %6d %7.3f",dnbMB,dnpathMB,srfpathMB);
     printf(" %8.6f %5.3f %9.7f\n",aMB,cMB,mMB);
  end
  if (pMC) then begin
     mMC = aMC * cMC;
     printf("     MC %4d %6d %7.3f",dnbMC,dnpathMC,srfpathMC);
     printf(" %8.6f %5.3f %9.7f\n",aMC,cMC,mMC);
  end
  if (pMD) then begin
     mMD = aMD * cMD;
     printf("     MD %4d %6d %7.3f",dnbMD,dnpathMD,srfpathMD);
     printf(" %8.6f %5.3f %9.7f\n",aMD,cMD,mMD);
  end
  if (pME) then begin
     mME = aME * cME;
     printf("     ME %4d %6d %7.3f",dnbME,dnpathME,srfpathME);
     printf(" %8.6f %5.3f %9.7f\n",aME,cME,mME);
  end
  if (pMF) then begin
     mMF = aMF * cMF;
     printf("     MF %4d %6d %7.3f",dnbMF,dnpathMF,srfpathMF);
     printf(" %8.6f %5.3f %9.7f\n",aMF,cMF,mMF);
  end
  if (pMG) then begin
     mMG = aMG * cMG;
     printf("     MG %4d %6d %7.3f",dnbMG,dnpathMG,srfpathMG);
     printf(" %8.6f %5.3f %9.7f\n",aMG,cMG,mMG);
  end
  printf("     SRFpath Model Parameter: %6.4f\n\n",p);
  printf("   CONVERSION EQUATIONS for Band XX:\n");
  printf("     mXX       = aXX * cXX\n");
  printf("     srftoaXX  = (dnXX - dnbXX) * aXX\n");
  printf("     srfpathXX = (dnpathXX - dnbXX) * aXX\n");
  printf("     srfapcXX  = srftoaXX - srfpathXX\n");
  printf("     srfsfcXX  = srfapcXX * cXX\n");
  printf("     srfsfcXX  = (dnXX - dnbXX - dnpathXX)");
  printf(" * mXX\n");
  printf("     srfitoaXX = round(100 * srftoaXX)\n");
  printf("     srfiapcXX = round(100 * srfapcXX)\n");
  printf("     srfisfcXX = round(100 * srfsfcXX)\n");
  printf("     dnXX      = dnpathXX + round");
  printf("(0.01 * srfiXX / mXX)\n");
# ------------------------------------------------------------
# GIVE REFERENCE FOR ANALYSIS HINTS: Refer to B37.
  printf("\nFOR ANALYSIS HINTS: Refer to B37.\n\n");
  printf("TO SAVE THE CONSOLE WINDOW TEXT AS A REPORT:\n");
  printf("   1. RIGHT-CLICK IN THE CONSOLE WINDOW.\n");
  printf("   2. SELECT THE Save As... OPTION.\n");
  printf("   3. NAVIGATE TO THE DESIRED LOCATION.\n");
  printf("   4. PROVIDE A REPORT NAME (or OVERWRITE).\n");
  printf("   5. CLICK OK.");