# ------------------------------------------------------------
# TASCAP.sml
# ------------------------------------------------------------
# SET WARNING LEVEL:
  $warnings 3
# ------------------------------------------------------------
# DEFINE PROCEDURE writeTitle:
# PURPOSE: WRITES TITLE & AUTHOR INFO TO CONSOLE WINDOW.
  proc writeTitle() begin
     printf("TASCAP.sml:\n");
     printf("         VERSION: November 15, 2005\n");
     printf("         PURPOSE: CONVERT SRFI DATA ");
     printf("TO TASSELED-CAP DATA\n");
     printf("                  WITH AUTOMATIC CONTRAST ");
     printf("LOOK-UP TABLES.\n");
     printf("         DETAILS: FAQs_by_Jack A & F\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
# ------------------------------------------------------------
# DECLARE VARIABLES
  numeric nbands,bp0,bp1,bp2,bp3,bp4,bp5,bp6,bp7,bp8,bp9;
  numeric i,j,nbandsm1,m;
  numeric itc,nTC,tc0,tc1,tc2,tc3,tc4,tc5,tc6,tc7,tc8,tc9;
  numeric lin,col,nlins,ncols,fsize;
  numeric tcn,tcnp1,isnull,lintcp,coltcp;
  array lintcp[10]; array coltcp[10];
  string bp0$,bp1$,bp2$,bp3$,bp4$,bp5$,bp6$,bp7$,bp8$;
  string stype$,dtype$,ttype$;
  dtype$ = "16-bit unsigned";
  ttype$ = "16-bit signed"; tcn = -32000; tcnp1 = tcn + 1;
  numeric d,d2,sd2,mag,f,sf2,ds,tc;
  numeric dot1,dot2,dot3,dot4,dot5,dot6,dot7,dot8,dot9;
  numeric c1DS0,c1DS1,c1DS2,c1DS3,c1DS4;
  numeric c1DS5,c1DS6,c1DS7,c1DS8;
  numeric c2DS0,c2DS1,c2DS2,c2DS3,c2DS4;
  numeric c2DS5,c2DS6,c2DS7,c2DS8;
  numeric c3DS0,c3DS1,c3DS2,c3DS3,c3DS4;
  numeric c3DS5,c3DS6,c3DS7,c3DS8;
  numeric pDS0,pDS1,pDS2,pDS3,pDS4,pDS5,pDS6,pDS7,pDS8;
  numeric c1TC1,c1TC2,c1TC3,c1TC4,c1TC5,c1TC6,c1TC7,c1TC8;
  numeric c2TC1,c2TC2,c2TC3,c2TC4,c2TC5,c2TC6,c2TC7,c2TC8;
  numeric c3TC1,c3TC2,c3TC3,c3TC4,c3TC5,c3TC6,c3TC7,c3TC8;
  numeric pTC1,pTC2,pTC3,pTC4,pTC5,pTC6,pTC7,pTC8;
  numeric hnmDS0,hnmDS1,hnmDS2,hnmDS3,hnmDS4;
  numeric hnmDS5,hnmDS6,hnmDS7,hnmDS8;
  numeric hnDS0,hnDS1,hnDS2,hnDS3,hnDS4;
  numeric hnDS5,hnDS6,hnDS7,hnDS8;
  numeric nvDS0,nvDS1,nvDS2,nvDS3,nvDS4;
  numeric nvDS5,nvDS6,nvDS7,nvDS8;
  numeric hnmTC1,hnmTC2,hnmTC3,hnmTC4;
  numeric hnmTC5,hnmTC6,hnmTC7,hnmTC8;
  numeric hnTC1,hnTC2,hnTC3,hnTC4;
  numeric hnTC5,hnTC6,hnTC7,hnTC8;
  numeric nvTC1,nvTC2,nvTC3,nvTC4;
  numeric nvTC5,nvTC6,nvTC7,nvTC8;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO CONTRAST LOOKUP PROCESS: 
  class CONTRAST smlContrast;
  string clut1$,clut2$,clut3$;
  array numeric v[131072];
  clut1$  = "exp"; clut2$  = "Exponential";
  clut3$  = "Exponential contrast table";
  smlContrast.OutputLowerLimit =   1;
  smlContrast.OutputUpperLimit = 255;
# ------------------------------------------------------------
# DEFINE STRING VARIABLES FOR POPUP WINDOWS:
  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$;
# ------------------------------------------------------------
# DECLARE ARRAYS
  array bp0[10]; array bp1[10]; array bp2[10]; array bp3[10];
  array bp4[10]; array bp5[10]; array bp6[10]; array bp7[10];
  array bp8[10]; array bp9[10]; array tc0[10]; array tc1[10];
  array tc2[10]; array tc3[10]; array tc4[10]; array tc5[10];
  array tc6[10]; array tc7[10]; array tc8[10]; array tc9[10];
  array numeric srfi[10];
# ------------------------------------------------------------
# CLEAR CONSOLE WINDOW & REQUEST REPOSITIONING:
  clear();
  p1$ = "CONSOLE-WINDOW ADJUSTMENT:\n";
  p2$ = "* REPOSITION the CONSOLE WINDOW.\n";
  p3$ = "* Then, CLICK OK.";
  p$  = p1$ + p2$ + p3$;
  PopupMessage(p$);
# ------------------------------------------------------------
# SELECT METHOD (OUT OF 3 CHOICES): m
  p1$ = "METHOD-NUMBER ENTRY:\n";
  p2$ = "   CHOICES:\n";
  p3$ = "      1: Process LANDSAT SRFItoa DATA (6 Bands)\n";
  p4$ = "         Using USGS TASCAP Coefficients.\n";
  p5$ = "      2: Process with PRE-DEFINED bp VALUES.\n";
  p6$ = "      3: Specify KEY PIXELS (lin,col) to Fetch\n";
  p7$ = "         bp Values from Input SRFT Rasters.\n";
  p8$ = "METHOD-NUMBER ENTERED:";
  p$  = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
  m = PopupNum(p$,1,1,3,0);
# ------------------------------------------------------------
# ASSIGN VALUES TO dtype$, ttype$, tcn, tcnp1, nbands & nTC:
  dtype$ = "16-bit unsigned";
  ttype$ = "16-bit signed"; tcn = -32000; tcnp1 = tcn + 1;
  nbands = 6; nTC = 4;
# ============================================================
# DEFAULT PARAMETERS FOR METHOD 1: Landsat 6-Band MS Data.
# Classic TASCAP for SRFI Rasters: BL, GL, RL, NA, MB, MC:
# WARNING: Only Applicable to Landsat SRFItoa Raster Values.
  if (m == 1) then begin
     bp0$ = "BLACK"; bp1$ = "Brightness";
     bp2$ = "Greenness"       ; bp3$ = "Wetness"   ;
     bp4$ = "Haze"            ;
     tc0[1]=     0.0; tc0[2]=     0.0; tc0[3]=     0.0;
     tc0[4]=     0.0; tc0[5]=     0.0; tc0[6]=     0.0;
     tc1[1]=  0.3561; tc1[2]=  0.3972; tc1[3]=  0.3904;
     tc1[4]=  0.6966; tc1[5]=  0.2286; tc1[6]=  0.1596;
     tc2[1]= -0.3344; tc2[2]= -0.3544; tc2[3]= -0.4556;
     tc2[4]=  0.6966; tc2[5]= -0.0242; tc2[6]= -0.2630;
     tc3[1]=  0.2626; tc3[2]=  0.2141; tc3[3]=  0.0926;
     tc3[4]=  0.0656; tc3[5]= -0.7629; tc3[6]= -0.5388;
     tc4[1]=  0.0805; tc4[2]= -0.0498; tc4[3]=  0.1950;
     tc4[4]= -0.1327; tc4[5]=  0.5752; tc4[6]= -0.7775;
  end
# ============================================================
# DEFAULT PARAMETERS FOR METHOD 2:
# Default TASCAP for SRFI Rasters: BL, GL, RL, NA, (MB, MC):
# METHOD 2 Requires No Input Key Inputs from User:
# NOTE: METHOD 2 can be Applied to 3 to 6 Bands (SRFI):
  if (m == 2) then begin
     bp0$="Dark Soil       "; bp1$="Bright Soil";
     bp2$="Dense Green Veg."; bp3$="Water"      ;
     bp4$="Yellow Veg."     ;
     t$  = "n-SPACE ORIGIN TYPE";
     p1$ = "SELECT n-SPACE ORIGIN TYPE:\n";
     p2$ = "     ENTER 'BLACK' or ACCEPT DEFAULT:\n";
     p3$ = "n-SPACE ORIGIN-TYPE ENTERED:";
     p$ = p1$ + p2$ + p3$;
     bp0$ = PopupString(p$,bp0$,t$);
     bp0[1]= 190; bp0[2]= 326; bp0[3]= 100;
     bp0[4]= 200; bp0[5]=142; bp0[6]=  276;
     if (bp0$ == "BLACK") then begin
        bp0[1]=0; bp0[2]=0; bp0[3]=0;
        bp0[4]=0; bp0[5]=0; bp0[6]=0;
     end
     bp1[1]= 945; bp1[2]=1157; bp1[3]=1746;
     bp1[4]=1998; bp1[5]=2688; bp1[6]=1967;
     bp2[1]= 400; bp2[2]= 600; bp2[3]= 300;
     bp2[4]=6000; bp2[5]=2000; bp2[6]= 600;
     bp3[1]= 400; bp3[2]= 500; bp3[3]= 200;
     bp3[4]= 189; bp3[5]= 105; bp3[6]=  70;
     bp4[1]= 450; bp4[2]=1000; bp4[3]=1400;
     bp4[4]=4500; bp4[5]=2500; bp4[6]= 800;
  end
# ============================================================
# METHOD 3:  User must Provide Raster lin & col Locations that
# TASCAP.sml Will Use to Fetch Related bp Values.
# ============================================================
# 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
# ------------------------------------------------------------
# func GlobalPtile(Rl,ptilel):
  func GlobalPtile(raster Rl, numeric ptilel,hnl,nvl) begin
     local numeric funcval,rl,roffl;
     local numeric nlinsl,ncolsl,linl,coll,il;
     local numeric countl,suml,flagl,cfql;
     roffl = -65535; countl = 0;
     if (hnl == 1) then begin 
        for each Rl begin
           rl = Rl;
           if (rl <> nvl) then countl = countl + 1;
        end
     end
     if (hnl == 0) then begin
        nlinsl = NumLins(Rl); ncolsl = NumCols(Rl);
        countl = nlinsl * ncolsl;
     end
     for il = 0 to 131071 v[il] = 0;
     for each Rl begin
        rl = Rl;
        if (hnl == 1) then begin
           if (rl <> nvl) then begin
              il = rl - roffl; v[il] = v[il] + 1;
           end
        end
        if (hnl == 0) then begin
           il = rl - roffl; v[il] = v[il] + 1;
        end
     end
     suml = v[0]; flagl = 0;
     for il=1 to 131071 begin
        if (flagl == 0) then begin
           suml = suml + v[il]; cfql = suml * 100 / countl;
           if (cfql >= ptilel) then begin
              flagl = 1; funcval = il + roffl;
           end
        end
     end
     return (funcval);
  end
# ------------------------------------------------------------
# DEFINE func expPower(dnLl,dnMl,dnHl):
  func expPower(dnLl,dnMl,dnHl) begin
     local numeric delLHl,delLMl,sr1l,sr2l,exppl;
     sr1l = 0.3; sr2l = 2 * (1 - sr1l);
     delLHl = dnHl - dnLl; delLMl = dnMl - dnLl;
     exppl = sr1l + sr2l * delLMl / delLHl;
     return (exppl);
  end
# ------------------------------------------------------------
# WRITE TITLE & AUTHOR INFORMATION:
  writeTitle();
# ------------------------------------------------------------
# WRITE METHOD TO CONSOLE WINDOW:
  if (m == 1) then begin
     printf("METHOD 1: Use USGS Coefficients\n");
     printf("          Applied to Landsat SRFItoa Data.\n\n");
  end
  if (m == 2) then begin
     printf("METHOD 2: Use Default TASCAP Coefficients\n");
     printf("          Applied to SRFIsfc Data.\n\n");
  end
  if (m == 3) then begin
     printf("METHOD 3: Use Customized TASCAP Coefficients\n");
     printf("          Applied to SRFIsfc.\n");
     printf("          Customized Names with Associated\n");
     printf("          SRFI Image LIN & COL Locations.\n\n");
  end
# ------------------------------------------------------------
# SPECIFY nbands & nTC PARAMETERS:
  if (m == 2) then begin
     p1$ = "NUMBER OF INPUT SRFI RASTERS ENTRY:\n";
     p2$ = "   RANGE: 4 to 6\n";
     p3$ = "NUMBER ENTERED:";
     p$  = p1$ + p2$ + p3$;
     nbands = PopupNum(p$,nbands,4,6,0);
     nbandsm1 = nbands - 1;
  end
  if (m == 2) then begin
     if (nTC > (nbands - 1)) then nTC = nbands - 1;
     p1$ = "NUMBER OF OUTPUT RASTER PAIRS ENTRY:\n";
     p2$ = "     RANGE: 2 to (NumInputBands - 1)\n";
     p3$ = "NUMBER ENTERED:";
     p$  = p1$ + p2$ + p3$;
     nTC = PopupNum(p$,nTC,2,nbandsm1,0);
  end
  if (m == 3) then begin
     p1$ = "NUMBER OF INPUT SRFI RASTERS ENTRY:\n";
     p2$ = "   RANGE: 3 to 9\n";
     p3$ = "NUMBER ENTERED:";
     p$  = p1$ + p2$ + p3$;
     nbands = PopupNum(p$,nbands,3,9,0);
     nbandsm1 = nbands - 1;
  end
  printf("NUMBERS OF INPUT RASTERS & OUTPUT RASTERS:\n");
  printf("   %1d Sets of SRFI INPUT Rasters\n",nbands);
  if (m == 3) then begin
     if (nTC > (nbands - 1)) then nTC = nbands - 1;
     p1$ = "NUMBER OF OUTPUT RASTER PAIRS ENTRY:\n";
     p2$ = "   RANGE: 2 to (NumInputBands - 1)\n";
     p3$ = "NUMBER ENTERED:";
     p$  = p1$ + p2$ + p3$;
     nTC = PopupNum(p$,nTC,2,nbandsm1,0);
  end
  printf("   %1d Pairs of TC & DS OUTPUT Rasters\n\n",nTC);
# ------------------------------------------------------------
# OPEN INPUT RASTERS:
  raster SRFI1,SRFI2,SRFI3,SRFI4,SRFI5;
  raster SRFI6,SRFI7,SRFI8,SRFI9;
  GetInputRaster(SRFI1);
  nlins = NumLins(SRFI1); ncols = NumCols(SRFI1);
  stype$ = RastType(SRFI1);
  checkHisto(SRFI1);
  GetInputRaster(SRFI2,nlins,ncols,stype$);
  checkHisto(SRFI2);
  GetInputRaster(SRFI3,nlins,ncols,stype$);
  checkHisto(SRFI3);
  if (nbands>3) then begin
     GetInputRaster(SRFI4,nlins,ncols,stype$);
     checkHisto(SRFI4);
  end
  if (nbands>4) then begin
     GetInputRaster(SRFI5,nlins,ncols,stype$);
     checkHisto(SRFI5);
  end
  if (nbands>5) then begin
     GetInputRaster(SRFI6,nlins,ncols,stype$);
     checkHisto(SRFI6);
  end
  if (nbands>6) then begin
     GetInputRaster(SRFI7,nlins,ncols,stype$);
     checkHisto(SRFI7);
  end
  if (nbands>7) then begin
     GetInputRaster(SRFI8,nlins,ncols,stype$);
     checkHisto(SRFI8);
  end
  if (nbands>8) then begin
     GetInputRaster(SRFI9,nlins,ncols,stype$);
     checkHisto(SRFI9);
  end
# ------------------------------------------------------------
# GET USER-PROVIDED LIN & COL PARAMETERS:
  if (m == 3) then begin
     t$   = "METHOD-3 PARAMETER";
     bp0$ = "bp0 NAME";
     p$   = "ENTER bp0 NAME:";
     bp0$ = PopupString(p$,bp0$,t$);
     if (bp0$ == "BLACK") then begin
        lintcp[0] = 1;
        coltcp[0] = 1;
     end
     p1$ = "METHOD-3 LINE & COLUMN PARAMETERS ENTRY:\n";
     if (bp0$ <> "BLACK") then begin
        p2$  = "FOR " + bp0$ + ", ENTER LIN:";
        p$   = p1$ + p2$;
        lintcp[0] = PopupNum(p$,1,1,100000,0);
        p2$  = "FOR " + bp0$ + ", ENTER COL:";
        p$   = p1$ + p2$;
        coltcp[0] = PopupNum(p$,1,1,100000,0);
     end
     bp1$ = "bp1 NAME";
     p$   = "ENTER bp1 NAME:";
     bp1$ = PopupString(p$,bp1$,t$);
     p2$  = "FOR " + bp1$ + ", ENTER LIN:";
     p$   = p1$ + p2$;
     lintcp[1] = PopupNum(p$,1,1,100000,0);
     p2$  = "FOR " + bp1$ + ", ENTER COL:";
     p$   = p1$ + p2$;
     coltcp[1] = PopupNum(p$,1,1,100000,0);
     bp2$ = "bp2 NAME";
     p$   = "ENTER bp2 NAME:";
     bp2$ = PopupString(p$,bp2$,t$);
     p2$  = "FOR " + bp2$ + ", ENTER LIN:";
     p$   = p1$ + p2$;
     lintcp[2] = PopupNum(p$,1,1,100000,0);
     p2$  = "FOR " + bp2$ + ", ENTER COL:";
     p$   = p1$ + p2$;
     coltcp[2] = PopupNum(p$,1,1,100000,0);
     if (nTC > 2) then begin
        bp3$ = "bp3 NAME";
        p$   = "ENTER bp3 NAME:";
        bp3$ = PopupString(p$,bp3$,t$);
        p2$  = "FOR " + bp3$ + ", ENTER LIN:";
        p$   = p1$ + p2$;
        lintcp[3] = PopupNum(p$,1,1,100000,0);
        p2$  = "FOR " + bp3$ + ", ENTER COL:";
        p$   = p1$ + p2$;
        coltcp[3] = PopupNum(p$,1,1,100000,0);
     end
     if (nTC > 3) then begin
        bp4$ = "bp4 NAME";
        p$   = "ENTER bp4 NAME:";
        bp4$ = PopupString(p$,bp4$,t$);
        p2$  = "FOR " + bp4$ + ", ENTER LIN:";
        p$   = p1$ + p2$;
        lintcp[4] = PopupNum(p$,1,1,100000,0);
        p2$  = "FOR " + bp4$ + ", ENTER COL:";
        p$   = p1$ + p2$;
        coltcp[4] = PopupNum(p$,1,1,100000,0);
     end
     if (nTC > 4) then begin
        bp5$ = "bp5 NAME";
        p$   = "ENTER bp5 NAME:";
        bp5$ = PopupString(p$,bp5$,t$);
        p2$  = "FOR " + bp5$ + ", ENTER LIN:";
        p$   = p1$ + p2$;
        lintcp[5] = PopupNum(p$,1,1,100000,0);
        p2$  = "FOR " + bp5$ + ", ENTER COL:";
        p$   = p1$ + p2$;
        coltcp[5] = PopupNum(p$,1,1,100000,0);
     end
     if (nTC > 5) then begin
        bp6$ = "bp6 NAME";
        p$   = "ENTER bp6 NAME:";
        bp6$ = PopupString(p$,bp6$,t$);
        p2$  = "FOR " + bp6$ + ", ENTER LIN:";
        p$   = p1$ + p2$;
        lintcp[6] = PopupNum(p$,1,1,100000,0);
        p2$  = "FOR " + bp6$ + ", ENTER COL:";
        p$   = p1$ + p2$;
        coltcp[6] = PopupNum(p$,1,1,100000,0);
     end
     if (nTC > 6) then begin
        bp7$ = "bp7 NAME";
        p$   = "ENTER bp7 NAME:";
        bp7$ = PopupString(p$,bp7$,t$);
        p2$  = "FOR " + bp7$ + ", ENTER LIN:";
        p$   = p1$ + p2$;
        lintcp[7] = PopupNum(p$,1,1,100000,0);
        p2$  = "FOR " + bp7$ + ", ENTER COL:";
        p$   = p1$ + p2$;
        coltcp[7] = PopupNum(p$,1,1,100000,0);
     end
     if (nTC > 7) then begin
        bp8$ = "bp8 NAME";
        p$   = "ENTER bp8 NAME:";
        bp8$ = PopupString(p$,bp8$,t$);
        p2$  = "FOR " + bp8$ + ", ENTER LIN:";
        p$   = p1$ + p2$;
        lintcp[8] = PopupNum(p$,1,1,100000,0);
        p2$  = "FOR " + bp8$ + ", ENTER COL:";
        p$   = p1$ + p2$;
        coltcp[8] = PopupNum(p$,1,1,100000,0);
     end
  end
# ------------------------------------------------------------
# GET bpn[nband] VALUES:
  for i=0 to nTC begin
     if (m == 3) then begin
        lin = lintcp[i]; col = coltcp[i];
        if (i==0) then begin
           bp0[1] = 0; bp0[2] = 0; bp0[3] = 0;
           bp0[4] = 0; bp0[5] = 0; bp0[6] = 0;
           bp0[7] = 0; bp0[8] = 0; bp0[9] = 0;
           if (bp0$ <> "BLACK") then begin
              bp0[1] = SRFI1[lin,col]; bp0[2] = SRFI2[lin,col];
              bp0[3] = SRFI3[lin,col];
              if (nbands > 3) then bp0[4] = SRFI4[lin,col];
              if (nbands > 4) then bp0[5] = SRFI5[lin,col];
              if (nbands > 5) then bp0[6] = SRFI6[lin,col];
              if (nbands > 6) then bp0[7] = SRFI7[lin,col];
              if (nbands > 7) then bp0[8] = SRFI8[lin,col];
              if (nbands > 8) then bp0[9] = SRFI9[lin,col];
           end
        end
        if (i==1) then begin
           bp1[1] = SRFI1[lin,col]; bp1[2] = SRFI2[lin,col];
           bp1[3] = SRFI3[lin,col];
           if (nbands > 3) then bp1[4] = SRFI4[lin,col];
           if (nbands > 4) then bp1[5] = SRFI5[lin,col];
           if (nbands > 5) then bp1[6] = SRFI6[lin,col];
           if (nbands > 6) then bp1[7] = SRFI7[lin,col];
           if (nbands > 7) then bp1[8] = SRFI8[lin,col];
           if (nbands > 8) then bp1[9] = SRFI9[lin,col];
        end
        if (i==2) then begin
           bp2[1] = SRFI1[lin,col]; bp2[2] = SRFI2[lin,col];
           bp2[3] = SRFI3[lin,col];
           if (nbands > 3) then bp2[4] = SRFI4[lin,col];
           if (nbands > 4) then bp2[5] = SRFI5[lin,col];
           if (nbands > 5) then bp2[6] = SRFI6[lin,col];
           if (nbands > 6) then bp2[7] = SRFI7[lin,col];
           if (nbands > 7) then bp2[8] = SRFI8[lin,col];
           if (nbands > 8) then bp2[9] = SRFI9[lin,col];
        end
        if (i==3) then begin
           bp3[1] = SRFI1[lin,col]; bp3[2] = SRFI2[lin,col];
           bp3[3] = SRFI3[lin,col];
           if (nbands > 3) then bp3[4] = SRFI4[lin,col];
           if (nbands > 4) then bp3[5] = SRFI5[lin,col];
           if (nbands > 5) then bp3[6] = SRFI6[lin,col];
           if (nbands > 6) then bp3[7] = SRFI7[lin,col];
           if (nbands > 7) then bp3[8] = SRFI8[lin,col];
           if (nbands > 8) then bp3[9] = SRFI9[lin,col];
        end
        if (i==4) then begin
           bp4[1] = SRFI1[lin,col]; bp4[2] = SRFI2[lin,col];
           bp4[3] = SRFI3[lin,col];
           if (nbands > 3) then bp4[4] = SRFI4[lin,col];
           if (nbands > 4) then bp4[5] = SRFI5[lin,col];
           if (nbands > 5) then bp4[6] = SRFI6[lin,col];
           if (nbands > 6) then bp4[7] = SRFI7[lin,col];
           if (nbands > 7) then bp4[8] = SRFI8[lin,col];
           if (nbands > 8) then bp4[9] = SRFI9[lin,col];
        end
        if (i==5) then begin
           bp5[1] = SRFI1[lin,col]; bp5[2] = SRFI2[lin,col];
           bp5[3] = SRFI3[lin,col];
           if (nbands > 3) then bp5[4] = SRFI4[lin,col];
           if (nbands > 4) then bp5[5] = SRFI5[lin,col];
           if (nbands > 5) then bp5[6] = SRFI6[lin,col];
           if (nbands > 6) then bp5[7] = SRFI7[lin,col];
           if (nbands > 7) then bp5[8] = SRFI8[lin,col];
           if (nbands > 8) then bp5[9] = SRFI9[lin,col];
        end
        if (i==6) then begin
           bp6[1] = SRFI1[lin,col]; bp6[2] = SRFI2[lin,col];
           bp6[3] = SRFI3[lin,col];
           if (nbands > 3) then bp6[4] = SRFI4[lin,col];
           if (nbands > 4) then bp6[5] = SRFI5[lin,col];
           if (nbands > 5) then bp6[6] = SRFI6[lin,col];
           if (nbands > 6) then bp6[7] = SRFI7[lin,col];
           if (nbands > 7) then bp6[8] = SRFI8[lin,col];
           if (nbands > 8) then bp6[9] = SRFI9[lin,col];
        end
        if (i==7) then begin
           bp7[1] = SRFI1[lin,col]; bp7[2] = SRFI2[lin,col];
           bp7[3] = SRFI3[lin,col];
           if (nbands > 3) then bp7[4] = SRFI4[lin,col];
           if (nbands > 4) then bp7[5] = SRFI5[lin,col];
           if (nbands > 5) then bp7[6] = SRFI6[lin,col];
           if (nbands > 6) then bp7[7] = SRFI7[lin,col];
           if (nbands > 7) then bp7[8] = SRFI8[lin,col];
           if (nbands > 8) then bp7[9] = SRFI9[lin,col];
        end
        if (i==8) then begin
           bp8[1] = SRFI1[lin,col]; bp8[2] = SRFI2[lin,col];
           bp8[3] = SRFI3[lin,col];
           if (nbands > 3) then bp8[4] = SRFI4[lin,col];
           if (nbands > 4) then bp8[5] = SRFI5[lin,col];
           if (nbands > 5) then bp8[6] = SRFI6[lin,col];
           if (nbands > 6) then bp8[7] = SRFI7[lin,col];
           if (nbands > 7) then bp8[8] = SRFI8[lin,col];
           if (nbands > 8) then bp8[9] = SRFI9[lin,col];
        end
        if (i==9) then begin
           bp9[1] = SRFI1[lin,col]; bp9[2] = SRFI2[lin,col];
           bp9[3] = SRFI3[lin,col];
           if (nbands > 3) then bp9[4] = SRFI4[lin,col];
           if (nbands > 4) then bp9[5] = SRFI5[lin,col];
           if (nbands > 5) then bp9[6] = SRFI6[lin,col];
           if (nbands > 6) then bp9[7] = SRFI7[lin,col];
           if (nbands > 7) then bp9[8] = SRFI8[lin,col];
           if (nbands > 8) then bp9[9] = SRFI9[lin,col];
        end
     end
  end
# ------------------------------------------------------------
# WRTIE BIOPHYSICAL VECTORS TO CONSOLE WINDOW (SRFI Units):
  printf("BIOPHYSICAL CONTROL MATERIALS:\n");
  printf("   BPO: %s\n",bp0$);
  if (m > 1) then begin
     printf("   BP0:");
     nbandsm1 = nbands - 1;
     for i=1 to nbandsm1 begin
        printf("%8d",bp0[i]);
     end
     printf("%8d\n\n",bp0[nbands]);
  end
  printf("   BP1: %s\n",bp1$);
  if (m > 1) then begin
     printf("   BP1:");
     for i=1 to nbandsm1 begin
        printf("%8d",bp1[i]);
     end
     printf("%8d\n\n",bp1[nbands]);
  end
  printf("   BP2: %s\n",bp2$);
  if (m > 1) then begin
     printf("   BP2:");
     for i=1 to nbandsm1 begin
        printf("%8d",bp2[i]);
     end
     printf("%8d\n\n",bp2[nbands]);
  end
  for itc=3 to nTC begin
    if (itc==3) then begin
       printf("   BP3: %s\n",bp3$);
       if (m > 1) then begin
          printf("   BP3:");
          for i=1 to nbandsm1 begin
             printf("%8d",bp3[i]);
          end
          printf("%8d\n\n",bp3[nbands]);
       end
    end
    if (itc==4) then begin
       printf("   BP4: %s\n",bp4$);
       if (m > 1) then begin
          printf("   BP4:");
          for i=1 to nbandsm1 begin
             printf("%8d",bp4[i]);
          end
          printf("%8d\n\n",bp4[nbands]);
       end
    end
    if (itc==5) then begin
       printf("   BP5: %s\n",bp5$);
       if (m > 1) then begin
          printf("   BP5:");
          for i=1 to nbandsm1 begin
             printf("%8d",bp5[i]);
          end
          printf("%8d\n\n",bp5[nbands]);
       end
    end
    if (itc==6) then begin
       printf("   BP6: %s\n",bp6$);
       if (m > 1) then begin
          printf("   BP6:");
          for i=1 to nbandsm1 begin
             printf("%8d",bp6[i]);
          end
          printf("%8d\n\n",bp6[nbands]);
       end
    end
    if (itc==7) then begin
       printf("   BP7: %s\n",bp7$);
       if (m > 1) then begin
          printf("   BP7:");
          for i=1 to nbandsm1 begin
             printf("%8d",bp7[i]);
          end
          printf("%8d\n\n",bp7[nbands]);
       end
    end
    if (itc==8) then begin
       printf("   BP8: %s\n",bp8$);
       if (m > 1) then begin
          printf("   BP8:");
          for i=1 to nbandsm1 begin
             printf("%8d",bp8[i]);
          end
          printf("%8d\n\n",bp8[nbands]);
       end
    end
  end 
# ------------------------------------------------------------
# COMPUTE TC OFFSETS & ROTATION COEFFICIENTS:
  printf("TC OFFSETS (tc0) COEFFICIENTS (tc1 etc.):\n");
# tc0:
  printf("  tc0:  ");
  for i=1 to nbands begin
     if (m > 1) then begin
        tc0[i] = bp0[i];
     end
     printf("%8d",tc0[i]);
  end
# tc1:
  printf("\n  tc1:  ");
  if (m > 1) then begin
     sd2 = 0;
     for i=1 to nbands begin
        d = bp1[i] - tc0[i];
        sd2 = sd2 + d * d;
     end
     mag = sqrt(sd2);
  end
  for i=1 to nbands begin
     if (m > 1) then begin
        d = bp1[i] - tc0[i];
        tc1[i] = d / mag;
     end
     printf("%8.4f",tc1[i]);
  end
# tc2:
  printf("\n  tc2:  ");
     if (m > 1) then begin
     dot1 = 0;
     for i=1 to nbands begin
        d = bp2[i] - tc0[i];
        dot1 = dot1 + tc1[i] * d;
     end
     sf2 = 0;
     for i=1 to nbands begin
        d = bp2[i] - tc0[i];
        f = d - dot1 * tc1[i];
        sf2 = sf2 + f * f;
     end
     mag = sqrt(sf2);
  end
  for i=1 to nbands begin
     if (m > 1) then begin
        d = bp2[i] - tc0[i];
        f = d - dot1 * tc1[i];
        tc2[i] = f / mag;
     end
     printf("%8.4f",tc2[i]);
  end
# tc3:
  if (nTC > 2) then begin
     printf("\n  tc3:  ");
     if (m > 1) then begin
        dot1 = 0; dot2 = 0;
        for i=1 to nbands begin
           d = bp3[i] - tc0[i];
           dot1 = dot1 + tc1[i] * d;
           dot2 = dot2 + tc2[i] * d;
        end
        sf2 = 0;
        for i=1 to nbands begin
           d = bp3[i] - tc0[i];
           f = d - dot1 * tc1[i];
           f = f - dot2 * tc2[i];
           sf2 = sf2 + f * f;
        end
        mag = sqrt(sf2);
     end
     for i=1 to nbands begin
        if (m > 1) then begin
           d = bp3[i] - tc0[i];
           f = d - dot1 * tc1[i];
           f = f - dot2 * tc2[i];
           tc3[i] = f / mag;
        end
        printf("%8.4f",tc3[i]);
     end
  end
# tc4:
  if (nTC > 3) then begin
     printf("\n  tc4:  ");
     if (m > 1) then begin
        dot1 = 0; dot2 = 0; dot3 = 0;
        for i=1 to nbands begin
           d = bp4[i] - tc0[i];
           dot1 = dot1 + tc1[i] * d;
           dot2 = dot2 + tc2[i] * d;
           dot3 = dot3 + tc3[i] * d;
        end
        sf2 = 0;
        for i=1 to nbands begin
           d = bp4[i] - tc0[i];
           f = d - dot1 * tc1[i];
           f = f - dot2 * tc2[i];
           f = f - dot3 * tc3[i];
           sf2 = sf2 + f * f;
        end
        mag = sqrt(sf2);
     end
     for i=1 to nbands begin
        if (m > 1) then begin
           d = bp4[i] - tc0[i];
           f = d - dot1 * tc1[i];
           f = f - dot2 * tc2[i];
           f = f - dot3 * tc3[i];
           tc4[i] = f / mag;
        end
        printf("%8.4f",tc4[i]);
     end
  end
# tc5:
  if (nTC > 4) then begin
     printf("\n  tc5:  ");
     dot1 = 0; dot2 = 0; dot3 = 0; dot4 = 0;
     for i=1 to nbands begin
        d = bp5[i] - tc0[i];
        dot1 = dot1 + tc1[i] * d;
        dot2 = dot2 + tc2[i] * d;
        dot3 = dot3 + tc3[i] * d;
        dot4 = dot4 + tc4[i] * d;
     end
     sf2 = 0;
     for i=1 to nbands begin
        d = bp5[i] - tc0[i];
        f = d - dot1 * tc1[i];
        f = f - dot2 * tc2[i];
        f = f - dot3 * tc3[i];
        f = f - dot4 * tc4[i];
        sf2 = sf2 + f * f;
     end
     mag = sqrt(sf2);
     for i=1 to nbands begin
        d = bp5[i] - tc0[i];
        f = d - dot1 * tc1[i];
        f = f - dot2 * tc2[i];
        f = f - dot3 * tc3[i];
        f = f - dot4 * tc4[i];
        tc5[i] = f / mag;
        printf("%8.4f",tc5[i]);
     end
  end
# tc6:
  if (nTC > 5) then begin
     printf("\n  tc6:  ");
     dot1 = 0; dot2 = 0; dot3 = 0; dot4 = 0; dot5 = 0;
     for i=1 to nbands begin
        d = bp6[i] - tc0[i];
        dot1 = dot1 + tc1[i] * d;
        dot2 = dot2 + tc2[i] * d;
        dot3 = dot3 + tc3[i] * d;
        dot4 = dot4 + tc4[i] * d;
        dot5 = dot5 + tc5[i] * d;
     end
     sf2 = 0;
     for i=1 to nbands begin
        d = bp6[i] - tc0[i];
        f = d - dot1 * tc1[i];
        f = f - dot2 * tc2[i];
        f = f - dot3 * tc3[i];
        f = f - dot4 * tc4[i];
        f = f - dot5 * tc5[i];
        sf2 = sf2 + f * f;
     end
     mag = sqrt(sf2);
     for i=1 to nbands begin
        d = bp6[i] - tc0[i];
        f = d - dot1 * tc1[i];
        f = f - dot2 * tc2[i];
        f = f - dot3 * tc3[i];
        f = f - dot4 * tc4[i];
        f = f - dot5 * tc5[i];
        tc6[i] = f / mag;
        printf("%8.4f",tc6[i]);
     end
  end
# tc7:
  if (nTC > 6) then begin
     printf("\n  tc7:  ");
     dot1 = 0; dot2 = 0; dot3 = 0; dot4 = 0; dot5 = 0;
     dot6 = 0;
     for i=1 to nbands begin
        d = bp7[i] - tc0[i];
        dot1 = dot1 + tc1[i] * d;
        dot2 = dot2 + tc2[i] * d;
        dot3 = dot3 + tc3[i] * d;
        dot4 = dot4 + tc4[i] * d;
        dot5 = dot5 + tc5[i] * d;
        dot6 = dot6 + tc6[i] * d;
     end
     sf2 = 0;
     for i=1 to nbands begin
        d = bp7[i] - tc0[i];
        f = d - dot1 * tc1[i];
        f = f - dot2 * tc2[i];
        f = f - dot3 * tc3[i];
        f = f - dot4 * tc4[i];
        f = f - dot5 * tc5[i];
        f = f - dot6 * tc6[i];
        sf2 = sf2 + f * f;
     end
     mag = sqrt(sf2);
     for i=1 to nbands begin
        d = bp7[i] - tc0[i];
        f = d - dot1 * tc1[i];
        f = f - dot2 * tc2[i];
        f = f - dot3 * tc3[i];
        f = f - dot4 * tc4[i];
        f = f - dot5 * tc5[i];
        f = f - dot6 * tc6[i];
        tc7[i] = f / mag;
        printf("%8.4f",tc7[i]);
     end
  end
# tc8:
  if (nTC > 7) then begin
     printf("\n  tc8:  ");
     dot1 = 0; dot2 = 0; dot3 = 0; dot4 = 0; dot5 = 0;
     dot6 = 0; dot7 = 0;
     for i=1 to nbands begin
        d = bp8[i] - tc0[i];
        dot1 = dot1 + tc1[i] * d;
        dot2 = dot2 + tc2[i] * d;
        dot3 = dot3 + tc3[i] * d;
        dot4 = dot4 + tc4[i] * d;
        dot5 = dot5 + tc5[i] * d;
        dot6 = dot6 + tc6[i] * d;
        dot7 = dot7 + tc7[i] * d;
     end
     sf2 = 0;
     for i=1 to nbands begin
        d = bp8[i] - tc0[i];
        f = d - dot1 * tc1[i];
        f = f - dot2 * tc2[i];
        f = f - dot3 * tc3[i];
        f = f - dot4 * tc4[i];
        f = f - dot5 * tc5[i];
        f = f - dot6 * tc6[i];
        f = f - dot7 * tc7[i];
        sf2 = sf2 + f * f;
     end
     mag = sqrt(sf2);
     for i=1 to nbands begin
        d = bp8[i] - tc0[i];
        f = d - dot1 * tc1[i];
        f = f - dot2 * tc2[i];
        f = f - dot3 * tc3[i];
        f = f - dot4 * tc4[i];
        f = f - dot5 * tc5[i];
        f = f - dot6 * tc6[i];
        f = f - dot7 * tc7[i];
        tc8[i] = f / mag;
        printf("%8.4f",tc8[i]);
     end
  end
# tc9:
  if (nTC > 6) then begin
     printf("\n  tc9:  ");
     dot1 = 0; dot2 = 0; dot3 = 0; dot4 = 0; dot5 = 0;
     dot6 = 0; dot7 = 0; dot8 = 0;
     for i=1 to nbands begin
        d = bp9[i] - tc0[i];
        dot1 = dot1 + tc1[i] * d;
        dot2 = dot2 + tc2[i] * d;
        dot3 = dot3 + tc3[i] * d;
        dot4 = dot4 + tc4[i] * d;
        dot5 = dot5 + tc5[i] * d;
        dot6 = dot6 + tc6[i] * d;
        dot7 = dot7 + tc7[i] * d;
        dot8 = dot8 + tc8[i] * d;
     end
     sf2 = 0;
     for i=1 to nbands begin
        d = bp9[i] - tc0[i];
        f = d - dot1 * tc1[i];
        f = f - dot2 * tc2[i];
        f = f - dot3 * tc3[i];
        f = f - dot4 * tc4[i];
        f = f - dot5 * tc5[i];
        f = f - dot6 * tc6[i];
        f = f - dot7 * tc7[i];
        f = f - dot8 * tc8[i];
        sf2 = sf2 + f * f;
     end
     mag = sqrt(sf2);
     for i=1 to nbands begin
        d = bp9[i] - tc0[i];
        f = d - dot1 * tc1[i];
        f = f - dot2 * tc2[i];
        f = f - dot3 * tc3[i];
        f = f - dot4 * tc4[i];
        f = f - dot5 * tc5[i];
        f = f - dot6 * tc6[i];
        f = f - dot7 * tc7[i];
        f = f - dot8 * tc8[i];
        tc9[i] = f / mag;
        printf("%8.4f",tc9[i]);
     end
  end
# ------------------------------------------------------------
# SETUP OUTPUT RASTERS:
  raster DS0,TC1,DS1,TC2,DS2,TC3,DS3,TC4,DS4,TC5,DS5;
  raster TC6,DS6,TC7,DS7,TC8,DS8;
  GetOutputRaster(DS0,nlins,ncols,dtype$); SetNull(DS0,0);
  CopySubobjects(SRFI1,DS0,"GEOREF");
  GetOutputRaster(TC1,nlins,ncols,ttype$); SetNull(TC1,tcn);
  CopySubobjects(SRFI1,TC1,"GEOREF");
  GetOutputRaster(DS1,nlins,ncols,dtype$); SetNull(DS1,0);
  CopySubobjects(SRFI1,DS1,"GEOREF");
  GetOutputRaster(TC2,nlins,ncols,ttype$); SetNull(TC2,tcn);
  CopySubobjects(SRFI1,TC2,"GEOREF");
  GetOutputRaster(DS2,nlins,ncols,dtype$); SetNull(DS2,0);
  CopySubobjects(SRFI1,DS2,"GEOREF");
  if (nTC>2) then begin
     GetOutputRaster(TC3,nlins,ncols,ttype$);
     SetNull(TC3,tcn);
     CopySubobjects(SRFI1,TC3,"GEOREF");
     GetOutputRaster(DS3,nlins,ncols,dtype$); SetNull(DS3,0);
     CopySubobjects(SRFI1,DS3,"GEOREF");
  end
  if (nTC>3) then begin
     GetOutputRaster(TC4,nlins,ncols,ttype$);
     SetNull(TC4,tcn);
     CopySubobjects(SRFI1,TC4,"GEOREF");
     GetOutputRaster(DS4,nlins,ncols,dtype$); SetNull(DS4,0);
     CopySubobjects(SRFI1,DS4,"GEOREF");
  end
  if (nTC>4) then begin
     GetOutputRaster(TC5,nlins,ncols,ttype$);
     SetNull(TC5,tcn);
     CopySubobjects(SRFI1,TC5,"GEOREF");
     GetOutputRaster(DS5,nlins,ncols,dtype$); SetNull(DS5,0);
     CopySubobjects(SRFI1,DS5,"GEOREF");
  end
  if (nTC>5) then begin
     GetOutputRaster(TC6,nlins,ncols,ttype$);
     SetNull(TC6,tcn);
     CopySubobjects(SRFI1,TC6,"GEOREF");
     GetOutputRaster(DS6,nlins,ncols,dtype$); SetNull(DS6,0);
     CopySubobjects(SRFI1,DS6,"GEOREF");
  end
  if (nTC>6) then begin
     GetOutputRaster(TC7,nlins,ncols,ttype$);
     SetNull(TC7,tcn);
     CopySubobjects(SRFI1,TC7,"GEOREF");
     GetOutputRaster(DS7,nlins,ncols,dtype$); SetNull(DS7,0);
     CopySubobjects(SRFI1,DS7,"GEOREF");
  end
  if (nTC>7) then begin
     GetOutputRaster(TC8,nlins,ncols,ttype$);
     SetNull(TC8,tcn);
     CopySubobjects(SRFI1,TC8,"GEOREF");
     GetOutputRaster(DS8,nlins,ncols,dtype$); SetNull(DS8,0);
     CopySubobjects(SRFI1,DS8,"GEOREF");
  end
# ------------------------------------------------------------
# COMPUTE DS & TC VALUES & ASSIGN TO DS & TC RASTERS:
  printf("\n\nPRODUCING OUTPUT RASTERS:\n");
  printf("   DS0");
  for each SRFI1 begin
     isnull = IsNull(SRFI1);
     if (isnull <> 1) then begin
        srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
        if (nbands>3) then srfi[4] = SRFI4;
        if (nbands>4) then srfi[5] = SRFI5;
        if (nbands>5) then srfi[6] = SRFI6;
        if (nbands>6) then srfi[7] = SRFI7;
        if (nbands>7) then srfi[8] = SRFI8;
        if (nbands>8) then srfi[9] = SRFI9;
        sd2 = 0;
        for i=1 to nbands begin
           d = srfi[i] - tc0[i];
           sd2 = sd2 + d * d;
        end
        ds = round(sqrt(sd2));
        if (ds < 1) then ds = 1;
        DS0 = ds;
     end
  end
  printf(" TC1 DS1");
  for each SRFI1 begin
     isnull = IsNull(SRFI1);
     if (isnull <> 1) then begin
        srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
        if (nbands>3) then srfi[4] = SRFI4;
        if (nbands>4) then srfi[5] = SRFI5;
        if (nbands>5) then srfi[6] = SRFI6;
        if (nbands>6) then srfi[7] = SRFI7;
        if (nbands>7) then srfi[8] = SRFI8;
        if (nbands>8) then srfi[9] = SRFI9;
        ds = DS0;
        tc = 0;
        for i=1 to nbands begin
           tc = tc + tc1[i] * (srfi[i] - tc0[i]);
        end
        ds = ds * ds - tc * tc;
        if (ds < 1) then ds = 1;
        ds = round(sqrt(ds));
        if (tc < tcnp1) then tc = tcnp1;
        TC1 = tc; DS1 = ds;
     end
  end
  printf(" TC2 DS2");
  for each SRFI1 begin
     isnull = IsNull(SRFI1);
     if (isnull <> 1) then begin
        srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
        if (nbands>3) then srfi[4] = SRFI4;
        if (nbands>4) then srfi[5] = SRFI5;
        if (nbands>5) then srfi[6] = SRFI6;
        if (nbands>6) then srfi[7] = SRFI7;
        if (nbands>7) then srfi[8] = SRFI8;
        if (nbands>8) then srfi[9] = SRFI9;
        ds = DS1; tc = 0;
        for i=1 to nbands begin
           tc = tc + tc2[i] * (srfi[i] - tc0[i]);
        end
        ds = ds * ds - tc * tc;
        if (ds < 1) then ds = 1;
        ds = round(sqrt(ds));
        if (tc < tcnp1) then tc = tcnp1;
        TC2 = tc; DS2 = ds;
     end
  end
  if (nTC>2) then begin
     printf(" TC3 DS3");
     for each SRFI1 begin
        isnull = IsNull(SRFI1);
        if (isnull <> 1) then begin
           srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
           if (nbands>3) then srfi[4] = SRFI4;
           if (nbands>4) then srfi[5] = SRFI5;
           if (nbands>5) then srfi[6] = SRFI6;
           if (nbands>6) then srfi[7] = SRFI7;
           if (nbands>7) then srfi[8] = SRFI8;
           if (nbands>8) then srfi[9] = SRFI9;
           ds = DS2; tc = 0;
           for i=1 to nbands begin
              tc = tc + tc3[i] * (srfi[i] - tc0[i]);
           end
           ds = ds * ds - tc * tc;
           if (ds < 1) then ds = 1;
           ds = round(sqrt(ds));
           if (tc < tcnp1) then tc = tcnp1;
           TC3 = tc; DS3 = ds;
        end
     end
  end
  if (nTC>3) then begin
     printf(" TC4 DS4");
     for each SRFI1 begin
        isnull = IsNull(SRFI1);
        if (isnull <> 1) then begin
           srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
           if (nbands>3) then srfi[4] = SRFI4;
           if (nbands>4) then srfi[5] = SRFI5;
           if (nbands>5) then srfi[6] = SRFI6;
           if (nbands>6) then srfi[7] = SRFI7;
           if (nbands>7) then srfi[8] = SRFI8;
           if (nbands>8) then srfi[9] = SRFI9;
           ds = DS3; tc = 0;
           for i=1 to nbands begin
              tc = tc + tc4[i] * (srfi[i] - tc0[i]);
           end
           ds = ds * ds - tc * tc;
           if (ds < 1) then ds = 1;
           ds = round(sqrt(ds));
           if (tc < tcnp1) then tc = tcnp1;
           TC4 = tc; DS4 = ds;
        end
     end
  end
  if (nTC>4) then begin
     printf(" TC5 DS5");
     for each SRFI1 begin
        isnull = IsNull(SRFI1);
        if (isnull <> 1) then begin
           srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
           if (nbands>3) then srfi[4] = SRFI4;
           if (nbands>4) then srfi[5] = SRFI5;
           if (nbands>5) then srfi[6] = SRFI6;
           if (nbands>6) then srfi[7] = SRFI7;
           if (nbands>7) then srfi[8] = SRFI8;
           if (nbands>8) then srfi[9] = SRFI9;
           ds = DS4; tc = 0;
           for i=1 to nbands begin
              tc = tc + tc5[i] * (srfi[i] - tc0[i]);
           end
           ds = ds * ds - tc * tc;
           if (ds < 1) then ds = 1;
           ds = round(sqrt(ds));
           if (tc < tcnp1) then tc = tcnp1;
           TC5 = tc; DS5 = ds;
        end
     end
  end
  if (nTC>5) then begin
     printf(" TC6 DS6");
     for each SRFI1 begin
        isnull = IsNull(SRFI1);
        if (isnull <> 1) then begin
           srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
           if (nbands>3) then srfi[4] = SRFI4;
           if (nbands>4) then srfi[5] = SRFI5;
           if (nbands>5) then srfi[6] = SRFI6;
           if (nbands>6) then srfi[7] = SRFI7;
           if (nbands>7) then srfi[8] = SRFI8;
           if (nbands>8) then srfi[9] = SRFI9;
           ds = DS5; tc = 0;
           for i=1 to nbands begin
              tc = tc + tc6[i] * (srfi[i] - tc0[i]);
           end
           ds = ds * ds - tc * tc;
           if (ds < 1) then ds = 1;
           ds = round(sqrt(ds));
           if (tc < tcnp1) then tc = tcnp1;
           TC6 = tc; DS6 = ds;
        end
     end
  end
  if (nTC>6) then begin
     printf(" TC7 DS7");
     for each SRFI1 begin
        isnull = IsNull(SRFI1);
        if (isnull <> 1) then begin
           srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
           if (nbands>3) then srfi[4] = SRFI4;
           if (nbands>4) then srfi[5] = SRFI5;
           if (nbands>5) then srfi[6] = SRFI6;
           if (nbands>6) then srfi[7] = SRFI7;
           if (nbands>7) then srfi[8] = SRFI8;
           if (nbands>8) then srfi[9] = SRFI9;
           ds = DS6; tc = 0;
           for i=1 to nbands begin
              tc = tc + tc7[i] * (srfi[i] - tc0[i]);
           end
           ds = ds * ds - tc * tc;
           if (ds < 1) then ds = 1;
           ds = round(sqrt(ds));
           if (tc < tcnp1) then tc = tcnp1;
           TC7 = tc; DS7 = ds;
        end
     end
  end
  if (nTC>7) then begin
     printf(" TC8 DS8");
     for each SRFI1 begin
        isnull = IsNull(SRFI1);
        if (isnull <> 1) then begin
           srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
           if (nbands>3) then srfi[4] = SRFI4;
           if (nbands>4) then srfi[5] = SRFI5;
           if (nbands>5) then srfi[6] = SRFI6;
           if (nbands>6) then srfi[7] = SRFI7;
           if (nbands>7) then srfi[8] = SRFI8;
           if (nbands>8) then srfi[9] = SRFI9;
           ds = DS7; tc = 0;
           for i=1 to nbands begin
              tc = tc + tc8[i] * (srfi[i] - tc0[i]);
           end
           ds = ds * ds - tc * tc;
           if (ds < 1) then ds = 1;
           ds = round(sqrt(ds));
           if (tc < tcnp1) then tc = tcnp1;
           TC8 = tc; DS8 = ds;
        end
     end
  end
  printf("\n\n");
# ------------------------------------------------------------
# CLOSE INPUT RASTERS:
  CloseRaster(SRFI1); CloseRaster(SRFI2); CloseRaster(SRFI3);
  if (nbands>3) then CloseRaster(SRFI4);
  if (nbands>4) then CloseRaster(SRFI5);
  if (nbands>5) then CloseRaster(SRFI6);
  if (nbands>6) then CloseRaster(SRFI7);
  if (nbands>7) then CloseRaster(SRFI8);
  if (nbands>8) then CloseRaster(SRFI9);
# ------------------------------------------------------------
# MAKE SUBOBJECTS:
  printf("MAKING SUBOBJECTS:\n");
  printf("CONTRAST LOOK-UP TABLE PARAMETERS:\n");
  printf("RASTER      LOW      HIGH     POWER\n");
  CreateHistogram(DS0,0); CreatePyramid(DS0);
  hnmDS0 = HasNullMask(DS0); hnDS0 = HasNull(DS0);
  if (hnDS0 == 1) then nvDS0 = NullValue(DS0);
  c1DS0 = GlobalPtile(DS0,0.5,hnDS0,nvDS0);
  c3DS0 = GlobalPtile(DS0,50,hnDS0,nvDS0);
  c2DS0 = GlobalPtile(DS0,99.5,hnDS0,nvDS0);
  pDS0 = expPower(c1DS0,c3DS0,c2DS0);
  smlContrast.InputLowerLimit = c1DS0;
  smlContrast.InputUpperLimit = c2DS0;
  smlContrast.Power           =  pDS0;
  smlContrast.Compute(DS0,clut1$,clut2$,clut3$);
  printf("  DS0 %9d %9d %9.2f\n",c1DS0,c2DS0,pDS0);
  CloseRaster(DS0);
  CreateHistogram(TC1,0); CreatePyramid(TC1);
  hnmTC1 = HasNullMask(TC1); hnTC1 = HasNull(TC1);
  if (hnTC1 == 1) then nvTC1 = NullValue(TC1);
  c1TC1 = GlobalPtile(TC1,0.5,hnTC1,nvTC1);
  c3TC1 = GlobalPtile(TC1,50,hnTC1,nvTC1);
  c2TC1 = GlobalPtile(TC1,99.5,hnTC1,nvTC1);
  pTC1 = expPower(c1TC1,c3TC1,c2TC1);
  smlContrast.InputLowerLimit = c1TC1;
  smlContrast.InputUpperLimit = c2TC1;
  smlContrast.Power           =  pTC1;
  smlContrast.Compute(TC1,clut1$,clut2$,clut3$);
  printf("  TC1 %9d %9d %9.2f\n",c1TC1,c2TC1,pTC1);
  CloseRaster(TC1);
  CreateHistogram(DS1,0); CreatePyramid(DS1);
  hnmDS1 = HasNullMask(DS1); hnDS1 = HasNull(DS1);
  if (hnDS1 == 1) then nvDS1 = NullValue(DS1);
  c1DS1 = GlobalPtile(DS1,0.5,hnDS1,nvDS1);
  c3DS1 = GlobalPtile(DS1,50,hnDS1,nvDS1);
  c2DS1 = GlobalPtile(DS1,99.5,hnDS1,nvDS1);
  pDS1 = expPower(c1DS1,c3DS1,c2DS1);
  smlContrast.InputLowerLimit = c1DS1;
  smlContrast.InputUpperLimit = c2DS1;
  smlContrast.Power           =  pDS1;
  smlContrast.Compute(DS1,clut1$,clut2$,clut3$);
  printf("  DS1 %9d %9d %9.2f\n",c1DS1,c2DS1,pDS1);
  CloseRaster(DS1);
  CreateHistogram(TC2,0); CreatePyramid(TC2);
  hnmTC2 = HasNullMask(TC2); hnTC2 = HasNull(TC2);
  if (hnTC2 == 1) then nvTC2 = NullValue(TC2);
  c1TC2 = GlobalPtile(TC2,0.5,hnTC2,nvTC2);
  c3TC2 = GlobalPtile(TC2,50,hnTC2,nvTC2);
  c2TC2 = GlobalPtile(TC2,99.5,hnTC2,nvTC2);
  pTC2 = expPower(c1TC2,c3TC2,c2TC2);
  smlContrast.InputLowerLimit = c1TC2;
  smlContrast.InputUpperLimit = c2TC2;
  smlContrast.Power           =  pTC2;
  smlContrast.Compute(TC2,clut1$,clut2$,clut3$);
  printf("  TC2 %9d %9d %9.2f\n",c1TC2,c2TC2,pTC2);
  CloseRaster(TC2);
  CreateHistogram(DS2,0); CreatePyramid(DS2);
  hnmDS2 = HasNullMask(DS2); hnDS2 = HasNull(DS2);
  if (hnDS2 == 1) then nvDS2 = NullValue(DS2);
  c1DS2 = GlobalPtile(DS2,0.5,hnDS2,nvDS2);
  c3DS2 = GlobalPtile(DS2,50,hnDS2,nvDS2);
  c2DS2 = GlobalPtile(DS2,99.5,hnDS2,nvDS2);
  pDS2 = expPower(c1DS2,c3DS2,c2DS2);
  smlContrast.InputLowerLimit = c1DS2;
  smlContrast.InputUpperLimit = c2DS2;
  smlContrast.Power           =  pDS2;
  smlContrast.Compute(DS2,clut1$,clut2$,clut3$);
  printf("  DS2 %9d %9d %9.2f\n",c1DS2,c2DS2,pDS2);
  CloseRaster(DS2);
  if (nTC>2) then begin
     CreateHistogram(TC3,0); CreatePyramid(TC3);
     hnmTC3 = HasNullMask(TC3); hnTC3 = HasNull(TC3);
     if (hnTC3 == 1) then nvTC3 = NullValue(TC3);
     c1TC3 = GlobalPtile(TC3,0.5,hnTC3,nvTC3);
     c3TC3 = GlobalPtile(TC3,50,hnTC3,nvTC3);
     c2TC3 = GlobalPtile(TC3,99.5,hnTC3,nvTC3);
     pTC3 = expPower(c1TC3,c3TC3,c2TC3);
     smlContrast.InputLowerLimit = c1TC3;
     smlContrast.InputUpperLimit = c2TC3;
     smlContrast.Power           =  pTC3;
     smlContrast.Compute(TC3,clut1$,clut2$,clut3$);
     printf("  TC3 %9d %9d %9.2f\n",c1TC3,c2TC3,pTC3);
     CloseRaster(TC3);
     CreateHistogram(DS3,0); CreatePyramid(DS3);
     hnmDS3 = HasNullMask(DS3); hnDS3 = HasNull(DS3);
     if (hnDS3 == 1) then nvDS3 = NullValue(DS3);
     c1DS3 = GlobalPtile(DS3,0.5,hnDS3,nvDS3);
     c3DS3 = GlobalPtile(DS3,50,hnDS3,nvDS3);
     c2DS3 = GlobalPtile(DS3,99.5,hnDS3,nvDS3);
     pDS3 = expPower(c1DS3,c3DS3,c2DS3);
     smlContrast.InputLowerLimit = c1DS3;
     smlContrast.InputUpperLimit = c2DS3;
     smlContrast.Power           =  pDS3;
     smlContrast.Compute(DS3,clut1$,clut2$,clut3$);
     printf("  DS3 %9d %9d %9.2f\n",c1DS3,c2DS3,pDS3);
     CloseRaster(DS3);
  end
  if (nTC>3) then begin
     CreateHistogram(TC4,0); CreatePyramid(TC4);
     hnmTC4 = HasNullMask(TC4); hnTC4 = HasNull(TC4);
     if (hnTC4 == 1) then nvTC4 = NullValue(TC4);
     c1TC4 = GlobalPtile(TC4,0.5,hnTC4,nvTC4);
     c3TC4 = GlobalPtile(TC4,50,hnTC4,nvTC4);
     c2TC4 = GlobalPtile(TC4,99.5,hnTC4,nvTC4);
     pTC4 = expPower(c1TC4,c3TC4,c2TC4);
     smlContrast.InputLowerLimit = c1TC4;
     smlContrast.InputUpperLimit = c2TC4;
     smlContrast.Power           =  pTC4;
     smlContrast.Compute(TC4,clut1$,clut2$,clut3$);
     printf("  TC4 %9d %9d %9.2f\n",c1TC4,c2TC4,pTC4);
     CloseRaster(TC4);
     CreateHistogram(DS4,0); CreatePyramid(DS4);
     hnmDS4 = HasNullMask(DS4); hnDS4 = HasNull(DS4);
     if (hnDS4 == 1) then nvDS4 = NullValue(DS4);
     c1DS4 = GlobalPtile(DS4,0.5,hnDS4,nvDS4);
     c3DS4 = GlobalPtile(DS4,50,hnDS4,nvDS4);
     c2DS4 = GlobalPtile(DS4,99.5,hnDS4,nvDS4);
     pDS4 = expPower(c1DS4,c3DS4,c2DS4);
     smlContrast.InputLowerLimit = c1DS4;
     smlContrast.InputUpperLimit = c2DS4;
     smlContrast.Power           =  pDS4;
     smlContrast.Compute(DS4,clut1$,clut2$,clut3$);
     printf("  DS4 %9d %9d %9.2f\n",c1DS4,c2DS4,pDS4);
     CloseRaster(DS4);
  end
  if (nTC>4) then begin
     CreateHistogram(TC5,0); CreatePyramid(TC5);
     hnmTC5 = HasNullMask(TC5); hnTC5 = HasNull(TC5);
     if (hnTC5 == 1) then nvTC5 = NullValue(TC5);
     c1TC5 = GlobalPtile(TC5,0.5,hnTC5,nvTC5);
     c3TC5 = GlobalPtile(TC5,50,hnTC5,nvTC5);
     c2TC5 = GlobalPtile(TC5,99.5,hnTC5,nvTC5);
     pTC5 = expPower(c1TC5,c3TC5,c2TC5);
     smlContrast.InputLowerLimit = c1TC5;
     smlContrast.InputUpperLimit = c2TC5;
     smlContrast.Power           =  pTC5;
     smlContrast.Compute(TC5,clut1$,clut2$,clut3$);
     printf("  TC5 %9d %9d %9.2f\n",c1TC5,c2TC5,pTC5);
     CloseRaster(TC5);
     CreateHistogram(DS5,0); CreatePyramid(DS5);
     hnmDS5 = HasNullMask(DS5); hnDS5 = HasNull(DS5);
     if (hnDS5 == 1) then nvDS5 = NullValue(DS5);
     c1DS5 = GlobalPtile(DS5,0.5,hnDS5,nvDS5);
     c3DS5 = GlobalPtile(DS5,50,hnDS5,nvDS5);
     c2DS5 = GlobalPtile(DS5,99.5,hnDS5,nvDS5);
     pDS5 = expPower(c1DS5,c3DS5,c2DS5);
     smlContrast.InputLowerLimit = c1DS5;
     smlContrast.InputUpperLimit = c2DS5;
     smlContrast.Power           =  pDS5;
     smlContrast.Compute(DS5,clut1$,clut2$,clut3$);
     printf("  DS5 %9d %9d %9.2f\n",c1DS5,c2DS5,pDS5);
     CloseRaster(DS5);
  end
  if (nTC>5) then begin
     CreateHistogram(TC6,0); CreatePyramid(TC6);
     hnmTC6 = HasNullMask(TC6); hnTC6 = HasNull(TC6);
     if (hnTC6 == 1) then nvTC6 = NullValue(TC6);
     c1TC6 = GlobalPtile(TC6,0.5,hnTC6,nvTC6);
     c3TC6 = GlobalPtile(TC6,50,hnTC6,nvTC6);
     c2TC6 = GlobalPtile(TC6,99.5,hnTC6,nvTC6);
     pTC6 = expPower(c1TC6,c3TC6,c2TC6);
     smlContrast.InputLowerLimit = c1TC6;
     smlContrast.InputUpperLimit = c2TC6;
     smlContrast.Power           =  pTC6;
     smlContrast.Compute(TC6,clut1$,clut2$,clut3$);
     printf("  TC6 %9d %9d %9.2f\n",c1TC6,c2TC6,pTC6);
     CloseRaster(TC6);
     CreateHistogram(DS6,0); CreatePyramid(DS6);
     hnmDS6 = HasNullMask(DS6); hnDS6 = HasNull(DS6);
     if (hnDS6 == 1) then nvDS6 = NullValue(DS6);
     c1DS6 = GlobalPtile(DS6,0.5,hnDS6,nvDS6);
     c3DS6 = GlobalPtile(DS6,50,hnDS6,nvDS6);
     c2DS6 = GlobalPtile(DS6,99.5,hnDS6,nvDS6);
     pDS6 = expPower(c1DS6,c3DS6,c2DS6);
     smlContrast.InputLowerLimit = c1DS6;
     smlContrast.InputUpperLimit = c2DS6;
     smlContrast.Power           =  pDS6;
     smlContrast.Compute(DS6,clut1$,clut2$,clut3$);
     printf("  DS6 %9d %9d %9.2f\n",c1DS6,c2DS6,pDS6);
     CloseRaster(DS6);
  end
  if (nTC>6) then begin
     CreateHistogram(TC7,0); CreatePyramid(TC7);
     hnmTC7 = HasNullMask(TC7); hnTC7 = HasNull(TC7);
     if (hnTC7 == 1) then nvTC7 = NullValue(TC7);
     c1TC7 = GlobalPtile(TC7,0.5,hnTC7,nvTC7);
     c3TC7 = GlobalPtile(TC7,50,hnTC7,nvTC7);
     c2TC7 = GlobalPtile(TC7,99.5,hnTC7,nvTC7);
     pTC7 = expPower(c1TC7,c3TC7,c2TC7);
     smlContrast.InputLowerLimit = c1TC7;
     smlContrast.InputUpperLimit = c2TC7;
     smlContrast.Power           =  pTC7;
     smlContrast.Compute(TC7,clut1$,clut2$,clut3$);
     printf("  TC7 %9d %9d %9.2f\n",c1TC7,c2TC7,pTC7);
     CloseRaster(TC7);
     CreateHistogram(DS7,0); CreatePyramid(DS7);
     hnmDS7 = HasNullMask(DS7); hnDS7 = HasNull(DS7);
     if (hnDS7 == 1) then nvDS7 = NullValue(DS7);
     c1DS7 = GlobalPtile(DS7,0.5,hnDS7,nvDS7);
     c3DS7 = GlobalPtile(DS7,50,hnDS7,nvDS7);
     c2DS7 = GlobalPtile(DS7,99.5,hnDS7,nvDS7);
     pDS7 = expPower(c1DS7,c3DS7,c2DS7);
     smlContrast.InputLowerLimit = c1DS7;
     smlContrast.InputUpperLimit = c2DS7;
     smlContrast.Power           =  pDS7;
     smlContrast.Compute(DS7,clut1$,clut2$,clut3$);
     printf("  DS7 %9d %9d %9.2f\n",c1DS7,c2DS7,pDS7);
     CloseRaster(DS7);
  end
  if (nTC>7) then begin
     CreateHistogram(TC8,0); CreatePyramid(TC8);
     hnmTC8 = HasNullMask(TC8); hnTC8 = HasNull(TC8);
     if (hnTC8 == 1) then nvTC8 = NullValue(TC8);
     c1TC8 = GlobalPtile(TC8,0.5,hnTC8,nvTC8);
     c3TC8 = GlobalPtile(TC8,50,hnTC8,nvTC8);
     c2TC8 = GlobalPtile(TC8,99.5,hnTC8,nvTC8);
     pTC8 = expPower(c1TC8,c3TC8,c2TC8);
     smlContrast.InputLowerLimit = c1TC8;
     smlContrast.InputUpperLimit = c2TC8;
     smlContrast.Power           =  pTC8;
     smlContrast.Compute(TC8,clut1$,clut2$,clut3$);
     printf("  TC8 %9d %9d %9.2f\n",c1TC8,c2TC8,pTC8);
     CloseRaster(TC8);
     CreateHistogram(DS8,0); CreatePyramid(DS8);
     hnmDS8 = HasNullMask(DS8); hnDS8 = HasNull(DS8);
     if (hnDS8 == 1) then nvDS8 = NullValue(DS8);
     c1DS8 = GlobalPtile(DS8,0.5,hnDS8,nvDS8);
     c3DS8 = GlobalPtile(DS8,50,hnDS8,nvDS8);
     c2DS8 = GlobalPtile(DS8,99.5,hnDS8,nvDS8);
     pDS8 = expPower(c1DS8,c3DS8,c2DS8);
     smlContrast.InputLowerLimit = c1DS8;
     smlContrast.InputUpperLimit = c2DS8;
     smlContrast.Power           =  pDS8;
     smlContrast.Compute(DS8,clut1$,clut2$,clut3$);
     printf("  DS8 %9d %9d %9.2f\n",c1DS8,c2DS8,pDS8);
     CloseRaster(DS8);
  end
  printf("\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.");