'=========================================================================== ' 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