# ------------------------------------------------------------
# TERCOR.sml
# ------------------------------------------------------------
# SET WARNING LEVEL: Refer to B5 & B6.
  $warnings 3
# ------------------------------------------------------------
# DEFINE VARIABLES RELATED TO TERCOR MODEL:
  numeric p,flr;
  numeric s,fCB,fBL,fGL,fYL,fRL,fRE;
  numeric fNA,fNB,fMA,fMB,fMC,fMD,fME,fMF,fMG;
  numeric srfi1,srfi,srfimax;
  flr = 0.12;
  p = 2.0;
  srfimax = 10000;
# ------------------------------------------------------------
# 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 & REQUEST REPOSITIONING:
# Refer to B7, B8 & B9.
  clear();
  p1$ = "CONSOLE-WINDOW ADJUSTMENT\n";
  p2$ = "* REPOSITION the CONSOLE WINDOW.\n";
  p3$ = "* Then, CLICK the OK Button.";
  p$  = p1$ + p2$ + p3$;
  PopupMessage(p$);
# ------------------------------------------------------------
# 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
# ------------------------------------------------------------
# DEFINE PROCEDURE writeTitle: Refer to B13 & B14.
# PURPOSE: WRITES TITLE & AUTHOR INFO TO CONSOLE WINDOW.
  proc writeTitle() begin
     printf("TERCOR.sml:\n");
     printf("         VERSION: November 15, 2005\n");
     printf("         PURPOSE: CORRECT SRFI RASTERS ");
     printf("FOR TERRAIN-SHADING EFFECTS\n");
     printf("         DETAILS: FAQs_by_Jack A & D\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
# ------------------------------------------------------------
# WRITE TITLE & AUTHOR INFORMATION: Refer to B14.
  writeTitle();
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO USER INPUTS:
# Refer to B15 & B16.
  string site$ = "Stockton, CA";
  numeric imager,vnir;
# ------------------------------------------------------------
# DECLARE VARIABLES DIRECTLY RELATED TO USER PROVIDED VALUES:
# Refer to B15 & B16.
  string imager$;
# ------------------------------------------------------------
# DECLARE STRING VARIABLE RELATED TO SRFI RASTERS:
# Refer to DXX. 
  string rtype$ = "16-bit unsigned";
  string stype$ =  "8-bit unsigned";
# ------------------------------------------------------------
# DECLARE STRING VARIABLE RELATED TO ALL RASTER TYPES:
# Refer to DXX. 
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO BOOLEAN ENABLERS:
# Refer to B15 & B16
  numeric pCB,pBL,pGL,pYL,pRL,pRE;
  numeric pNA,pNB,pMA,pMB,pMC,pMD,pME,pMF,pMG,iip,isp,ipp;
# ------------------------------------------------------------
# DECLARE VARIBLES RELATED TO FOR-EACH LOOPS:
  numeric i,nlins,ncols,maxi;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO PVI & PBI: Refer to B15 & B16.
  numeric pfac,angrad,sinang,cosang;
  numeric pbioff,pvioff,srfiNRint,bslineslope,maxpvipbi;
  numeric srfiRL,srfiNR,tsrfiNR,pbif,pbi,pvif,pvi;
# ------------------------------------------------------------
# DECLARE LIST OF POSSIBLE INPUT RASTERS:
  raster SRFICB1,SRFIBL1,SRFIGL1,SRFIYL1,SRFIRL1,SRFIRE1;
  raster SRFINA1,SRFINB1,SRFIMA1,SRFIMB1,SRFIMC1,SRFIMD1;
  raster SRFIME1,SRFIMF1,SRFIMG1;
  raster PBI1,PVI1,SHADING;
# ------------------------------------------------------------
# DECLARE LIST OF POSSIBLE OUTPUT RASTERS:
  raster SRFICB,SRFIBL,SRFIGL,SRFIYL,SRFIRL,SRFIRE;
  raster SRFINA,SRFINB,SRFIMA,SRFIMB,SRFIMC,SRFIMD,SRFIME;
  raster SRFIMF,SRFIMG;
  raster PBI,PVI;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO WAVELENGTH-BASED MODELS:
# Refer to B15 & B16.
  numeric wLenCB,wLenBL,wLenGL,wLenYL,wLenRL,wLenRE;
  numeric wLenNA,wLenNB,wLenMA,wLenMB,wLenMC,wLenMD,wLenME;
  numeric wLenMF,wLenMG;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO SUN ANGLES:
  numeric sunelevang,sunaziang;
# ------------------------------------------------------------
# GET NAME OF SITE FROM THE USER: Refer to B17.
  t$  = "SITE NAME";
  p1$ = "SITE-NAME ENTRY\n";
  p2$ = "* ENTER the SITE NAME.\n";
  p3$ = "* Then, CLICK the OK Button.\n\n";
  p4$ = "SITE-NAME ENTERED:";
  p$  = p1$ + p2$ + p3$ + p4$;
  site$ = PopupString(p$,site$,t$);
  printf("           SITE NAME: %s\n",site$);
# ------------------------------------------------------------
# GET IMAGER NUMBER FROM THE USER: Refer to B17.
  p1$  = "IMAGER-NUMBER:\n";
  p2$  = "  IMAGER\n";
  p3$  = "  NUMBER:  SYSTEM NAME\n";
  p4$  = "       1:  QuickBird 2 MS\n";
  p5$  = "       2:  Ikonos 2 MS\n";
  p6$  = "       3:  OrbView 3 MS\n";
  p7$  = "       4:  Landsat 7 ETM+\n";
  p8$  = "       5:  Landsat 5 TM\n";
  p9$  = "       6:  Landsat 5 MSS\n";
  p10$ = "       7:  Landsat 4 TM\n";
  p11$ = "       8:  Landsat 4 MSS\n";
  p12$ = "       9:  Landsat 3 MSS\n";
  p13$ = "      10:  Landsat 2 MSS\n";
  p14$ = "      11:  Landsat 1 MSS\n";
  p15$ = "      12:  Terra ASTER\n";
  p16$ = "      13:  Terra MODIS\n";
  p17$ = "      14:  Aqua MODIS\n";
  p18$ = "* Either ACCEPT the Default NUMBER,\n";
  p19$ = "* Or, SELECT a Different NUMBER.\n";
  p20$ = "* Then, CLICK the OK Button.\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$,4,1,14,0);
# ------------------------------------------------------------
# GENERATE IMAGER-SPECIFIC PARAMETERS: Refer to B25.
  if (imager == 1) then begin
     imager$ = "QuickBird 2 MS";
     pBL=1; pNA=1;
     wLenBL=0.482; wLenGL=0.548; wLenRL=0.654; wLenNA=0.809;
  end
  if (imager == 2) then begin
     imager$ = "Ikonos 2 MS"; pBL=1; pNA=1;
     wLenBL=0.480; wLenGL=0.551; wLenRL=0.665; wLenNA=0.805;
  end
  if (imager == 3) then begin
     imager$ = "OrbView 3 MS"; pBL=1; pNA=1;
     wLenBL=0.482; wLenGL=0.548; wLenRL=0.654; wLenNA=0.809;
  end
  if (imager == 4) then begin
     imager$ = "Landsat 7 ETM+"; pBL=1; pNA=1;
     pMB=1; pMC=1;
     wLenBL=0.482; wLenGL=0.565; wLenRL=0.660; wLenNA=0.825;
     wLenMB=1.650; wLenMC=2.220;
  end
  if (imager == 5) then begin
     imager$ = "Landsat 5 TM";
     pBL=1; pNA=1; pMB=1; pMC=1;
     wLenBL=0.482; wLenGL=0.565; wLenRL=0.660; wLenNA=0.825;
     wLenMB=1.650; wLenMC=2.220;
  end
  if (imager == 6) then begin
     imager$ = "Landsat 5 MSS"; pRE=1; pNB=1;
     wLenGL=0.550; wLenRL=0.650; wLenRE=0.750;
     wLenNB=0.950;
  end
  if (imager == 7) then begin
     imager$ = "Landsat 4 TM"; 
     pBL=1; pNA=1; pMB=1; pMC=1;
     wLenBL=0.482; wLenGL=0.548; wLenRL=0.660; wLenNA=0.825;
     wLenMB=1.650; wLenMC=2.220;
  end
  if (imager == 8) then begin
     imager$ = "Landast 4 MSS"; pRE=1; pNB=1;
     wLenGL=0.550; wLenRL=0.650; wLenRE=0.750;
     wLenNB=0.950;
  end
  if (imager == 9) then begin
     imager$ = "Landsat 3 MSS"; pRE=1; pNB=1;
     wLenGL=0.550; wLenRL=0.650; wLenRE=0.750;
     wLenNB=0.950;
  end
  if (imager == 10) then begin
     imager$ = "Landsat 2 MSS"; pRE=1; pNB=1;
     wLenGL=0.550; wLenRL=0.650; wLenRE=0.750;
     wLenNB=0.950;
  end
  if (imager == 11) then begin
     imager$ = "Landsat 1 MSS"; pRE=1; pNB=1;
     wLenGL=0.550; wLenRL=0.650; wLenRE=0.750;
     wLenNB=0.950;
  end
  if (imager == 12) then begin
     imager$ = "Terra ASTER"; pNA=1; pMB=1; pMC=1;
     pMD=1; pME=1; pMF=1; pMG=1;
     p1$ = "VNIR or ALL OPTION:\n";
     p2$ = "   OPTION 1: VNIR ONLY\n";
     p3$ = "   OPTION 2: ALL 9 BANDS\n";
     p4$ = "OPTION ENTERED:";
     p$ = p1$ + p2$ + p3$ + p4$;
     vnir = PopupNum(p$,2,1,2,0);
     if (vnir == 1) then begin
        pMB=0; pMC=0; pMD=0; pME=0; pMF=0; pMG=0;
     end
     wLenGL=0.560; wLenRL=0.660; wLenNA=0.820;
     wLenMB=1.650; wLenMC=2.165; wLenMD=2.205; wLenME=2.260;
     wLenMF=2.330; wLenMG=2.395;
 end
  if (imager == 13) then begin
     imager$ = "Terra MODIS"; 
     pBL=1; pNA=1; pMA=1; pMB=1; pMC=1;
     wLenBL=0.469; wLenGL=0.555; wLenRL=0.645; wLenNA=0.8585;
     wLenMA=1.240; wLenMB=1.640; wLenMC=2.013;
  end
  if (imager == 14) then begin
     imager$ = "Aqua MODIS";
     pBL=1; pNA=1; pMA=1; pMB=1; pMC=1;
     wLenBL=0.469; wLenGL=0.555; wLenRL=0.645; wLenNA=0.8585;
     wLenMA=1.240; wLenMB=1.640; wLenMC=2.013;
  end
  printf("              IMAGER: %s\n\n",imager$);
# ------------------------------------------------------------
# GET fGL VALUE FROM USER: Refer to D2.
  p1$ = "GL SKY SPECTRAL-IRRADIANCE FRACTION:\n";
  p2$ = "   fGL ROLE: Determines f for Other Bands)\n";
  p3$ = "  fGL RANGE: 0.20 to 0.90\n";
  p4$ = "* Either ACCEPT the Default FRACTION,\n";
  p5$ = "* Or, ENTER a Different FRACTION.\n";
  p6$ = "* Then, CLICK the OK Button.\n\n";
  p7$ = "FRACTION ENTERED:";
  p$  = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$;
  fGL = PopupNum(p$,0.40,0.10,1.0,2);
# ------------------------------------------------------------
# ASSIGN VALUES TO f-FACTORS: Refer to D2.
  printf("SKY SPECTRAL-IRRADIANCE FRACTIONS, fBAND:\n");
  printf("   BAND   wLen  fBAND\n");
  if (pCB) then begin
     fCB = flr + (fGL - flr) * (wLenGL / wLenCB)^p;
     printf("     CB%7.3f%7.3f  (MODEL)\n",wLenCB,fCB);
  end
  if (pBL) then begin
     fBL = flr + (fGL - flr) * (wLenGL / wLenBL)^p;
     printf("     BL%7.3f%7.3f  (MODEL)\n",wLenBL,fBL);
  end
  printf("     GL%7.3f%7.3f  (INPUT)\n",wLenGL,fGL);
  if (pYL) then begin
     fYL = flr + (fGL - flr) * (wLenGL / wLenYL)^p;
     printf("     YL%7.3f%7.3f  (MODEL)\n",wLenYL,fYL);
  end
  fRL = flr + (fGL - flr) * (wLenGL / wLenRL)^p;
  printf("     RL%7.3f%7.3f  (MODEL)\n",wLenRL,fRL);
  if (pRE) then begin
     fRE = flr + (fGL - flr) * (wLenGL / wLenRE)^p;
     printf("     RE%7.3f%7.3f  (MODEL)\n",wLenRE,fRE);
  end
  if (pNA) then begin
     fNA = flr + (fGL - flr) * (wLenGL / wLenNA)^p;
     printf("     NA%7.3f%7.3f  (MODEL)\n",wLenNA,fNA);
  end
  if (pNB) then begin
     fNB = flr + (fGL - flr) * (wLenGL / wLenNB)^p;
     printf("     NB%7.3f%7.3f  (MODEL)\n",wLenNB,fNB);
  end
  if (pMA) then begin
     fMA = flr + (fGL - flr) * (wLenGL / wLenMA)^p;
     printf("     MA%7.3f%7.3f  (MODEL)\n",wLenMA,fMA);
  end
  if (pMB) then begin
     fMB = flr + (fGL - flr) * (wLenGL / wLenMB)^p;
     printf("     MB%7.3f%7.3f  (MODEL)\n",wLenMB,fMB);
  end
  if (pMC) then begin
     fMC = flr + (fGL - flr) * (wLenGL / wLenMC)^p;
     printf("     MC%7.3f%7.3f  (MODEL)\n",wLenMC,fMC);
  end
  if (pMD) then begin
     fMD = flr + (fGL - flr) * (wLenGL / wLenMD)^p;
     printf("     MD%7.3f%7.3f  (MODEL)\n",wLenMD,fMD);
  end
  if (pME) then begin
     fME = flr + (fGL - flr) * (wLenGL / wLenME)^p;
     printf("     ME%7.3f%7.3f  (MODEL)\n",wLenME,fME);
  end
  if (pMF) then begin
     fMF = flr + (fGL - flr) * (wLenGL / wLenMF)^p;
     printf("     MF%7.3f%7.3f  (MODEL)\n",wLenMC,fMC);
  end
  if (pMG) then begin
     fMG = flr + (fGL - flr) * (wLenGL / wLenMG)^p;
     printf("     MC%7.3f%7.3f  (MODEL)\n",wLenMC,fMC);
  end
  printf("\n  fBAND MODEL EXPONENT: %6.2f\n",p);
  printf("     fBAND MODEL FLOOR: %6.2f\n",flr);
# ------------------------------------------------------------
# GET SUN ELEVATION ANGLE FROM THE USER: Refer to B17.
  p1$ = "SUN-ELEVATION-ANGLE ENTRY\n";
  p2$ = "   UNITS: degrees (above the horizon)\n";
  p3$ = "  FORMAT: NN.NN\n";
  p4$ = "   RANGE: 10.00  to 90.00 degrees\n";
  p5$ = "* Either ACCEPT the Default ANGLE,\n";
  p6$ = "* Or, ENTER a Different ANGLE.\n";
  p7$ = "* Then, CLICK the OK Button.\n\n";
  p8$ = "SUN-ELEVATION-ANGLE (degrees) ENTERED:";
  p$  = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
  sunelevang = PopupNum(p$,45.18,10,90,2);
  printf("       SUN ELEV. ANGLE: %6.2f deg\n",sunelevang);
# ------------------------------------------------------------
# GET SUN AZIMUTH ANGLE FROM THE USER: Refer to EXX.
  p1$ = "SUN-AZIMUTH-ANGLE:\n";
  p2$ = "   UNITS: degrees (relative to true north)\n";
  p3$ = "  FORMAT: NNN.NN\n";
  p4$ = "   RANGE: 0.00  to 360.00 degrees\n";
  p5$ = "* Either ACCEPT the Default ANGLE,\n";
  p6$ = "* Or, ENTER a Different ANGLE.\n";
  p7$ = "* Then, CLICK the OK Button.\n\n";
  p8$ = "ANGLE ENTERED:";
  p$  = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
  sunaziang = PopupNum(p$,149.41,0,360,2);
  printf("        SUN AZI. ANGLE: %6.2f deg\n\n",sunaziang);
# ------------------------------------------------------------
# OPEN INPUT RASTERS: Refer to EXX.
  printf("INPUT RASTERS:\n");
  if (pCB) then begin
     printf(" SRFICB");
     GetInputRaster(SRFICB1);
     nlins = NumLins(SRFICB1); ncols = NumCols(SRFICB1);
     checkHisto(SRFICB1);
  end
  if (pCB) then begin
     if (pBL) then begin
        printf(" SRFIBL");
        GetInputRaster(SRFIRL1,nlins,ncols,rtype$);
     checkHisto(SRFIRL1);
     end
  end
  else begin
     if (pBL) then begin
        printf(" SRFIBL");
        GetInputRaster(SRFIBL1);
        nlins = NumLins(SRFIBL1); ncols = NumCols(SRFIBL1);
        checkHisto(SRFIBL1);
     end
  end
  printf(" SRFIGL");
  if (pBL) then begin
     GetInputRaster(SRFIGL1,nlins,ncols,rtype$);
     checkHisto(SRFIGL1);
  end
  else begin
     GetInputRaster(SRFIGL1);
     nlins = NumLins(SRFIGL1); ncols = NumCols(SRFIGL1);
     checkHisto(SRFIGL1);
  end
  if (pYL) then begin
     printf(" SRFIYL");
     GetInputRaster(SRFIYL1,nlins,ncols,rtype$);
     checkHisto(SRFIYL1);
  end
  printf(" SRFIRL");
  GetInputRaster(SRFIRL1,nlins,ncols,rtype$);
  checkHisto(SRFIRL1);
  if (pRE) then begin
     printf(" SRFIRE");
     GetInputRaster(SRFIRE1,nlins,ncols,rtype$);
     checkHisto(SRFIRE1);
  end
  if (pNA) then begin
     printf(" SRFINA");
     GetInputRaster(SRFINA1,nlins,ncols,rtype$);
     checkHisto(SRFINA1);
  end
  if (pNB) then begin
     printf(" SRFINB");
     GetInputRaster(SRFINB1,nlins,ncols,rtype$);
     checkHisto(SRFINB1);
  end
  if (pMA) then begin
     printf(" SRFIMA");
     GetInputRaster(SRFIMA1,nlins,ncols,rtype$);
     checkHisto(SRFIMA1);
  end
  if (pMB) then begin
     printf(" SRFIMB");
     GetInputRaster(SRFIMB1,nlins,ncols,rtype$);
     checkHisto(SRFIMB1);
  end
  if (pMC) then begin
     printf(" SRFIMC");
     GetInputRaster(SRFIMC1,nlins,ncols,rtype$);
     checkHisto(SRFIMC1);
  end
  if (pMD) then begin
     printf(" SRFIMD");
     GetInputRaster(SRFIMD1,nlins,ncols,rtype$);
     checkHisto(SRFIMD1);
  end
  if (pME) then begin
     printf(" SRFIME");
     GetInputRaster(SRFIME1,nlins,ncols,rtype$);
     checkHisto(SRFIME1);
  end
  if (pMF) then begin
     printf(" SRFIMF");
     GetInputRaster(SRFIMF1,nlins,ncols,rtype$);
     checkHisto(SRFIMF1);
  end
  if (pMG) then begin
     printf(" SRFIMG");
     GetInputRaster(SRFIMG1,nlins,ncols,rtype$);
     checkHisto(SRFIMG1);
  end
  printf(" PVI");
  GetInputRaster(PVI1,nlins,ncols,rtype$);
  checkHisto(PVI1);
  printf(" PBI");
  GetInputRaster(PBI1,nlins,ncols,rtype$);
     checkHisto(PBI1);
  printf(" SHADING");
  GetInputRaster(SHADING,nlins,ncols,stype$);
  checkHisto(SHADING);
# ------------------------------------------------------------
# SET UP OUTPUT RASTERS:
  printf("\n\nOUTPUT RASTERS:\n");
  if (pCB) then begin
     printf(" SRFICB");
     GetOutputRaster(SRFICB,nlins,ncols,rtype$);
     SetNull(SRFICB,0);
     CopySubobjects(SRFICB1,SRFICB,"GEOREF");
     CopySubobjects(SRFICB1,SRFICB,"CONTAB");
  end
  if (pBL) then begin
     printf(" SRFIBL");
     GetOutputRaster(SRFIBL,nlins,ncols,rtype$);
     SetNull(SRFIBL,0);
     CopySubobjects(SRFIBL1,SRFIBL,"GEOREF");
     CopySubobjects(SRFIBL1,SRFIBL,"CONTAB");
  end
  printf(" SRFIGL");
  GetOutputRaster(SRFIGL,nlins,ncols,rtype$);
  SetNull(SRFIGL,0);
  CopySubobjects(SRFIGL1,SRFIGL,"GEOREF");
  CopySubobjects(SRFIGL1,SRFIGL,"CONTAB");
  if (pYL) then begin
     printf(" SRFIYL");
     GetOutputRaster(SRFIYL,nlins,ncols,rtype$);
     SetNull(SRFIYL,0);
     CopySubobjects(SRFIYL1,SRFIYL,"GEOREF");
     CopySubobjects(SRFIYL1,SRFIYL,"CONTAB");
  end
  printf(" SRFIRL");
  GetOutputRaster(SRFIRL,nlins,ncols,rtype$);
  SetNull(SRFIRL,0);
  CopySubobjects(SRFIRL1,SRFIRL,"GEOREF");
  CopySubobjects(SRFIRL1,SRFIRL,"CONTAB");
  if (pRE) then begin
     printf(" SRFIRE");
     GetOutputRaster(SRFIRE,nlins,ncols,rtype$);
     SetNull(SRFIRE,0);
     CopySubobjects(SRFIRE1,SRFIRE,"GEOREF");
     CopySubobjects(SRFIRE1,SRFIRE,"CONTAB");
  end
  if (pNA) then begin
     printf(" SRFINA");
     GetOutputRaster(SRFINA,nlins,ncols,rtype$);
     SetNull(SRFINA,0); 
     CopySubobjects(SRFINA1,SRFINA,"GEOREF");
     CopySubobjects(SRFINA1,SRFINA,"CONTAB");
  end
  if (pNB) then begin
     printf(" SRFINB");
     GetOutputRaster(SRFINB,nlins,ncols,rtype$);
     SetNull(SRFINB,0);
     CopySubobjects(SRFINB1,SRFINB,"GEOREF");
     CopySubobjects(SRFINB1,SRFINB,"CONTAB");
  end
  if (pMA) then begin
     printf(" SRFIMA");
     GetOutputRaster(SRFIMA,nlins,ncols,rtype$);
     SetNull(SRFIMA,0);
     CopySubobjects(SRFIMA1,SRFIMA,"GEOREF");
     CopySubobjects(SRFIMA1,SRFIMA,"CONTAB");
  end
  if (pMB) then begin
     printf(" SRFIMB");
     GetOutputRaster(SRFIMB,nlins,ncols,rtype$);
     SetNull(SRFIMB,0);
     CopySubobjects(SRFIMB1,SRFIMB,"GEOREF");
     CopySubobjects(SRFIMB1,SRFIMB,"CONTAB");
  end
  if (pMC) then begin
     printf(" SRFIMC");
     GetOutputRaster(SRFIMC,nlins,ncols,rtype$);
     SetNull(SRFIMC,0);
     CopySubobjects(SRFIMC1,SRFIMC,"GEOREF");
     CopySubobjects(SRFIMC1,SRFIMC,"CONTAB");
  end
  if (pMD) then begin
     printf(" SRFIMD");
     GetOutputRaster(SRFIMD,nlins,ncols,rtype$);
     SetNull(SRFIMD,0);
     CopySubobjects(SRFIMD1,SRFIMD,"GEOREF");
     CopySubobjects(SRFIMD1,SRFIMD,"CONTAB");
  end
  if (pME) then begin
     printf(" SRFIME");
     GetOutputRaster(SRFIME,nlins,ncols,rtype$);
     SetNull(SRFIME,0);
     CopySubobjects(SRFIME1,SRFIME,"GEOREF");
     CopySubobjects(SRFIME1,SRFIME,"CONTAB");
  end
  if (pMF) then begin
     printf(" SRFIMF");
     GetOutputRaster(SRFIMF,nlins,ncols,rtype$);
     SetNull(SRFIMF,0);
     CopySubobjects(SRFIMF1,SRFIMF,"GEOREF");
     CopySubobjects(SRFIMF1,SRFIMF,"CONTAB");
  end
  if (pMG) then begin
     printf(" SRFIMG");
     GetOutputRaster(SRFIMG,nlins,ncols,rtype$);
     SetNull(SRFIMG,0);
     CopySubobjects(SRFIMG1,SRFIMG,"GEOREF");
     CopySubobjects(SRFIMG1,SRFIMG,"CONTAB");
  end
  printf(" PVI");
  GetOutputRaster(PVI,nlins,ncols,rtype$);
  SetNull(PVI,0);
  CopySubobjects(PVI1,PVI,"GEOREF");
  CopySubobjects(PVI1,PVI,"CONTAB");
  CloseRaster(PVI1);
  printf(" PBI");
  GetOutputRaster(PBI,nlins,ncols,rtype$);
  SetNull(PBI,0);
  CopySubobjects(PBI1,PBI,"GEOREF");
  CopySubobjects(PBI1,PBI,"CONTAB");
  CloseRaster(PBI1);
# ------------------------------------------------------------
# PRODUCE OUTPUT RASTERS:  Refer to EXX.
  printf("\n\n");
  printf("CORRECTING SRFT RASTERS:\n");
  array numeric tc[256];
  for i=0 to 254 begin
     tc[i] = i / 180;
  end
# CB DATA:
  if (pCB) then begin
     printf(" SRFICB");
     for each SRFICB1 begin
        iip = 1; isp = 1;
        if (IsNull(SRFICB1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1   = SRFICB1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fCB)+fCB));
           if (srfi < 1)    then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFICB = srfi;
        end
     end
     CreateHistogram(SRFICB,0); CreatePyramid(SRFICB,0);
     CloseRaster(SRFICB1); CloseRaster(SRFICB);
  end
# BL DATA:
  if (pBL) then begin
     printf(" SRFIBL");
     for each SRFIBL1 begin
        iip = 1; isp = 1;
        if (IsNull(SRFIBL1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1   = SRFIBL1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fBL)+fBL));
           if (srfi < 1)    then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFIBL = srfi;
        end
     end
     CreateHistogram(SRFIBL,0); CreatePyramid(SRFIBL,0);
     CloseRaster(SRFIBL1); CloseRaster(SRFIBL);
  end
# GL DATA: 
  printf(" SRFIGL");
  for each SRFIGL1 begin
     iip = 1; isp = 1;
     if (IsNull(SRFIGL1)) then iip = 0;
     if (IsNull(SHADING)) then isp = 0;
     ipp = iip * isp;
     if (ipp) then begin
        srfi1 = SRFIGL1; s = SHADING; s = tc[s];
        srfi = round(srfi1/(s*(1-fGL)+fGL));
        if (srfi < 1)    then srfi = 1;
        if (srfi > srfimax) then srfi = srfimax;
        SRFIGL = srfi;
     end
  end
  CreateHistogram(SRFIGL,0); CreatePyramid(SRFIGL,0);
  CloseRaster(SRFIGL1); CloseRaster(SRFIGL);
# YL DATA:
  if (pYL) then begin
     printf(" SRFIYL");
     for each SRFIYL1 begin
        iip = 1; isp = 1;
        if (IsNull(SRFIYL1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1 = SRFIYL1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fYL)+fYL));
           if (srfi < 1)    then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFIYL = srfi;
        end
     end
     CreateHistogram(SRFIYL,0); CreatePyramid(SRFIYL,0);
     CloseRaster(SRFIYL1); CloseRaster(SRFIYL);
  end
# RL DATA:
  printf(" SRFIRL");
  for each SRFIRL1 begin
     iip = 1; isp = 1;
     if (IsNull(SRFIRL1)) then iip = 0;
     if (IsNull(SHADING)) then isp = 0;
     ipp = iip * isp;
     if (ipp) then begin
        srfi1 = SRFIRL1; s = SHADING; s = tc[s];
        srfi = round(srfi1/(s*(1-fRL)+fRL));
        if (srfi < 1)    then srfi = 1;
        if (srfi > srfimax) then srfi = srfimax;
        SRFIRL = srfi;
     end
  end
  CreateHistogram(SRFIRL,0); CreatePyramid(SRFIRL,0);
  CloseRaster(SRFIRL1);
# RE DATA:
  if (pRE) then begin
     printf(" SRFIRE");
     for each SRFIRE1 begin 
        iip = 1; isp = 1;
        if (IsNull(SRFIRE1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1 = SRFIRE1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fRE)+fRE));
           if (srfi < 1)    then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFIRE = srfi;
        end
     end
     CreateHistogram(SRFIRE,0); CreatePyramid(SRFIRE,0);
     CloseRaster(SRFIRE1); CloseRaster(SRFIRE);
   end
# NA DATA:
  if (pNA) then begin
     printf(" SRFINA");
     for each SRFINA1 begin
        iip = 1; isp = 1;
        if (IsNull(SRFINA1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1 = SRFINA1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fNA)+fNA));
           if (srfi < 1)    then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFINA = srfi;
        end
     end
     CreateHistogram(SRFINA,0); CreatePyramid(SRFINA,0);
     CloseRaster(SRFINA1);
  end
# NB DATA:
  if (pNB) then begin
     printf(" SRFINB");
     for each SRFINB1 begin
        iip = 1; isp = 1;
        if (IsNull(SRFINB1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1 = SRFINB1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fNB)+fNB));
           if (srfi < 1)    then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFINB = srfi;
        end
     end
     CreateHistogram(SRFINB,0); CreatePyramid(SRFINB,0);
     CloseRaster(SRFINB1);
  end
# MA DATA:
  if (pMA) then begin
     printf(" SRFIMA");
     for each SRFIMA1 begin
        iip = 1; isp = 1;
        if (IsNull(SRFIMA1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1 = SRFIMA1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fMA)+fMA));
           if (srfi < 1)    then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFIMA = srfi;
        end
     end
     CreateHistogram(SRFIMA,0); CreatePyramid(SRFIMA,0);
     CloseRaster(SRFIMA1); CloseRaster(SRFIMA);
  end
# MB DATA:
  if (pMB) then begin
     printf(" SRFIMB");
     for each SRFIMB1 begin
        iip = 1; isp = 1;
        if (IsNull(SRFIMB1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1 = SRFIMB1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fMB)+fMB));
           if (srfi < 1)    then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFIMB = srfi;
        end
     end
     CreateHistogram(SRFIMB,0); CreatePyramid(SRFIMB,0);
     CloseRaster(SRFIMB1); CloseRaster(SRFIMB);
  end
# MC DATA:
  if (pMC) then begin
     printf(" SRFIMC");
     for each SRFIMC1 begin
        iip = 1; isp = 1;
        if (IsNull(SRFIMC1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1 = SRFIMC1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fMC)+fMC));
           if (srfi < 1) then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFIMC = srfi;
        end
     end
     CreateHistogram(SRFIMC,0); CreatePyramid(SRFIMC,0);
     CloseRaster(SRFIMC1); CloseRaster(SRFIMC);
  end
# MD DATA:
  if (pMD) then begin
     printf(" SRFIMD");
     for each SRFIMD1 begin
        iip = 1; isp = 1;
        if (IsNull(SRFIMD1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1 = SRFIMD1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fMD)+fMD));
           if (srfi < 1) then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFIMD = srfi;
        end
     end
     CreateHistogram(SRFIMD,0); CreatePyramid(SRFIMD,0);
     CloseRaster(SRFIMD1); CloseRaster(SRFIMD);
  end
# ME DATA:
  if (pME) then begin
     printf(" SRFIME");
     for each SRFIME1 begin
        iip = 1; isp = 1;
        if (IsNull(SRFIME1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1 = SRFIME1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fME)+fME));
           if (srfi < 1) then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFIME = srfi;
        end
     end
     CreateHistogram(SRFIME,0); CreatePyramid(SRFIME,0);
     CloseRaster(SRFIME1); CloseRaster(SRFIME);
  end
# MF DATA:
  if (pMF) then begin
     printf(" SRFIMF");
     for each SRFIMF1 begin
        iip = 1; isp = 1;
        if (IsNull(SRFIMF1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1 = SRFIMF1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fMF)+fMF));
           if (srfi < 1) then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFIMF = srfi;
        end
     end
     CreateHistogram(SRFIMF,0); CreatePyramid(SRFIMF,0);
     CloseRaster(SRFIMF1); CloseRaster(SRFIMF);
  end
# MG DATA:
  if (pMG) then begin
     printf(" SRFIMG");
     for each SRFIMG1 begin
        iip = 1; isp = 1;
        if (IsNull(SRFIMG1)) then iip = 0;
        if (IsNull(SHADING)) then isp = 0;
        ipp = iip * isp;
        if (ipp) then begin
           srfi1 = SRFIMG1; s = SHADING; s = tc[s];
           srfi = round(srfi1/(s*(1-fMG)+fMG));
           if (srfi < 1) then srfi = 1;
           if (srfi > srfimax) then srfi = srfimax;
           SRFIMG = srfi;
        end
     end
     CreateHistogram(SRFIMG,0); CreatePyramid(SRFIMG,0);
     CloseRaster(SRFIMG1); CloseRaster(SRFIMG);
  end
  CloseRaster(SHADING);
  printf("\n\n");
# ------------------------------------------------------------
# PRODUCE VALUES FOR PVI & PBI RASTERS:
  printf("PRODUCING PVI & PBI RASTERS.");
  srfiNRint = 254; bslineslope = 1.086;
  pbioff = 100; pvioff = 1000;
  pfac = 0.2723659; maxpvipbi = 3000;
  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 begin
           srfiNR = SRFINB;
        end
        tsrfiNR = srfiNR - srfiNRint;
        pbif = srfiRL * cosang - srfiNR * sinang;
        pbi  = round(pfac * pbif)+ pbioff;
        pvif = srfiRL * sinang + srfiNR * cosang;
        pvi  = round(pfac * pvif) + pvioff;
        if (pbi < 1) then pbi = 1;
        if (pbi > maxpvipbi) then pbi = maxpvipbi;
        if (pvi < 1) then pvi = 1;
        if (pvi > maxpvipbi) then pvi = maxpvipbi;
        PBI = pbi; PVI = pvi;
     end
  end
  CreateHistogram(PBI,0); CreateHistogram(PVI,0);
  CreatePyramid(PBI,0); CreatePyramid(PVI,0);
  CloseRaster(PBI); CloseRaster(PVI);
  if (pNA) then CloseRaster(SRFINA);
  if (pNB) then CloseRaster(SRFINB);
  CloseRaster(SRFIRL);
# ------------------------------------------------------------
  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.");