'===========================================================================
' Subject: 3D CONTOUR PLOTTER                 Date: 09-23-99 (08:21)       
'  Author: Robert L. Roach                    Code: QB, QBasic, PDS        
'  Origin: rroach@cvd-pro.com               Packet: ALGOR.ABC
'===========================================================================
REM------------ PROGRAM CONTOUR ----------------
REM
REM     Linear interpolation version
REM
REM---------------------------------------------
        DEFINT I-N
        III = 81: JJJ = 81: NCM = 101
        DIM X(III, JJJ), Y(III, JJJ), Z(III, JJJ)
        DIM F1(III, JJJ), FP(NCM)
        DIM IT(III, JJJ)
        DIM XF(101), YF(101), ZF(101)
        DIM SUF$(3)
       
REM---------------------------------------------
        DG$ = CHR$(248)
        PH$ = CHR$(237)
        TH$ = CHR$(233)
        NCASE = 1

REM----------- SCREEN SELECT ---------------
        CLS
        LOCATE 2, 21: PRINT "CONTOUR PACKAGE - R.L. Roach"
        LOCATE 5, 26: PRINT "Which Monitor are you using?"
        LOCATE 7, 26: PRINT "  1 - VGA/Color (Default)"
        LOCATE 8, 26: PRINT "  2 - EGA/Color"
        LOCATE 9, 26: PRINT "  3 - Monochrome"
        LOCATE 10, 26: PRINT "  4 - Other"
        LOCATE 12, 26: INPUT "Enter Choice: ", IC
        IF IC = 0 THEN IC = 1
        MT = 9
        LOCATE 14, 18
        IF IC = 1 THEN MT = 12
        IF IC = 2 THEN MT = 9
        IF IC = 3 THEN MT = 10
        IF IC = 4 THEN INPUT "Enter screen # (if 9 or 10 don't work, try 2) ", MT

        FIRSTPL$ = "TRUE"
25      LOCATE 16, 22: PRINT "Which background color do you want?"
        LOCATE 17, 22: PRINT "  1 - Black (Default)"
        LOCATE 18, 22: PRINT "  2 - White (Use this one for hard copies)"
        LOCATE 19, 22: INPUT "Enter choice: ", IBC
        IF IBC = 0 THEN IBC = 1
        IF IBC = 1 THEN
          BCKGRCLR = 0
          CLR1 = 1
          CLR2 = 11
          CLR3 = 5
'          CLR1 = 11
'          CLR2 = 1
'          CLR3 = 1
REM          CLR4 = 14
          CLR4 = 10
        ELSEIF IBC = 2 THEN
          BCKGRCLR = 16
          CLR1 = 0
          CLR2 = 0
          CLR3 = 0
          CLR4 = 0
        ELSE
          FOR ROW = 16 TO 19
            LOCATE ROW, 22: PRINT STRING$(50, 32)
          NEXT
          GOTO 25
        END IF

        SCREEN MT
REM        COLOR 15

30 REM------------- USEFUL STUFF ---------------
        PI = 3.141592654#
       
        SUF$(1) = "st"
        SUF$(2) = "nd"
        SUF$(3) = "rd"

        ISS = 0
        XF = 1
        YF = 1
        X0 = 360
        Y0 = 150
        
        PH$ = CHR$(237)
        TH$ = CHR$(233)

REM------------ DEFAULT VIEW ANGLE -------------
        PHI = 15                'Angle around z-axis from x-axis
        THETA = 25              'Angle up from x-y plane (elevation)
        RADIUS = 200
     
REM------------ OTHER DEFAULTS --------------
        NCD = 35               ' Default number of contour lines
        KK = NCD

REM--------------- DEFINE F1(I,J) -------------------
        CLS
        FM = -1000
        FL = 1000
40      LOCATE 1, 15: PRINT "Where is F(I,J)?"
        LOCATE 2, 15: PRINT "  1 - Internal function"
        LOCATE 3, 15: PRINT "  2 - Data file (imax < 201, jmax < 61)"
        LOCATE 4, 15: INPUT "Enter choice: ", ICF
        FOR ROW = 1 TO 4
          LOCATE ROW, 15
          PRINT STRING$(45, 32)
        NEXT
        IF ICF = 1 GOTO 50
        IF ICF = 2 GOTO 70
        GOTO 40

50 REM----------- INTERNAL FUNCTION ---------------
        PHI = -125              'Angle around z-axis from x-axis
        THETA = 15              'Angle up from x-y plane (elevation)
       
        FIRST$ = "TRUE"
        II = 61
        JJ = 61
        LOCATE 16, 10: PRINT "Generating nifty 3D internal function..."
        LOCATE 18, 12: PRINT USING "JJ = ##"; JJ
        DTH = 2 * PI / (3 * (JJ - 1))
        TH = PI / 6 - DTH
        FOR J = 1 TO JJ
          LOCATE 19, 12: PRINT USING " J = ##"; J
          TH = TH + DTH
          ST = SIN(TH)
          CT = COS(TH)
          FJ = (J - 1) / (JJ - 1)
          FOR I = 1 TO II
            FI = (I - 1) / (II - 1)
            XX = -3 + 6 * (FI + .1 * FJ)
            YY = -3 + 6 * (FJ - .1 * FI)
            E1 = 4 * EXP(-((XX - 1) ^ 2 + (YY - 2) ^ 2))
            E2 = 2 * EXP(-((XX + 1) ^ 2 + (YY + 2) ^ 2))
            E3 = -2 * EXP(-((XX + 1) ^ 2 + (YY - 2) ^ 2))
            E4 = -3 * EXP(-((XX - 1) ^ 2 + (YY + 2) ^ 2))
            F1(I, J) = E1 + E2 + E3 + E4
            IF F1(I, J) > FM THEN FM = F1(I, J)
            IF F1(I, J) < FL THEN FL = F1(I, J)
          NEXT
        NEXT
       
        FOR ROW = 16 TO 19
          LOCATE ROW, 10: PRINT STRING$(50, 32)
        NEXT

REM----------------- DEFINE COORDS ---------------
        XMX = -1000
        YMX = -1000
        ZMX = -1000
        XML = 1000
        YML = 1000
        ZML = 1000
      
        FOR J = 1 TO JJ
          FJ = (J - 1) / (JJ - 1)
          FOR I = 1 TO II
            FI = (I - 1) / (II - 1)
            X(I, J) = .5 * (1 - 2 * FI)
            Y(I, J) = .5 * (1 - 2 * FJ)
            Z(I, J) = EXP(-5 * (FI ^ 2 + FJ ^ 2))
REM            Z(I, J) = 0
            IF X(I, J) > XMX THEN XMX = X(I, J)
            IF X(I, J) < XML THEN XML = X(I, J)
            IF Y(I, J) > YMX THEN YMX = Y(I, J)
            IF Y(I, J) < YML THEN YML = Y(I, J)
            IF Z(I, J) > ZMX THEN ZMX = Z(I, J)
            IF Z(I, J) < ZML THEN ZML = Z(I, J)
          NEXT
        NEXT
       
        IF (FIRSTONE$ = "TRUE") THEN
          JMIN = 1
          JMAX = JJ
          IMIN = 1
          IMAX = II
          FIRSTONE$ = "FALSE"
        END IF
      
        GOTO 180
       
70 REM------------ READ F1(I,J) FROM FILE -------------
        IF NCASE > 1 THEN
          IF UCASE$(YNSF$) = "Y" GOTO 180
          GOTO 71
        END IF

        FIRST$ = "TRUE"         'New data - used to get function min & max

REM--------- READ RESTART FOR TIME STEP INFO
REM        OPEN "RESTART" FOR INPUT AS #1
REM        INPUT #1, ISTP, TT
REM        CLOSE #1

REM--------- NOW GET THE DATA
71      LOCATE 1, 15: PRINT "Is your data:"
        LOCATE 2, 15: PRINT "  1 - 2D (ie. just f(x,y))"
        LOCATE 3, 15: PRINT "  2 - 2D (ie. x, y, and f(x,y))"
        LOCATE 4, 15: PRINT "  3 - 3D (ie. x, y, z, and f(x,y,z))"
        LOCATE 5, 15: PRINT "  4 - Never mind (quit)"
        LOCATE 6, 15: INPUT "Enter choice: ", ICD
        FOR ROW = 1 TO 6
          LOCATE ROW, 15
          PRINT STRING$(42, 32)
        NEXT
72      IF ICD = 4 THEN
          GOTO 40
       
REM----------------------------------------------------
REM             Just F(i,j)
REM----------------------------------------------------
        ELSEIF ICD = 1 THEN
          PHI = 90
          THETA = 90
          RADIUS = 200
73        LOCATE 1, 5: PRINT STRING$(40, 32)
          LOCATE 1, 5: INPUT "Enter data file name: ", NAM$
          OPEN NAM$ FOR INPUT AS #1
          INPUT #1, II, JJ
          IF JJ > JJJ THEN JJ = JJJ
          LOCATE 1, 15: PRINT USING "  Grid is ### x ###..."; II; JJ
          FOR J = 1 TO JJ
            LOCATE 1, 48: PRINT USING "J = ###"; J
            FJ = (J - 1) / (JJ - 1)
            FOR I = 1 TO II
              INPUT #1, F1(I, J)
              FI = (I - 1) / (II - 1)
              X(I, J) = .5 * (1 - 2 * FI)
              Y(I, J) = .5 * (1 - 2 * FJ)
              IF F1(I, J) > FM THEN FM = F1(I, J)
              IF F1(I, J) < FL THEN FL = F1(I, J)
              IF X(I, J) > XMX THEN XMX = X(I, J)
              IF X(I, J) < XML THEN XML = X(I, J)
              IF Y(I, J) > YMX THEN YMX = Y(I, J)
              IF Y(I, J) < YML THEN YML = Y(I, J)
              Z(I, J) = 0
            NEXT
          NEXT
          CLOSE #1
          ZMX = .5 * (XMX + YMX)
          ZML = .5 * (XML + YML)
          IMIN = 1
          IMAX = II
          JMIN = 1
          JMAX = JJ

REM----------------------------------------------------
REM             X(i,j), Y(i,j), F(i,j)
REM----------------------------------------------------
        ELSEIF ICD = 2 THEN
          PHI = -90
          THETA = 90
          RADIUS = 200
'          ON ERROR GOTO 75
75        LOCATE 1, 5: PRINT STRING$(40, 32)
          LOCATE 1, 15: INPUT "  Enter data file name: ", NAM$
          LOCATE 1, 15: PRINT STRING$(45, 32)
          OPEN NAM$ FOR INPUT AS #1
         
          INPUT #1, II, JJ
          DELX = 1 / (II - 2)
          DELY = 1 / (JJ - 2)
'          IF JJ > JJJ THEN JJ = JJJ
          LOCATE 1, 15: PRINT USING "  Grid is ### x ###..."; II; JJ
          JZ = 0
          JZZ = 1
          FOR J = 1 TO JJ
            LOCATE 2, 18: PRINT USING "J = ###"; J
            FOR I = 1 TO II
              INPUT #1, X(I, J), Y(I, J), F1(I, J)
              Z(I, J) = 0
              IF NCASE > 1 AND YNSF$ = "Y" GOTO 80
              IF X(I, J) > XMX THEN XMX = X(I, J)
              IF X(I, J) < XML THEN XML = X(I, J)
              IF Y(I, J) > YMX THEN YMX = Y(I, J)
              IF Y(I, J) < YML THEN YML = Y(I, J)
              IF F1(I, J) > FM THEN
                FM = F1(I, J)
                IMXX = I
                JMXX = JZ
              END IF
              IF F1(I, J) < FL THEN FL = F1(I, J)
80          NEXT
          NEXT
          CLOSE #1
         
          FOR ROW = 1 TO 2
            LOCATE ROW, 15
            PRINT STRING$(20, 32)
          NEXT
          JMIN = 1
          JMAX = JJ
          IMIN = 1
          IMAX = II
          IF NCASE > 1 AND YNSF$ = "Y" GOTO 180
          ZMX = .5 * (XMX + YMX)
          ZML = .5 * (XML + YML)
       
REM----------------------------------------------------
REM          X(i,j), Y(i,j), Z(i,j), F(i,j)
REM----------------------------------------------------
        ELSEIF ICD = 3 THEN
          ON ERROR GOTO 82
82        LOCATE 1, 5: PRINT STRING$(40, 32)
          LOCATE 1, 5: INPUT "  Enter data file name: ", NAM$
          OPEN NAM$ FOR INPUT AS #1
          INPUT #1, II, JJ
          IF JJ > JJJ THEN JJ = JJJ
          LOCATE 1, 5: PRINT STRING$(40, 32)
          LOCATE 1, 15: PRINT USING "  Grid is ### x ###..."; II; JJ
          FOR J = 1 TO JJ
            LOCATE 1, 48: PRINT USING "J = ###"; J
            FOR I = 1 TO II
              INPUT #1, X(I, J), Y(I, J), Z(I, J), F1(I, J)
              IF X(I, J) > XMX THEN XMX = X(I, J)
              IF X(I, J) < XML THEN XML = X(I, J)
              IF Y(I, J) > YMX THEN YMX = Y(I, J)
              IF Y(I, J) < YML THEN YML = Y(I, J)
              IF Z(I, J) > ZMX THEN ZMX = Z(I, J)
              IF Z(I, J) < ZML THEN ZML = Z(I, J)
              IF F1(I, J) > FM THEN FM = F1(I, J)
              IF F1(I, J) < FL THEN FL = F1(I, J)
            NEXT
          NEXT
          CLOSE #1
          IMIN = 1
          IMAX = II
          JMIN = 1
          JMAX = JJ
        ELSE
          GOTO 70
        END IF

180 REM---------------- PLOT OPTIONS ------------------
REM------- VIEW ANGLE
        IF ICD = 3 OR ICF = 1 THEN
          LOCATE 1, 5: PRINT STRING$(74, 32)
          LOCATE 1, 5
          PRINT USING "Present view angle (&=####&  &=####&  r=#####):"; PH$, PHI, DG$, TH$, THETA, DG$, RADIUS
          LOCATE 2, 5: PRINT STRING$(74, 32)
          LOCATE 2, 5
          PRINT USING " where:  &=angle around vertical axis  &=azimuth angle"; PH$, TH$
          LOCATE 3, 5: PRINT STRING$(74, 32)
          LOCATE 3, 5
          INPUT "Do you want a new view angle (Y/N)"; YN$
          IF UCASE$(YN$) = "Y" THEN
            LOCATE 4, 5: PRINT STRING$(74, 32)
            LOCATE 4, 5
            INPUT "Enter new view angle (Phi,Theta,Radius): ", PHI, THETA, RADIUS
          END IF
        END IF

REM------- INDEX LIMITS
        IF (FIRSTPL$ = "TRUE") THEN
          JMIN = 1
          JMAX = JJ
          IMIN = 1
          IMAX = II
          FIRSTPL$ = "FALSE"
        ELSE
          IF (JMAX > JJ) THEN JMAX = JJ
          IF (IMAX > II) THEN IMAX = II
        END IF
       
        IF NCASE = 1 THEN CLS
       
REM--------- PLOT A SUBREGION
'        IF NCASE > 1 GOTO 182
        FOR ROW = 1 TO 6
          LOCATE ROW, 1
          PRINT STRING$(75, 32)
        NEXT
        LOCATE 1, 1
        PRINT USING "  II = ###    JJ = ###"; II; JJ
        PRINT USING " Current Imin, Imax = ###  ###"; IMIN; IMAX
        PRINT USING " Current Jmin, Jmax = ###  ###"; JMIN; JMAX
        INPUT " Do you want to plot a sub-region (y/n)"; YSR$
        IF YSR$ = "Y" OR YSR$ = "y" THEN
          INPUT "Enter minimum and maximum I indices:", IMIN, IMAX
          INPUT "Enter minimum and maximum J indices:", JMIN, JMAX
          FIRST$ = "FALSE"
          FMX = -1E+10
          FMN = 1E+10
          FOR J = JMIN TO JMAX
            FOR I = IMIN TO IMAX
              IF F1(I, J) > FMX THEN FMX = F1(I, J)
              IF F1(I, J) < FMN THEN FMN = F1(I, J)
            NEXT
          NEXT
          IF NCASE = 1 THEN CLS
        END IF
        IF FIRST$ = "TRUE" OR YNSF$ <> "Y" THEN
          FIRST$ = "FALSE"
          FMX = -1E+10
          FMN = 1E+10
          FOR J = JMIN TO JMAX
            FOR I = IMIN TO IMAX
              IF F1(I, J) > FMX THEN FMX = F1(I, J)
              IF F1(I, J) < FMN THEN FMN = F1(I, J)
            NEXT
          NEXT
          IF NCASE = 1 THEN CLS
        END IF
       
REM--------- NUMBER OF CONTOUR LEVELS
'        IF NCASE > 1 GOTO 182
        FOR ROW = 1 TO 6
          LOCATE ROW, 1
          PRINT STRING$(75, 32)
        NEXT
        LOCATE 1, 1
        PRINT USING "   Current Fmin = ###.#####  Fmax = ###.#####"; FL; FM
        PRINT USING "   Window  Fmin = #.####^^^^  Fmax = #.####^^^^"; FMN; FMX
        PRINT "  "
        PRINT USING "   Current number of contours = ### (max = 101)"; KK
        KKK = KK
        INPUT "   How many contours do you want (Default = current):"; KK
        IF KK = 0 THEN KK = KKK
       
REM--------- SPACING OF CONTOUR LEVELS
        TOL = 1E-10
        FOR ROW = 1 TO 5
          LOCATE ROW, 1
          PRINT STRING$(75, 32)
        NEXT
        LOCATE 1, 1
        PRINT "   What values do you wish to use for the contours:"
        PRINT "     1 - Equal increment between current min & max"
        PRINT "     2 - Equal increment between window min & max (Default)"
        PRINT "     3 - Range of values of F1(I,J)"
        PRINT "     4 - Specific values of F1(I,J)"
        INPUT "   Enter choice: ", ICV
       
REM------------- Equal increment between grid min & max
182     IF ICV = 0 THEN ICV = 2
        IF ICV = 1 THEN
           FOR K = 1 TO KK
             FP(K) = FL + K * (FM - FL) / (KK + 1)
             IF ABS(FP(K)) < TOL THEN FP(K) = 1E-10
           NEXT
REM------------- Equal increment between window min & max
         ELSEIF ICV = 2 THEN
           FM = FMX
           FL = FMN
           FOR K = 1 TO KK
             FP(K) = FL + K * (FM - FL) / (KK + 1)
             IF ABS(FP(K)) < TOL THEN FP(K) = 1E-10
           NEXT

REM------------- Equal increment between given min & max
         ELSEIF ICV = 3 THEN
           FOR ROW = 1 TO 6
             LOCATE ROW, 2: PRINT STRING$(65, 32)
           NEXT
           LOCATE 1, 1
           PRINT USING "   Current Fmin = #.####^^^^  Fmax = #.####^^^^"; FL; FM
           PRINT USING "   Window  Fmin = #.####^^^^  Fmax = #.####^^^^"; FMN; FMX
           PRINT "  "
           INPUT "  Enter Fmin, Fmax: ", FL, FM
           FOR K = 1 TO KK
             FP(K) = FL + K * (FM - FL) / (KK + 1)
             IF ABS(FP(K)) < TOL THEN FP(K) = .0001
           NEXT

REM------------- Specific contour values
         ELSE
           FOR ROW = 1 TO 6
             LOCATE ROW, 2: PRINT STRING$(65, 32)
           NEXT
           LOCATE 1, 5: PRINT USING "Function Max = #.###^^^^  Function Min = #.###^^^^"; FM, FL
           IF NCASE > 1 THEN
             LOCATE 2, 5: INPUT "Use same levels as before (Y/N)"; YN$
             IF UCASE$(YN$) = "Y" GOTO 181
           END IF
           FOR K = 1 TO KK
             IF K < 4 THEN
               SUFF$ = SUF$(K)
             ELSE
               SUFF$ = "th"
             END IF
             LOCATE 3, 1: PRINT STRING$(50, 32)
             LOCATE 3, 1
             PRINT USING "  Enter ###& contour value: "; K; SUFF$
             IROW = CSRLIN: ICOL = 29
             LOCATE IROW - 1, ICOL: INPUT " ", FP(K)
           NEXT
        END IF
        IF NCASE > 1 AND YNSF$ = "Y" GOTO 184
        FOR ROW = 1 TO 8
          LOCATE ROW, 1
          PRINT STRING$(60, 32)
        NEXT
       
181 REM--------------- GRID ON SCREEN ------------------
        IF NCASE > 1 AND YNSF$ = "Y" GOTO 184
        LOCATE 1, 1
        PRINT "  "
        INPUT "   Want the grid drawn (Y/N) (Default = Y)"; YDG$
        LOCATE 2, 1: PRINT STRING$(60, 32)

        IF YDG$ = "" THEN YDG$ = "Y"
        IF UCASE$(YDG$) = "N" THEN
          LOCATE 2, 1
          INPUT " Do you want any lines drawn (Y/N)"; YNDL$
        END IF
        IF NCASE = 1 THEN CLS
       
REM----------- USEFUL VIEW ANGLE STUFF -----------
REM        IF FIRST$ = "FALSE" THEN
REM          LOCATE 5, 5
REM          PRINT USING "Current viewpoint (& & r) = #### #### ####.##"; PH$, TH$, PHI, THETA, RADIUS
REM          LOCATE 6, 5
REM          INPUT "Want new ones (Y/N)"; YN$
REM          IF UCASE$(YN$) = "Y" THEN
REM            LOCATE 7, 5
REM            PRINT USING "Enter new & (deg) & (deg) and r:"; PH$, TH$
REM            LOCATE 7, 38
REM            INPUT " "; PHI, THETA, RADIUS
REM          END IF
REM          CLS
REM        END IF

184     RV = RADIUS

        PH = PI * PHI / 180
        TH = PI * THETA / 180
        CP = COS(PH)
        SP = SIN(PH)
        CT = COS(TH)
        ST = SIN(TH)
         
        XV = RV * CP * CT
        YV = RV * SP * CT
        ZV = RV * ST

REM---------- GET WINDOW -----------
        IF NCASE > 1 GOTO 250

        IWMIN = 1
        IWMAX = II
        JWMIN = 1
        JWMAX = JJ
       
'        IF IMIN > 1 OR IMAX < II OR JMIN > 1 OR JMAX < JJ THEN
          LOCATE 1, 20: PRINT "Getting window for the domain..."
          LOCATE 2, 20: PRINT "How do you want the window to be computed?"
          LOCATE 3, 20: PRINT "  1 - Use entire domain"
          LOCATE 4, 20: PRINT "  2 - Use subdomain (default)"
          LOCATE 5, 20: PRINT "  3 - Use current domain"
          LOCATE 6, 20: PRINT "  4 - User specified subdomain"
          LOCATE 7, 20: INPUT "Enter choice: ", ICWD
         
          IF ICWD = 3 GOTO 200
          IF ICWD = 0 THEN ICWD = 2

          IF ICWD = 2 THEN
            IWMIN = IMIN
            IWMAX = IMAX
            JWMIN = JMIN
            JWMAX = JMAX
            CLS
          END IF
'        END IF
       
        XWL = 10000
        YWL = 10000
        XWM = -10000
        YWM = -10000

        FOR J = JWMIN TO JWMAX
          FOR I = IWMIN TO IWMAX
            XP = X(I, J)
            YP = Y(I, J)
            ZP = Z(I, J)
            GOSUB 2000
            IF XX < XWL THEN XWL = XX
            IF XX > XWM THEN XWM = XX
            IF YY < YWL THEN YWL = YY
            IF YY > YWM THEN YWM = YY
          NEXT
        NEXT
      
REM--------- GRID EXTREMES
        XDL = XWL
        XDM = XWM
        YDL = YWL
        YDM = YWM

REM--------- WINDOW
        FMF = 1
        IF ICD = 3 THEN FMF = 1.3
        DXW = FMF * (XWM - XWL)
        DYW = FMF * (YWM - YWL)
        XWH = .5 * (XWM + XWL)
        YWH = .5 * (YWM + YWL)
        IF DYW > .75 * DXW THEN
          DXW = DYW / .75
        ELSE
          DYW = .75 * DXW
        END IF
       
        XWL = XWH - .6 * DXW
        XWM = XWH + .6 * DXW
        YWL = YWH - .8 * DYW
        YWM = YWH + .6 * DYW

        IF ICWD = 4 THEN
          LOCATE 8, 10
          PRINT USING "Present window: (#.####^^^^ #.####^^^^)-(#.####^^^^ #.####^^^^)"; XWL; YWL; XWM; YWM
          LOCATE 9, 10
          PRINT USING "Grid extremes: (#.####^^^^ #.####^^^^)-(#.####^^^^ #.####^^^^)"; XDL, YDL, XDM, YDM
         
          LOCATE 10, 10
          INPUT "Enter new window (Xmin,Ymin,Xmax,Ymax): ", XWL, YWL, XWM, YWM
          LOCATE 12, 10
          INPUT "Would you like the window to be unstretched (Y/N)"; YNS$
         
          IF UCASE$(YNS$) = "Y" THEN
            DXW = XWM - XWL
            DYW = YWM - YWL
            XWH = .5 * (XWM + XWL)
            YWH = .5 * (YWM + YWL)
            IF DYW > .75 * DXW THEN
              DXW = DYW / .75
              XWL = XWH - .5 * DXW
              XWM = XWH + .5 * DXW
            ELSE
              DYW = .75 * DXW
              YWL = YWH - .55 * DYW
              YWM = YWH + .5 * DYW
            END IF
            CLS
          END IF
        END IF

200 REM----------- WINDOW --------------
        WINDOW (XWL, YWL)-(XWM, YWM)
        LINE (XWL, YWL)-(XWM, YWM), BCKGRCLR, BF
        LOCATE 1, 5
        PRINT USING "Fmin = #.#####^^^^  Fmax = #.#####^^^^ (### ###)"; FL, FM, IMXX, JMXX
        IF ICF = 2 THEN
          LOCATE 2, 10
          PRINT USING "File = &"; NAM$
        END IF
REM        PALETTE 0, 3
       
250 REM--------------- DRAW GRID --------------------
        JSTEP = 1
        ISTEP = 1

        IF UCASE$(YDG$) = "N" THEN
          JSTEP = JMAX - 1
          ISTEP = IMAX - 1
        END IF

        IF UCASE$(YNDL$) = "Y" THEN
251       GOSUB 1800
          LOCATE 5, 60: INPUT "Done (Y/N)"; YN$
          IF UCASE$(YN$) = "N" GOTO 251
          LOCATE 5, 60: PRINT STRING$(18, 32)
          GOTO 490
        END IF
            
        FOR J = JMIN TO JMAX STEP JSTEP
          FOR I = IMIN TO IMAX
            XP = X(I, J)
            YP = Y(I, J)
            ZP = Z(I, J)
            GOSUB 2000
            IF I = IMIN THEN PSET (XX, YY)
            LINE -(XX, YY), CLR1
          NEXT
        NEXT
       
        FOR I = IMIN TO IMAX STEP ISTEP
          FOR J = JMIN TO JMAX
            XP = X(I, J)
            YP = Y(I, J)
            ZP = Z(I, J)
            GOSUB 2000
            IF J = JMIN THEN PSET (XX, YY)
            LINE -(XX, YY), CLR1
          NEXT
        NEXT

490 REM----------------------------------------------------
REM                 DRAW CONTOURS
REM--------------------------------------------------------
        NC = 0
        K = 0
540     K = K + 1
        IF K > KK GOTO 1620
        LOCATE 1, 62: PRINT USING "F = #.####^^^^"; FP(K)
        LOCATE 2, 65: PRINT USING "Contour ###"; K

REM---------- ZERO OUT EXCLUSION ARRAY
        FOR J = JMIN TO JMAX
          FOR I = IMIN TO IMAX
            IT(I, J) = 0
          NEXT
        NEXT
       
        II1 = IMIN
        START$ = "TRUE"
        NS = 0

640 REM---------------- FIND START OF FP CONTOUR -------------
        FOR I = II1 TO IMAX - 1
          FOR J = JMIN + 1 TO JMAX
          
REM------- Check vertical segments for value of F1(I,J)
REM            IF I < IAA AND J <= JAA GOTO 850
REM            IF I >= IBB AND J <= JBB GOTO 850
            IF IT(I, J - 1) = 1 GOTO 850                'Already been here?
            IT(I, J - 1) = 1                            'New one
            Q1 = FP(K) - F1(I, J - 1)
            Q2 = F1(I, J) - FP(K)
            Q3 = F1(I, J) - F1(I, J - 1)
            IF Q1 ^ 2 + Q2 ^ 2 > Q3 ^ 2 GOTO 830
            FR = Q1 / Q3
            NC = NC + 1
            I1 = I
            I2 = I + 1
            J1 = J - 1
            J2 = J
            II1 = I1
            XP = X(I, J - 1) + FR * (X(I, J) - X(I, J - 1))
            YP = Y(I, J - 1) + FR * (Y(I, J) - Y(I, J - 1))
            ZP = Z(I, J - 1) + FR * (Z(I, J) - Z(I, J - 1))
            GOSUB 2000
            PSET (XX, YY)

REM----- CHECK TO SEE IF PREVIOUS HORIZONTAL SEGMENT WAS A STARTING POINT
            IF I > IMIN THEN
              Q1 = FP(K) - F1(I - 1, J)
              Q2 = F1(I, J) - FP(K)
              Q3 = F1(I, J) - F1(I - 1, J)
              IF Q1 ^ 2 + Q2 ^ 2 < Q3 ^ 2 THEN
                FR = Q1 / Q3
                XP = X(I - 1, J) + FR * (X(I, J) - X(I - 1, J))
                YP = Y(I - 1, J) + FR * (Y(I, J) - Y(I - 1, J))
                ZP = Z(I - 1, J) + FR * (Z(I, J) - Z(I - 1, J))
                GOSUB 2000
                CLR = CLR2
                IF FP(K) < 0 THEN CLR = CLR3
                IF ABS(FP(K)) < .000001 THEN CLR = CLR4
                LINE -(XX, YY), CLR
              END IF
            END IF
REM            IF I = IBB AND J < JBB GOTO 850
REM------------ END OF CHECK
            IJ1 = J1
            IJ2 = J2
            I = I1
            J = J1
            JD = 1
            ILOC = 3
            GOTO 860

REM------- Check horizontal segments for value of F1(I,J)
830 REM        IF I = IBB GOTO 850
            Q1 = FP(K) - F1(I, J - 1)
            Q2 = F1(I + 1, J - 1) - FP(K)
            Q3 = F1(I + 1, J - 1) - F1(I, J - 1)
            IF Q1 ^ 2 + Q2 ^ 2 > Q3 ^ 2 GOTO 850
            FR = Q1 / Q3
            NC = NC + 1
            I1 = I
            I2 = I + 1
            J1 = J - 1
            J2 = J
            XP = X(I, J - 1) + FR * (X(I + 1, J - 1) - X(I, J - 1))
            YP = Y(I, J - 1) + FR * (Y(I + 1, J - 1) - Y(I, J - 1))
            ZP = Z(I, J - 1) + FR * (Z(I + 1, J - 1) - Z(I, J - 1))
            GOSUB 2000
            PSET (XX, YY)
            II1 = I1
            II2 = I2
            IJ1 = J1
            I = I1
            J = J1
            ISS = 0
            ILOC = 1
            JD = 1
            GOTO 1150
850       NEXT
        NEXT
        GOTO 540
860     ISS = 0
        ILOC = 3
        JD = 1
910     ISS = ISS + 1
        IF ISS > 2 GOTO 540

REM----------------- HORIZONTAL LINE SEGMENT
        Q1 = FP(K) - F1(I1, J)
        Q2 = F1(I2, J) - FP(K)
        Q3 = F1(I2, J) - F1(I1, J)
        IF Q1 ^ 2 + Q2 ^ 2 > Q3 ^ 2 GOTO 1150
        FR = Q1 / (Q3 + 1E-10)
        ISS = 0
        IT(I1, J) = 1
        XP = X(I1, J) + FR * (X(I2, J) - X(I1, J))
        YP = Y(I1, J) + FR * (Y(I2, J) - Y(I1, J))
        ZP = Z(I1, J) + FR * (Z(I2, J) - Z(I1, J))
        GOSUB 2000
        CLR = CLR2
        IF FP(K) < 0 THEN CLR = CLR3
        IF ABS(FP(K)) < .000001 THEN CLR = CLR4
        LINE -(XX, YY), CLR
REM        DY = YP - YL
        IF (ILOC = 3) THEN JD = -JD
        J1 = J1 + JD
        J2 = J2 + JD
        I = I - JD
        ILOC = 1
'        XL = XP
'        YL = YP
'        ZL = ZP
REM        IF (I2 <= IAA) AND (J1 < JAA) GOTO 640
REM        IF (I2 > IBB) AND (J1 < JBB) GOTO 640
        IF (J2 > JMAX) OR (J1 < JMIN) GOTO 640
        IF (I1 = II1) AND (J1 = IJ1) AND (I2 = II2) GOTO 640
1150    ISS = ISS + 1
        IF ISS > 2 GOTO 540

REM----------- DIAGONAL LINE SEGMENT ------------
        Q1 = FP(K) - F1(I1, J2)
        Q2 = F1(I2, J1) - FP(K)
        Q3 = F1(I2, J1) - F1(I1, J2)
        IF Q1 ^ 2 + Q2 ^ 2 > Q3 ^ 2 GOTO 1350
        FR = Q1 / (Q3 + 1E-20)
        ISS = 0
        ILOC = 2
        XP = X(I1, J2) + FR * (X(I2, J1) - X(I1, J2))
        YP = Y(I1, J2) + FR * (Y(I2, J1) - Y(I1, J2))
        ZP = Z(I1, J2) + FR * (Z(I2, J1) - Z(I1, J2))
        IF JD = 1 THEN J = J2: I = I2
        IF JD = -1 THEN J = J1: I = I1
1350    ISS = ISS + 1

REM------------- VERTICAL LINE SEGMENT ---------
        Q1 = FP(K) - F1(I, J1)
        Q2 = F1(I, J2) - FP(K)
        Q3 = F1(I, J2) - F1(I, J1)
        IF Q1 ^ 2 + Q2 ^ 2 > Q3 ^ 2 GOTO 910
        FR = Q1 / (Q3 + 1E-10)
        ISS = 0
        IT(I, J1) = 1
        XP = X(I, J1) + FR * (X(I, J2) - X(I, J1))
        YP = Y(I, J1) + FR * (Y(I, J2) - Y(I, J1))
        ZP = Z(I, J1) + FR * (Z(I, J2) - Z(I, J1))
        GOSUB 2000
        CLR = CLR2
        IF FP(K) < 0 THEN CLR = CLR3
        IF ABS(FP(K)) < .000001 THEN CLR = CLR4
        LINE -(XX, YY), CLR
REM        DX = XP - XL
        IF ILOC = 1 THEN JD = -JD
        I1 = I1 + JD
        I2 = I2 + JD
        J = J - JD
        ILOC = 3
REM        IF (I2 <= IAA) AND (J1 < JAA) GOTO 640
REM        IF (I2 > IBB) AND (J1 < JBB) GOTO 640
        IF (I2 > IMAX) OR (I1 < IMIN) GOTO 640
        IF (J < JMIN) OR (J > JMAX) GOTO 640
        IF (I1 = II1) AND (J2 = IJ2) GOTO 640       ' STOP IF COMPLETE A LOOP
'        XL = XP
'        YL = YP
'        ZL = ZP
        GOTO 910

1620 REM--------- HOLD IT ------------
        LOCATE 28, 77: PRINT ".."
        DO: LOOP WHILE INKEY$ = ""
        FIRST$ = "FALSE"
        LOCATE 28, 30: PRINT STRING$(40, 32)
        LOCATE 28, 30: INPUT "Go again (Y/N)"; YN$
        IF UCASE$(YN$) = "Y" THEN
          LOCATE 28, 30: INPUT "Same function (Y/N)"; YNSF$
          LOCATE 28, 25: INPUT "Overlay new contours on old (Y/N)"; YN$
          LOCATE 28, 25: PRINT STRING$(40, 32)
          IF UCASE$(YN$) = "Y" THEN
            NCASE = NCASE + 1
            GOTO 180
          END IF
          IF UCASE$(YNSF$) = "Y" GOTO 180
          NCASE = 1
          CLS
          GOTO 25
        END IF

REM--------------- SO LONG -------------
        END



1800 REM----------- DRAW LINE SUBROUTINE ----------------
        LOCATE 5, 60: PRINT "Draw line...  "
        LOCATE 6, 60: INPUT "Enter I1,I2: ", IL1, IL2
        LOCATE 7, 60: INPUT "Enter J1,J2: ", JL1, JL2
        IF IL1 < IMIN THEN IL1 = IMIN
        IF IL2 > IMAX THEN IL2 = IMAX
        IF JL1 < JMIN THEN JL1 = JMIN
        IF JL2 > JMAX THEN JL2 = JMAX

        LOCATE 8, 60: INPUT "Line color?: ", CLRL

        LOCATE 9, 60: PRINT "How shall I do it?"
        LOCATE 10, 60: PRINT " 1 - straight line"
        LOCATE 11, 60: PRINT " 2 - follow grid"
1850    LOCATE 12, 60: INPUT "Enter choice: ", IC

        IF IC = 1 THEN
          XP = X(IL1, JL1)
          YP = Y(IL1, JL1)
          ZP = Z(IL1, JL1)
          GOSUB 2000
          PSET (XX, YY)

          XP = X(IL2, JL2)
          YP = Y(IL2, JL2)
          ZP = Z(IL2, JL2)
          GOSUB 2000
          LINE -(XX, YY), CLRL
        ELSEIF IC = 2 THEN
          XP = X(IL1, JL1)
          YP = Y(IL1, JL1)
          ZP = Z(IL1, JL1)
          GOSUB 2000
          PSET (XX, YY)
          FOR J = JL1 TO JL2
            FOR I = IL1 TO IL2
              XP = X(I, J)
              YP = Y(I, J)
              ZP = Z(I, J)
              GOSUB 2000
              LINE -(XX, YY), CLRL
            NEXT
          NEXT
        ELSE
          LOCATE 12, 74: PRINT STRING$(5, 32)
          GOTO 1850
        END IF

        FOR ROW = 5 TO 12
          LOCATE ROW, 60
          PRINT STRING$(19, 32)
        NEXT

REM-----------------------------------------------
        RETURN




2000 REM---------- SPACECOORDS --------------------
REM
REM     This subroutine computes the 3d coords
REM
REM------------------------------------------------
            RP = (XP * XV + YP * YV + ZP * ZV) / RV
            XR = RP * CP * CT
            YR = RP * SP * CT
            ZR = RP * ST
            DX = (XP - XR) * (RV + RP) / RV
            DY = (YP - YR) * (RV + RP) / RV
            DZ = (ZP - ZR) * (RV + RP) / RV
            XX = -DX * SP + DY * CP
            YY = -(DX * CP + DY * SP) * ST + DZ * CT

            RETURN