'ORB2.BAS
'second module of ORBWEAVE.BAS
'by Peter Bird, UCLA, August 1993-March 2000
'(c)Copyright 1993, 2000 by Peter Bird and the Regents of the University of California.
DECLARE FUNCTION ATAN2D# (Y#, X#)
DECLARE FUNCTION ATAN2F! (Y!, X!)
DECLARE FUNCTION GETBIT% (NUMBER%, PLACE%)
DECLARE FUNCTION Inside% (TempVec!(), NPlanes%, Normals!())
DECLARE FUNCTION Principal! (angle!)
DECLARE SUB ABG2LonLat (Ins, TempVec!(), Outs, Lon!, Lat!)
DECLARE SUB ABG2XY (Ins, TempVec!(), Outs, Visible%, X!, Y!)
DECLARE SUB Blanker (R1%, R2%, C1%, C2%)
DECLARE SUB Cross (Ins, V1!(), V2!(), Outs, VCross!())
DECLARE SUB Delay18th (n%)
DECLARE SUB DrawArc (Ins, k1%, k2%, NodeABG!(), NumSeg%)
DECLARE SUB DropNode (Ins, AKA%, n%, NFL%, NUMEL%, Mods, EQCM!(), NFAKEN%, NMemo%(), NODEF%(), NODES%(), NREALN%, NUMNOD%, NodeABG!())
DECLARE SUB FindFault (Ins, X!, Y!, NFL%, NODEF%(), NodeABG!(), Tolerance!, Outs, ie%)
DECLARE SUB Finish (FileNum%)
DECLARE SUB Flipped (Ins, i%, NODES%(), NodeABG!(), Outs, EMemo%())
DECLARE SUB Icosahedron (Ins, NSlice%)
DECLARE SUB IsJoined (Ins, i%, NFL%, NODEF%(), Outs, FMemo%())
DECLARE SUB NEXTto (Ins, i%, j%, NFL%, NODEF%(), NODES%(), NUMEL%, Outs, KFAULT%, jf%, KELE%, je%)
DECLARE SUB OnAnyE (Ins!, n%, n1%, NLast%, NODES%(), Outs!, ie%, je%)
DECLARE SUB OnAnyF (Ins!, n%, n1%, NLast%, NODEF%(), Outs!, ie%, je%)
DECLARE SUB Pixels (Ins!, X!, Y!, Outs!, col%, row%)
DECLARE SUB RaiseOne (Ins, VCol%(), VRow%())
DECLARE SUB Scaler ()
DECLARE SUB SETBIT (NUMBER%, PLACE%, ONOFF%)
DECLARE SUB Steradians (Ins, V1!(), V2!(), V3!(), Outs, SolidAngle!)
DECLARE SUB Unitize (Mods!, X#, Y#, Z#)
DECLARE SUB UnitVec (Mods, TempVec!())
DECLARE SUB WaitForKey (a$)
DECLARE SUB XandY (Ins!, col%, row%, Outs!, X!, Y!)
DECLARE SUB XY2ABG (Ins, X!, Y!, Outs, Outside%, TempVec!())
DECLARE SUB XY2LonLat (Ins, Which%, X!, Y!, Outs, Outside%, Lon!, Lat!)
'--------------------------------------------------------------
OPTION BASE 1
'--------------------------------------------------------------------
DIM SHARED Letters(2, 14) AS STRING * 1 'Letter associated with (col,line) in menu
'---------------------------------------------------------------------
'Following COMMON gets parameters and constant-length arrays
'(but not dynamic, variable-length arrays) to the subprograms
'DrawNode, DrawElement, DrawFault in the fastest way, without
'having to put them on the stack and take them off in every call:
DIM WinOut!(3), WinRight!(3), WinUp!(3)
DIM Cart2!(3, 3)
DIM FIPoint!(4)
DIM Icon1%(30), Icon2%(30)
DIM Scales!(2, 2)
DIM UnScale!(2, 2)
COMMON SHARED Cart2!(), Colored%, CursorOn%, DoIcon%, FEP%, FIPoint!(), HiColor%, HiCol%, HiRow%, Icon1%(), Icon2%(), LastPutC%, LastPutR%, R2C!, Scales!(), UnScale!(), WindowHeight!, WinOut!(), WinRight!(), WinUp!()
'---------------------------------------------------------------------

SUB ABG2LonLat (Ins, TempVec!(), Outs, Lon!, Lat!)
'converts any vectory in the (Alpha,Beta,Gamma) system
'   (where Alpha^2+Beta^2+Gamma^2 = 1), and
'   Alpha points toward Lat=0, Lon=0, and
'   Beta  points toward Lat=0, Lon=90E, and
'   Gamma points toward North Pole.
'to N longitude and E latitude, in radians.
Lon! = ATAN2F!(TempVec!(2), TempVec!(1))
EquaPart! = SQR(TempVec!(1) * TempVec!(1) + TempVec!(2) * TempVec!(2))
Lat! = ATAN2F!(TempVec!(3), EquaPart!)
END SUB

SUB ABG2XY (Ins, TempVec!(), Outs, Visible%, X!, Y!)
'Converts any unit vector (Alpha, Beta, Gamma),
'  where Alpha points to Lat=0, Lon=0,
'        Beta  points to Lat=0, Lon=90E,
'  and   Gamma points to North Pole,
'  to the (x,y) system of the current window
'  (which is independent of window rotation).
dot! = TempVec!(1) * WinOut!(1) + TempVec!(2) * WinOut!(2) + TempVec!(3) * WinOut!(3)
IF dot < 0! THEN
	Visible% = 0 'false
	X! = 0!
	Y! = 0!
ELSE
	Visible% = NOT 0 'true
	X! = TempVec!(1) * WinRight!(1) + TempVec!(2) * WinRight!(2) + TempVec!(3) * WinRight!(3)
	Y! = TempVec!(1) * WinUp!(1) + TempVec!(2) * WinUp!(2) + TempVec!(3) * WinUp!(3)
END IF
END SUB

SUB AddNode (Ins, MXNODE%, NFL%, NUMEL%, TempVec!(), Mods, EQCM!(), NMemo%(), NODEF%(), NODES%(), NREALN%, NUMNOD%, N1000%, NodeABG!(), Outs, n%)
'Add a real node at TempVec!()=(alpha,beta,gamma);
'bump all fake nodes up if necessary.
IF NUMNOD% = MXNODE% THEN
	LOCATE 1, 1: PRINT "CANNOT ADD ANY MORE NODES; LIMIT OF "; MXNODE%; " REACHED."
	EXIT SUB
END IF
NUMNOD% = NUMNOD% + 1
NREALN% = NREALN% + 1
IF NUMNOD% > NREALN% THEN 'move up all fake nodes
	FOR i% = NUMNOD% TO (NREALN% + 1) STEP -1
		FOR k% = 1 TO 4
			EQCM!(k%, i%) = EQCM!(k%, i% - 1)
		NEXT k%
		NMemo%(i%) = NMemo%(i% - 1)
		NodeABG(1, i%) = NodeABG!(1, i% - 1)
		NodeABG(2, i%) = NodeABG!(2, i% - 1)
		NodeABG(3, i%) = NodeABG!(3, i% - 1)
	NEXT i%
	FOR i% = 1 TO NUMEL%
		FOR j% = 1 TO 3
			IF NODES%(j%, i%) >= NREALN% THEN NODES%(j%, i%) = NODES%(j%, i%) + 1
		NEXT j%
	NEXT i%
	FOR i% = 1 TO NFL%
		FOR j% = 1 TO 4
			IF NODEF%(j%, i%) >= NREALN% THEN NODEF%(j%, i%) = NODEF%(j%, i%) + 1
		NEXT j%
	NEXT i%
END IF
IF NREALN% > N1000% THEN  'adjust round number N1000%
	IF (CLNG(N1000%) + 1000&) <= 32767& THEN
		N1000% = N1000% + 1000
	ELSE
		N1000% = (CLNG(N1000%) + 32767&) / 2
	END IF
END IF
n% = NREALN%
FOR k% = 1 TO 4
	EQCM!(k%, n%) = 0!
NEXT k%
NodeABG!(1, n%) = TempVec!(1)
NodeABG!(2, n%) = TempVec!(2)
NodeABG!(3, n%) = TempVec!(3)
NMemo%(n%) = 0
END SUB

FUNCTION ATAN2D# (Y#, X#)
'2-argument inverse tangent, in double precision.
'Returns answer between -Pi# and +Pi# radians.
IF (Y# = 0#) AND (X# = 0#) THEN
   'answer technically undefined; just return 0.
	ATAN2D# = 0#
ELSE
	IF X# > 0# THEN
	    'ordinary ATN function will work fine
		Argument# = Y# / X#
		ATAN2D# = ATN(Argument#)
	ELSEIF X# = 0# THEN
		IF Y# > 0# THEN
			ATAN2D# = 1.570796327#
		ELSE
			ATAN2D# = -1.570796327#
		END IF
	ELSE
	    'x < 0
	    Argument# = Y# / X#
	    Prince# = ATN(Argument#)
	    IF Prince# > 0# THEN
		    ATAN2D# = Prince# - 3.141592654#
	    ELSE
		    ATAN2D# = Prince# + 3.141592654#
	    END IF
	END IF
END IF
END FUNCTION

FUNCTION ATAN2F! (Y!, X!)
'arc-tangent function with two arguments, returning result from -Pi to +Pi.
IF (Y! = 0!) AND (X! = 0!) THEN
   'answer technically undefined; just return 0.
	ATAN2F! = 0!
ELSE
	IF X! > 0! THEN
	    'ordinary ATN function will work fine
		ATAN2F! = ATN(Y! / X!)
	ELSEIF X! = 0! THEN
		IF Y! > 0! THEN
			ATAN2F! = 1.570795
		ELSE
			ATAN2F! = -1.570795
		END IF
	ELSE
	    'x < 0
	    Prince! = ATN(Y! / X!)
	    IF Prince! > 0! THEN
		    ATAN2F! = Prince! - 3.141591
	    ELSE
		    ATAN2F! = Prince! + 3.141591
	    END IF
	END IF
END IF
END FUNCTION

SUB Blanker (R1%, R2%, C1%, C2%)
FOR i% = R1% TO R2%
	LOCATE i%, C1%
	FOR j% = C1% TO C2%
		PRINT " ";
	NEXT j%
NEXT i%
END SUB

SUB CIRCUIT (Ins, MXBN%, NFL%, NUMEL%, Mods, EQCM!(), NFAKEN%, NMemo%(), NODEF%(), NODES%(), NREALN%, NUMNOD%, NodeABG!(), Outs, FAILED%, NCOND%, NODCON%(), Work, FMemo%(), EMemo%())
'
'MAKE A LIST, "NODCON%()" OF THE "NCOND%" NODES THAT ARE ON THE BOUNDARY
'AND REQIRE BOUNDARY CONDITIONS; THESE ARE IN COUNTERCLOCKWISE ORDER.
'IN WORKSPACE, MAKES AND USES LOGICAL LISTS OF WHICH ELEMENT SIDES ARE
'EXTERNAL: "EMemo%(NUMEL%)" AND "FMemo%(NFL%)".
'
'Technical Note:
'   Outside this routine, the integer arrays NMemo%() and EMemo%()
'   have only the values 0 (False) or 1 (True, in an unconventional
'   form using only 1 bit; warning: this True differs from NOT False%!).
'   This routine will use only the higher (2, 4, and 8) bits,
'   leaving the 2^0 bit unchanged.
'   Integer array FMemo%() may have values 0, 1, 2, or 3 on entry;
'   only the higher (2^2 and 2^3) bits are used in this routine.
'   At line CleanUp below, the high bits are returned to 0.
'
DIM Per0!(3), TempVec!(3)
	  DropLoose% = 0  'False
	  CALL Blanker(16, 24, 2, 79)
	  LOCATE 17, 14
	  PRINT "Checking that all nodes belong to some element..."
	  i% = 1
LoopOni:
		  a$ = INKEY$: IF LEN(a$) > 0 THEN GOTO StopShort  'impatient?
		  LOCATE 17, 63: PRINT i%
		  CALL OnAnyE(Ins, i%, 1, NUMEL%, NODES%(), Outs, ie%, je%)
		  IF ie% = 0 THEN
			  CALL OnAnyF(Ins, i%, 1, NFL%, NODEF%(), Outs, ie%, je%)
			  IF ie% = 0 THEN
				  IF DropLoose% THEN
					CALL DropNode(Ins, 0, i%, NFL%, NUMEL%, Mods, EQCM!(), NFAKEN%, NMemo%(), NODEF%(), NODES%(), NREALN%, NUMNOD%, NodeABG!())
					i% = i% - 1
				  ELSE
					BEEP
					CALL Blanker(16, 24, 2, 79)
					LOCATE 17, 31
					PRINT "BAD GRID TOPOLOGY."
					LOCATE 19, 21
					PRINT "SOME NODES DO NOT BELONG TO ANY ELEMENT."
					LOCATE 21, 19
					PRINT "(You must either delete them or connect them.)"
					LOCATE 23, 22
					INPUT ; "DO YOU WISH TO DELETE ALL UNUSED NODES"; a$
					IF UCASE$(LEFT$(a$, 1)) = "Y" OR UCASE$(LEFT$(a$, 2)) = "OK" THEN
						DropLoose% = NOT 0
						CALL Blanker(16, 24, 2, 79)
						LOCATE 17, 14
						PRINT "Checking that all nodes belong to some element..."
						GOTO LoopOni
					ELSE
						CALL Blanker(16, 24, 2, 79)
						LOCATE 17, 19
						PRINT "PERIMETER TEST CANNOT BE COMPLETED."
						FAILED% = NOT 0
						GOTO Cleanup
					END IF
				  END IF
			  END IF
		  END IF
	  IF i% < NUMNOD% THEN i% = i% + 1: GOTO LoopOni
	  CALL Blanker(17, 17, 2, 79)
	  LOCATE 17, 19
	  PRINT "Making a list of all boundary nodes:"
	  NCOND% = 0   'prepare to count exterior nodes
	  FOR i% = 1 TO NUMEL%
		  FOR j% = 1 TO 3
			  a$ = INKEY$: IF LEN(a$) > 0 THEN GOTO StopShort
			  CALL NEXTto(Ins, i%, j%, NFL%, NODEF%(), NODES%(), NUMEL%, Outs, KFAULT%, jf9%, KELE%, je9%)
			  IF KELE% > 0 THEN
'                    (ORDINARY INTERIOR SIDE)
				  CALL SETBIT(EMemo%(i%), j%, 0)
			  ELSEIF KFAULT% = 0 THEN
'                    (EXTERIOR SIDE)
				  CALL SETBIT(EMemo%(i%), j%, 1)
				  n1% = NODES%((j% MOD 3) + 1, i%)
				  N2% = NODES%(((j% + 1) MOD 3) + 1, i%)
				  IF GETBIT%(NMemo%(n1%), 1) THEN 'I know about it.
				  ELSE
					  NCOND% = NCOND% + 1  'another exterior node
					  CALL SETBIT(NMemo%(n1%), 1, 1)
					  CALL SETBIT(NMemo%(n1%), 2, 1)
				  END IF
				  IF GETBIT%(NMemo%(N2%), 1) THEN
				  ELSE
					  NCOND% = NCOND% + 1
					  CALL SETBIT(NMemo%(N2%), 1, 1)
					  CALL SETBIT(NMemo%(N2%), 2, 1)
				  END IF
			  ELSE
'                    (TRIANGULAR ELEMENT HAS AN EXTERIOR FAULT ELEMENT
'                     ADJACENT TO IT)
				  CALL SETBIT(EMemo%(i%), j%, 0)
				  N2% = NODES%(((j% + 1) MOD 3) + 1, i%)
				  IF NODEF%(1, KFAULT%) = N2% THEN
					  CALL SETBIT(FMemo%(KFAULT%), 3, 1) 'side 2 external
					  FOR k% = 3 TO 4
						  n% = NODEF%(k%, KFAULT%)
						  IF GETBIT%(NMemo%(n%), 1) THEN 'already know this
						  ELSE
							  NCOND% = NCOND% + 1  'another external node
							  CALL SETBIT(NMemo%(n%), 1, 1)
							  CALL SETBIT(NMemo%(n%), 2, 1)
						  END IF
					  NEXT k%
				  ELSE
					  CALL SETBIT(FMemo%(KFAULT%), 2, 1) 'side 1 external
					  FOR k% = 1 TO 2
						  n% = NODEF%(k%, KFAULT%)
						  IF GETBIT%(NMemo%(n%), 1) THEN  'already know this
						  ELSE
							  NCOND% = NCOND% + 1  'another external node
							  CALL SETBIT(NMemo%(n%), 1, 1)
							  CALL SETBIT(NMemo%(n%), 2, 1)
						  END IF
					  NEXT k%
				  END IF
			  END IF
		  NEXT j%
		  LOCATE 18, 19
		  Ipc% = CINT(100! * CSNG(i%) / CSNG(NUMEL%))
		  PRINT USING "##### nodes found; ###% of search completed."; NCOND%; Ipc%
	  NEXT i%
	  IF NCOND% = 0 THEN
		  LOCATE 20, 19
		  PRINT "No boundary nodes found; grid covers entire sphere."
		  CALL Delay18th(54)
		  FAILED% = 0 'False
		  GOTO Cleanup
	  END IF
	  IF NCOND% > MXBN% THEN
		  CALL Blanker(16, 24, 2, 79)
		  LOCATE 16, 4
		  PRINT "NUMBER OF BOUNDARY NODES ("; NCOND%; ") EXCEEDS LIMIT OF "; MXBN%; ".";
		  LOCATE 18, 4
		  PRINT "TEST CANNOT BE COMPLETED."
		  BEEP
		  FAILED% = NOT 0  'True
		  GOTO Cleanup
	  ELSE
		  FAILED% = 0  'False
		  LOCATE 20, 15
		  PRINT "Linking these nodes in order to form the perimeter:"
	  END IF
	  IF NUMNOD% > NREALN% THEN
		  FOR i% = NREALN% + 1 TO NUMNOD%
			  IF GETBIT%(NMemo%(i%), 1) THEN
			  ELSE
				  CALL Blanker(16, 24, 2, 79)
				  LOCATE 16, 4
				  PRINT "BAD GRID TOPOLOGY; FAKE NODES ARE NOT PERMITTED IN THE INTERIOR."
				  LOCATE 18, 4
				  PRINT "Save, eXit, reNUMBER, restart DRAWGRID, and reLoad."
				  LOCATE 20, 4
				  PRINT "(When you reNUMBER, there is an option to convert fake nodes to real.)"
				  BEEP
				  FAILED% = NOT 0  'True
				  GOTO Cleanup
			  END IF
		  NEXT i%
	  END IF
'   BEGIN THE LOOP WITH LOWEST-NUMBERED BOUNDARY NODE
	  FOR i% = 1 TO NUMNOD%
		  IF GETBIT%(NMemo%(i%), 1) THEN GOTO L831
	  NEXT i%
L831:  NODCON%(1) = i%
	  Per0!(1) = NodeABG!(1, i%) 'coordinates at which search started
	  Per0!(2) = NodeABG!(2, i%)
	  Per0!(3) = NodeABG!(3, i%)
	  NDONE% = 1
	  NLEFT% = NCOND%
'      BEGINNING OF INDEFINATE LOOP WHICH TRACES AROUND THE PERIMETER.
'      EACH TIME, IT PROGRESSES BY ONE OF 3 STEPS:
'      -1 NODE AT A TIME ALONG A NEW TRIANGLE SIDE, OR
'      -1 NODE AT A TIME ALONG A NEW FAULT ELEMENT SIDE, OR
'      -BY FINDING ANOTHER NODE WHICH SHARES THE SAME LOCATION.
'
'      BEGINNING OF MAIN INDEFINATE LOOP:
L840:  node% = NODCON%(NDONE%)
		  a$ = INKEY$: IF LEN(a$) > 0 THEN GOTO StopShort  'Impatient?
		  'check that we have not already counted this node!
		  IF NDONE% > 1 THEN
			  FOR j% = 1 TO NDONE% - 1
				  IF node% = NODCON%(j%) THEN
					  CALL Blanker(16, 24, 2, 79)
					  LOCATE 17, 31
					  PRINT "BAD GRID TOPOLOGY."
					  LOCATE 18, 10
					  PRINT "CAME BACK TO INITIAL NODE BEFORE LINKING ALL BOUNDARY NODES;"
					  LOCATE 19, 26
					  PRINT "GRID IS NOT SIMPLY-CONNECTED;"
					  LOCATE 20, 3
					  PRINT "IT HAS EITHER SOME GAP(S) OR OVERLAP(S) INSIDE, OR SOME ISLAND(S) OUTSIDE."
					  BEEP
					  LOCATE 22, 7
					  PRINT "A special basemap file, BADNODES.DIG, is being created to show the"
					  LOCATE 23, 7
					  PRINT "locations of the 'extra' external nodes;"
					  LOCATE 24, 7
					  INPUT ; "enter desired symbol-size, in degrees of arc: ", BoxSize!
					  BoxHalf! = .5 * BoxSize!
					  CALL Blanker(24, 24, 7, 78)
					  LOCATE 23, 7
					  PRINT "locations of the 'extra' external nodes ..... (Working)"
					  OPEN "BADNODES.DIG" FOR OUTPUT AS #4
						  FOR i% = 1 TO NUMNOD%
							  IF GETBIT%(NMemo%(i%), 2) THEN  'it is external
								  node% = i%
								  OK% = 0  'false
								  FOR k% = 1 TO NDONE% - 1
									  IF node% = NODCON%(k%) THEN  'it is linked; OK
										  OK% = NOT 0 'true
									  END IF
								  NEXT k%
								  TempVec!(1) = NodeABG!(1, i%): TempVec!(2) = NodeABG!(2, i%): TempVec!(3) = NodeABG!(3, i%)
								  CALL ABG2LonLat(Ins, TempVec!(), Outs, Lon!, Lat!)
								  Lat! = Lat! * 57.2958
								  Lon! = Lon! * 57.2958
								  PRINT #4, USING "Boundary node ##### at (Lon +###.###_,Lat +##.###)"; node%; Lon!; Lat!
								  IF OK% THEN
									  PRINT #4, USING " +###.###^^^^_, +##.###^^^^"; Lon! - BoxHalf!; Lat! + BoxHalf!
									  PRINT #4, USING " +###.###^^^^_, +##.###^^^^"; Lon! - BoxHalf!; Lat! - BoxHalf!
									  PRINT #4, USING " +###.###^^^^_, +##.###^^^^"; Lon! + BoxHalf!; Lat! - BoxHalf!
									  PRINT #4, USING " +###.###^^^^_, +##.###^^^^"; Lon! + BoxHalf!; Lat! + BoxHalf!
									  PRINT #4, USING " +###.###^^^^_, +##.###^^^^"; Lon! - BoxHalf!; Lat! + BoxHalf!
									  PRINT #4, "*** end of line segment ***"
								  ELSE
									  PRINT #4, USING " +###.###^^^^_, +##.###^^^^"; Lon! - BoxHalf!; Lat!
									  PRINT #4, USING " +###.###^^^^_, +##.###^^^^"; Lon! + BoxHalf!; Lat!
									  PRINT #4, "*** end of line segment ***"
									  PRINT #4, USING "Boundary node #####_, continued"; node%
									  PRINT #4, USING " +###.###^^^^_, +##.###^^^^"; Lon!; Lat! + BoxHalf!
									  PRINT #4, USING " +###.###^^^^_, +##.###^^^^"; Lon!; Lat! - BoxHalf!
									  PRINT #4, "*** end of line segment ***"
								  END IF
							  END IF
						  NEXT i%
					  CLOSE #4
					  LOCATE 23, 7
					  PRINT "locations of the 'extra' external nodes ..... DONE.  Press Enter <Ù"
					  FAILED% = NOT 0  'True
					  GOTO Cleanup
				  END IF
			  NEXT j%
		  END IF
		  LOCATE 21, 15
		  Ipc% = CINT(100! * CSNG(NDONE%) / CSNG(NCOND%))
		  PRINT USING "##### nodes linked; ###% of perimeter completed."; NDONE%; Ipc%
		  Alpha! = NodeABG!(1, node%)
		  Beta! = NodeABG!(2, node%)
		  Gamma! = NodeABG!(3, node%)
'           LOOK FOR AN ADJACENT TRIANGULAR ELEMENT USING THIS NODE.
		  FOR i% = 1 TO NUMEL%
			  FOR j% = 1 TO 3
				  IF GETBIT%(EMemo%(i%), j%) THEN
					  n1% = NODES%((j% MOD 3) + 1, i%)
					  IF n1% = node% GOTO L846
				  END IF
			  NEXT j%
		  NEXT i%
		  GOTO L850
L846:       N2% = NODES%(((j% + 1) MOD 3) + 1, i%)
		  NDONE% = NDONE% + 1
		  IF NDONE% <= NCOND% THEN NODCON%(NDONE%) = N2%
		  CALL SETBIT(NMemo%(N2%), 1, 0)
		  NLEFT% = NLEFT% - 1
		  IF NLEFT% > 0 THEN
			  GOTO L840
		  ELSE
			  GOTO L870
		  END IF
'           ELSE, LOOK FOR AN ADJACENT FAULT ELEMENT USING THIS NODE.
L850:       FOR i% = 1 TO NFL%
			  IF GETBIT%(FMemo%(i%), 2) THEN 'side 1 is external
				  IF NODEF%(1, i%) = node% THEN
					  N2% = NODEF%(2, i%)
					  GOTO L856
				  END IF
			  ELSEIF GETBIT%(FMemo%(i%), 3) THEN 'side 2 is external
				  IF NODEF%(3, i%) = node% THEN
					  N2% = NODEF%(4, i%)
					  GOTO L856
				  END IF
			  END IF
		  NEXT i%
		  GOTO L860
L856:       NDONE% = NDONE% + 1
		  IF NDONE% <= NCOND% THEN NODCON%(NDONE%) = N2%
		  CALL SETBIT(NMemo%(N2%), 1, 0)
		  NLEFT% = NLEFT% - 1
		  IF NLEFT% > 0 THEN
			  GOTO L840
		  ELSE
			  GOTO L870
		  END IF
'           ELSE, LOOK FOR ANOTHER EXTERIOR NODE AT SAME LOCATION.
L860:       FOR i% = 1 TO NUMNOD%
			  IF GETBIT%(NMemo%(i%), 1) AND (i% <> node%) THEN 'different external
				  IF (NodeABG!(1, i%) = Alpha!) AND (NodeABG!(2, i%) = Beta!) AND (NodeABG!(3, i%) = Gamma!) THEN GOTO L867
			  END IF
		  NEXT i%
		  CALL Blanker(16, 24, 2, 79)
		  LOCATE 16, 10: PRINT "BAD GRID TOPOLOGY: WHILE TRACING PERIMETER,"
		  LOCATE 17, 10: PRINT "COULD NOT FIND ANY WAY TO CONTINUE,"
		  LOCATE 18, 10: PRINT "EITHER THROUGH SHARED BOUNDARY ELEMENTS, OR"
		  LOCATE 19, 10: PRINT "THROUGH OTHER BOUNDARY NODES SHARING THE SAME POSITION."
		  CALL ABG2LonLat(Ins, Per0!(), Outs, Lon!, Lat!)
		  Lon! = Lon! * 57.2958: Lat! = Lat! * 57.2958
		  LOCATE 21, 10: PRINT "Search began at (Lon"; Lon!; ", Lat"; Lat!; ") with node "; NODCON%(1);
		  LOCATE 22, 10: PRINT "circled around the perimeter in the counterclockwise direction,";
		  TempVec!(1) = Alpha!: TempVec!(2) = Beta!: TempVec!(3) = Gamma!
		  CALL ABG2LonLat(Ins, TempVec!(), Outs, Lon!, Lat!)
		  Lon! = Lon! * 57.2958: Lat! = Lat! * 57.2958
		  LOCATE 23, 10: PRINT "and failed at (Lon "; Lon!; ", Lat"; Lat!; ") with node"; node%; ".";
		  BEEP
		  FAILED% = NOT 0 'True
		  GOTO Cleanup
L867:       NDONE% = NDONE% + 1
		  IF NDONE% <= NCOND% THEN NODCON%(NDONE%) = i%
		  CALL SETBIT(NMemo%(i%), 1, 0)
		  NLEFT% = NLEFT% - 1
	  IF NLEFT% > 0 GOTO L840
'      END OF INDEFINATE LOOP WHICH TRACES AROUND PERIMETER.
L870:  LOCATE 21, 15
	  PRINT USING "##### nodes linked; 100% of perimeter completed."; NCOND%
'-----------------------------------------------------------------
Cleanup:
		'restore work arrays to their original bit values (0 or 1):
		FOR i% = 1 TO NUMNOD%
			NMemo%(i%) = NMemo%(i%) MOD 2
		NEXT i%
		FOR i% = 1 TO NUMEL%
			EMemo%(i%) = EMemo%(i%) MOD 2
		NEXT i%
		FOR i% = 1 TO NFL%
			FMemo%(i%) = FMemo%(i%) MOD 4
		NEXT i%
		EXIT SUB
'-----------------------------------------
StopShort:
FAILED% = NOT 0 'True
CALL Blanker(16, 24, 2, 79)
LOCATE 18, 21
PRINT "PERIMETER TEST WAS INTERRUPTED BY USER."
LOCATE 20, 24
PRINT "Press any key to return to menu."
GOTO Cleanup
'------------------------------------------
END SUB

SUB Cross (Ins, V1!(), V2!(), Outs, VCross!())
'vector cross product: V1 x V2 = VCross
VCross!(1) = V1!(2) * V2!(3) - V1!(3) * V2!(2)
VCross!(2) = V1!(3) * V2!(1) - V1!(1) * V2!(3)
VCross!(3) = V1!(1) * V2!(2) - V1!(2) * V2!(1)
END SUB

SUB Delay18th (n%)
'delay for n%/18 seconds
Start! = TIMER
StopAt! = Start! + CSNG(n%) / 18!
TimeNow! = TIMER
DO WHILE TimeNow! < StopAt!
	TimeNow! = TIMER
LOOP
EXIT SUB
END SUB

SUB DrawBase (UseColor%, ColOldBase%, RowOldBase%, LastEnd%)
'Add a line to the basemap, assuming file open and initialized.
DIM TempVec!(3)
GET #3, , TempVec!(1)
GET #3, , TempVec!(2)
GET #3, , TempVec!(3)
IF LastEnd% THEN
	IF TempVec!(1) > 1.1 THEN
		'consecutive end marks; do nothing
		LastEnd% = NOT 0 'true
	ELSE  'start new segment with single point, if visible
		CALL ABG2XY(Ins, TempVec!(), Outs, Visible%, XBASE!, YBASE!)
		IF Visible% THEN
			LastEnd% = 0     'false
			CALL Pixels(Ins, XBASE!, YBASE!, Outs, ColBase%, RowBase%)
			IF Colored% THEN
				PSET (ColBase%, RowBase%), UseColor%
			ELSE
				PSET (ColBase%, RowBase%)
			END IF
		ELSE
			LastEnd% = NOT 0 'true
		END IF
	END IF
ELSE
	IF TempVec!(1) > 1.1 THEN
		'segment is over; do nothing
		LastEnd% = NOT 0   'true
	ELSE
		CALL ABG2XY(Ins, TempVec!(), Outs, Visible%, XBASE!, YBASE!)
		IF Visible% THEN
			LastEnd% = 0      'false
			CALL Pixels(Ins, XBASE!, YBASE!, Outs, ColBase%, RowBase%)
			IF Colored% THEN
				LINE (ColOldBase%, RowOldBase%)-(ColBase%, RowBase%), UseColor%
			ELSE
				LINE (ColOldBase%, RowOldBase%)-(ColBase%, RowBase%)
			END IF
		ELSE
			LastEnd% = NOT 0  'true
		END IF
	END IF
END IF
ColOldBase% = ColBase%
RowOldBase% = RowBase%
END SUB

SUB DrawFault (FDIP!(), i%, NODEF%(), UseColor%, NodeABG!())
DIM ECol%(2), ERow%(2), TempVec!(3), Xt!(2), Yt!(2)
k1% = NODEF%(1, i%)
TempVec!(1) = NodeABG!(1, k1%): TempVec!(2) = NodeABG!(2, k1%): TempVec!(3) = NodeABG!(3, k1%)
CALL ABG2XY(Ins, TempVec!(), Outs, Visible%, Xt!(1), Yt!(1))
IF Visible% THEN  ELSE EXIT SUB
R21! = Xt!(1) * Xt!(1) + Yt!(1) * Yt!(1)
k2% = NODEF%(2, i%)
TempVec!(1) = NodeABG!(1, k2%): TempVec!(2) = NodeABG!(2, k2%): TempVec!(3) = NodeABG!(3, k2%)
CALL ABG2XY(Ins, TempVec!(), Outs, Visible%, Xt!(2), Yt!(2))
IF Visible% THEN  ELSE EXIT SUB
R22! = Xt!(2) * Xt!(2) + Yt!(2) * Yt!(2)
IF (R21! > R2C!) AND (R23! > R2C!) THEN EXIT SUB
CALL Pixels(Ins, Xt!(1), Yt!(1), Outs, ECol%(1), ERow%(1))
CALL Pixels(Ins, Xt!(2), Yt!(2), Outs, ECol%(2), ERow%(2))
IF Colored% THEN
	LINE (ECol%(1), ERow%(1))-(ECol%(2), ERow%(2)), UseColor%
	IF ABS(ECol%(1) - ECol%(2)) > ABS(ERow%(1) - ERow%(2)) THEN
		LINE (ECol%(1), ERow%(1) + 1)-(ECol%(2), ERow%(2) + 1), UseColor%
		LINE (ECol%(1), ERow%(1) - 1)-(ECol%(2), ERow%(2) - 1), UseColor%
	ELSE
		LINE (ECol%(1) + 1, ERow%(1))-(ECol%(2) + 1, ERow%(2)), UseColor%
		LINE (ECol%(1) - 1, ERow%(1))-(ECol%(2) - 1, ERow%(2)), UseColor%
	END IF
ELSE
	IF UseColor% THEN
		LINE (ECol%(1), ERow%(1))-(ECol%(2), ERow%(2))
	ELSE
		LINE (ECol%(1), ERow%(1))-(ECol%(2), ERow%(2)), 0
	END IF
END IF
'plot 4 dip symbols, where they won't overlap the nodes
DX! = Xt!(2) - Xt!(1)
DY! = Yt!(2) - Yt!(1)
Size! = .1 * SQR(DX! * DX! + DY! * DY!)
FOR m% = 1 TO 4
	s! = FIPoint!(m%)
	F1! = 1! - s!
	F2! = s!
	Dip! = FDIP!(1, i%) * F1! + FDIP!(2, i%) * F2!
	ADip! = 90! - ABS(Dip! - 90!)
	IF ADip! < 76! THEN
		X! = F1! * Xt!(1) + F2! * Xt!(2)
		Y! = F1! * Yt!(1) + F2! * Yt!(2)
		DF1DS! = -1!
		DF2DS! = 1!
		DXDS! = DF1DS! * Xt!(1) + DF2DS! * Xt!(2)
		DYDS! = DF1DS! * Yt!(1) + DF2DS! * Yt!(2)
		Arg! = ATAN2F!(DYDS!, DXDS!)
		IF Dip! < 90! THEN
			Normal! = Arg! - 1.5714
		ELSE
			Normal! = Arg! + 1.5714
		END IF
		IF ADip! > 55! THEN 'steep fault; use simple tick
			CALL Pixels(Ins, X!, Y!, Outs, ColA%, RowA%)
			XP! = X! + Size! * COS(Normal!)
			YP! = Y! + Size! * SIN(Normal!)
			CALL Pixels(Ins, XP!, YP!, Outs, ColB%, RowB%)
			IF Colored% THEN
				LINE (ColA%, RowA%)-(ColB%, RowB%), UseColor%
			ELSE
				IF UseColor% THEN
					LINE (ColA%, RowA%)-(ColB%, RowB%)
				ELSE
					LINE (ColA%, RowA%)-(ColB%, RowB%), 0
				END IF
			END IF
		ELSEIF ADip! > 35! THEN 'intermediate dip; use open box
			XP! = X! + .5 * Size! * COS(Arg!)
			YP! = Y! + .5 * Size! * SIN(Arg!)
			CALL Pixels(Ins, XP!, YP!, Outs, ColA%, RowA%)
			XP! = XP! + Size! * COS(Normal!)
			YP! = YP! + Size! * SIN(Normal!)
			CALL Pixels(Ins, XP!, YP!, Outs, ColB%, RowB%)
			IF Colored% THEN
				LINE (ColA%, RowA%)-(ColB%, RowB%), UseColor%
			ELSE
				IF UseColor% THEN
					LINE (ColA%, RowA%)-(ColB%, RowB%)
				ELSE
					LINE (ColA%, RowA%)-(ColB%, RowB%), 0
				END IF
			END IF
			XP! = XP! - Size! * COS(Arg!)
			YP! = YP! - Size! * SIN(Arg!)
			CALL Pixels(Ins, XP!, YP!, Outs, ColA%, RowA%)
			IF Colored% THEN
				LINE -(ColA%, RowA%), UseColor%
			ELSE
				IF UseColor% THEN
					LINE -(ColA%, RowA%)
				ELSE
					LINE -(ColA%, RowA%), 0
				END IF
			END IF
			XP! = XP! - Size! * COS(Normal!)
			YP! = YP! - Size! * SIN(Normal!)
			CALL Pixels(Ins, XP!, YP!, Outs, ColA%, RowA%)
			IF Colored% THEN
				LINE -(ColA%, RowA%), UseColor%
			ELSE
				IF UseColor% THEN
					LINE -(ColA%, RowA%)
				ELSE
					LINE -(ColA%, RowA%), 0
				END IF
			END IF
		ELSE ' shallow dip; triangle (filled if <= 27.5 degrees)
			XP! = X! + .5 * Size! * COS(Arg!)
			YP! = Y! + .5 * Size! * SIN(Arg!)
			CALL Pixels(Ins, XP!, YP!, Outs, ColA%, RowA%)
			XP! = X! + Size! * COS(Normal!)
			YP! = Y! + Size! * SIN(Normal!)
			CALL Pixels(Ins, XP!, YP!, Outs, ColB%, RowB%)
			IF Colored% THEN
				LINE (ColA%, RowA%)-(ColB%, RowB%), UseColor%
			ELSE
				IF UseColor% THEN
					LINE (ColA%, RowA%)-(ColB%, RowB%)
				ELSE
					LINE (ColA%, RowA%)-(ColB%, RowB%), 0
				END IF
			END IF
			XP! = X! - .5 * Size! * COS(Arg!)
			YP! = Y! - .5 * Size! * SIN(Arg!)
			CALL Pixels(Ins, XP!, YP!, Outs, ColC%, RowC%)
			IF Colored% THEN
				LINE -(ColC%, RowC%), UseColor%
			ELSE
				IF UseColor% THEN
					LINE -(ColC%, RowC%)
				ELSE
					LINE -(ColC%, RowC%), 0
				END IF
			END IF
			IF ADip! <= 27.5 THEN ' subduction zone: paint tick
			    ColP% = (ColA% + ColB% + ColC%) / 3
			    RowP% = (RowA% + RowB% + RowC%) / 3
			    IF Colored% THEN
				    PAINT (ColP%, RowP%), UseColor%, UseColor%
			    ELSE
				    PAINT (ColP%, RowP%)
			    END IF
			END IF
		END IF
	END IF
NEXT m%
END SUB

SUB DrawNode (n%, UseColor%, NodeABG!())
DIM TempVec!(3)
TempVec!(1) = NodeABG!(1, n%): TempVec!(2) = NodeABG!(2, n%): TempVec!(3) = NodeABG!(3, n%)
CALL ABG2XY(Ins, TempVec!(), Outs, Visible%, X!, Y!)
IF Visible% THEN
	IF (X! * X! + Y! * Y!) < R2C! THEN
		CALL Pixels(Ins, X!, Y!, Outs, col%, row%)
		IF Colored% THEN
			CIRCLE (col%, row%), 3, UseColor%
		ELSE
			IF UseColor% THEN
				CIRCLE (col%, row%), 3
			ELSE
				CIRCLE (col%, row%), 3, 0
			END IF
		END IF
	END IF
END IF
END SUB

SUB DropNode (Ins, AKA%, n%, NFL%, NUMEL%, Mods, EQCM!(), NFAKEN%, NMemo%(), NODEF%(), NODES%(), NREALN%, NUMNOD%, NodeABG!())
'eliminate node n%; if it occurs in element/fault lists, replace with AKA%
'(AKA% is the OLD alias, before any renumbering.)
NUMNOD% = NUMNOD% - 1
IF n% <= NREALN% THEN
	NREALN% = NREALN% - 1
ELSEIF NFAKEN% > 0 THEN
	NFAKEN% = NFAKEN% - 1
ELSE
	BEEP
	EXIT SUB
END IF
FOR k% = n% TO NUMNOD%
	FOR i% = 1 TO 4
		EQCM!(i%, k%) = EQCM!(i%, k% + 1)
	NEXT i%
	NMemo%(k%) = NMemo%(k% + 1)
	NodeABG!(1, k%) = NodeABG!(1, k% + 1)
	NodeABG!(2, k%) = NodeABG!(2, k% + 1)
	NodeABG!(3, k%) = NodeABG!(3, k% + 1)
NEXT k%
FOR k% = 1 TO 3
	FOR j% = 1 TO NUMEL%
		IF NODES%(k%, j%) = n% THEN NODES%(k%, j%) = AKA%
		IF NODES%(k%, j%) > n% THEN NODES%(k%, j%) = NODES%(k%, j%) - 1
	NEXT j%
NEXT k%
FOR k% = 1 TO 4
	FOR j% = 1 TO NFL%
		IF NODEF%(k%, j%) = n% THEN NODEF%(k%, j%) = AKA%
		IF NODEF%(k%, j%) > n% THEN NODEF%(k%, j%) = NODEF%(k%, j%) - 1
	NEXT j%
NEXT k%
END SUB

SUB Exists (Ins, n1%, NUMNOD%, X!, Y!, NodeABG!(), Tolerance!, Outs, n%)
'Check whether a node already exists at location and return number (or 0)
'Testing begins with node N1, to allow for repeat calls, finding more than 1.
'Tolerance is in radians.
DIM TempVec!(3)
CALL XY2ABG(Ins, X!, Y!, Outs, Outside%, TempVec!())
IF Outside% THEN
	n% = 0
ELSE
	R2Min! = 3E+38
	FOR i% = n1% TO NUMNOD%
		DA! = TempVec!(1) - NodeABG!(1, i%)
		DB! = TempVec!(2) - NodeABG!(2, i%)
		DG! = TempVec!(3) - NodeABG!(3, i%)
		R2! = DA! * DA! + DB! * DB! + DG! * DG!
		IF R2! < R2Min! THEN
			R2Min! = R2!
			iClose% = i%
		END IF
	NEXT i%
	RMin! = SQR(R2Min!)
	IF RMin! <= Tolerance! THEN
		n% = iClose%
	ELSE
		n% = 0
	END IF
END IF
END SUB

SUB FindFault (Ins, X!, Y!, NFL%, NODEF%(), NodeABG!(), Tolerance!, Outs, ie%)
'Returns number of fault element whose midpoint is near projected point (x,y),
'if any.  Otherwise, returns 0.
DIM TempVec!(3)
R2Min! = 3.3E+38
FOR i% = 1 TO NFL%
	k1% = NODEF%(1, i%)
	TempVec!(1) = NodeABG!(1, k1%): TempVec!(2) = NodeABG!(2, k1%): TempVec!(3) = NodeABG!(3, k1%)
	CALL ABG2XY(Ins, TempVec!(), Outs, Visible%, X1!, Y1!)
	IF Visible% THEN
		k2% = NODEF%(2, i%)
		TempVec!(1) = NodeABG!(1, k2%): TempVec!(2) = NodeABG!(2, k2%): TempVec!(3) = NodeABG!(3, k2%)
		CALL ABG2XY(Ins, TempVec!(), Outs, Visible%, X2!, Y2!)
		IF Visible% THEN
			XC! = .5 * (X1! + X2!)
			YC! = .5 * (Y1! + Y2!)
			DX! = X! - XC!
			DY! = Y! - YC!
			R2! = DX! * DX! + DY! * DY!
			IF R2! < R2Min! THEN
				R2Min! = R2!: ie% = i%
			END IF
		END IF
	END IF
NEXT i%
IF R2Min! > Tolerance! * Tolerance! THEN ie% = 0
END SUB

SUB Finish (FileNum%)
'guarantees that a whole line has been read from file #FileNum%,
'and that there are no remaining blanks, comments, etc. to confuse
'the next INPUT#FileNum instruction.
n& = SEEK(FileNum%)         'find position in file
SEEK #FileNum%, n& - 1&     'back up one byte
a$ = INPUT$(1, #FileNum%)   'get another copy of last byte read
IF ASC(a$) <> 10 THEN LINE INPUT #FileNum%, Dummy$
END SUB

SUB FirstScreen (Ins, BestMode%, Blue%, Foreground%, Green%, HotColor%, Lines%, PixelRatio!, Red%, White%, Yellow%, Outs, PageNow%)
'prints introductory screen, in graphics mode, using
'  multiple colors if possible.
SCREEN BestMode%: CLS : PageNow% = 0
'------------------------------- text portion ----------------------
IF Colored% THEN COLOR HotColor%
LOCATE 1, 1
PRINT "                  OrbWeaver v1.H";
LOCATE 2, 1
PRINT "an interactive editor for finite element grids of";
LOCATE 3, 1
PRINT "3-node spherical triangular elements and 4-node";
LOCATE 4, 1
PRINT "great-circle faults on a sphere, as used in";
LOCATE 5, 1
PRINT "my programs Shells, Restore, and NeoKinema.";
LOCATE 6, 1
PRINT "Capacity: up to 7600 nodes, 11400 elements,";
LOCATE 7, 1
PRINT "   and 1900 faults (running in 590 K of RAM).";
LOCATE 8, 1
PRINT "by Peter Bird";
LOCATE 9, 1
PRINT "   Department of Earth and Space Sciences";
LOCATE 10, 1
PRINT "   University of California";
LOCATE 11, 1
PRINT "   Los Angeles, CA 90095-1567";
LOCATE 12, 1
PRINT "(c)Copyright 1993-2002 by Peter Bird and the";
LOCATE 13, 1
PRINT "   Regents of the University of California.";
LOCATE Lines% - 10, 1
PRINT "REQUIRES:  Microsoft-compatible mouse & driver program"
PRINT "           (memory-resident), a graphics card, and a graphics monitor."
PRINT "DESIRABLE: Math processor, speed of 10O MHz or greater, color VGA or EGA,"
PRINT "           graphics card with enough memory to hold two graphics screens,"
PRINT "           and a hard disk."
PRINT "NOTES: (1) You never need to wait for a drawing to be completed;"
PRINT "           enter your commands as soon as you can see enough to work."
PRINT "       (2) Faults and element sides are really arcs of great circles;"
PRINT "           however, they will be shown as straight lines for drawing speed.";
'------------------------------ graphics of globe and sample elements -----------------
WINDOW (-2, -1.5)-(2, 1.5)  'upper right quarter is now 0<x<2, 0<y<1.5
IF Colored% THEN COLOR Blue%
'longitude lines
FOR Phi! = 10! TO 90! STEP 10!
	aspect! = PixelRatio! / COS((Phi! / 180!) * 3.14159)
	IF aspect! >= 1! THEN
		CIRCLE (1.25, .75), .74 * PixelRatio!, , , , aspect!
	ELSE
		radius! = .74 * COS((Phi! / 180!) * 3.14159)
		CIRCLE (1.25, .75), radius!, , , , aspect!
	END IF
NEXT Phi!
'latitude lines
FOR Theta! = 10! TO 170! STEP 10!
	angle! = (Theta! / 180!) * 3.14159
	sa! = SIN(angle!)
	ca! = COS(angle!)
	X1! = 1.25 - .74 * sa!
	X2! = 1.25 + .74 * sa!
	Y1! = .75 + .74 * ca!
	Y2! = Y1!
	LINE (X1!, Y1!)-(X2!, Y2!)
NEXT Theta!
'generic grid of elements
X1! = 1.25
Y1! = .75
FOR Arg! = 0! TO 360! STEP 60!
	angle! = (Arg! / 180!) * 3.14159
	sa! = SIN(angle!)
	ca! = COS(angle!)
	X2! = 1.25 + .4 * ca!
	Y2! = .75 + .4 * sa!
	IF Colored% THEN COLOR Green%
	LINE (X1!, Y1!)-(X2!, Y2!)
	IF Colored% THEN
		IF (Arg! = 60!) OR (Arg! = 240!) THEN COLOR Red%
	END IF
	LINE (X2!, Y2!)-(1.25, .75)
	IF Colored% THEN COLOR Yellow%
	CIRCLE (X2!, Y2!), .03
	X1! = X2!
	Y1! = Y2!
NEXT Arg!
IF Colored% THEN COLOR Yellow%
CIRCLE (1.25, .75), .03
'outer circle
IF Colored% THEN COLOR White%
CIRCLE (1.25, .75), .74
WINDOW   'restore original coordinates
'--------------------------------------------------------------------------------------
IF Colored% THEN COLOR Foreground%
LOCATE Lines%, 60
PRINT "Press any key ...";
CALL WaitForKey(a$)

END SUB

SUB Flipped (Ins, i%, NODES%(), NodeABG!(), Outs, EMemo%())
'sets the 2^0 bit of EMemo% to 1 if element in inverted; otherwise to 0.
DIM V12!(3), V23!(3), V31!(3), VCross!(3)
n1% = NODES%(1, i%): N2% = NODES%(2, i%): N3% = NODES%(3, i%)
V12!(1) = NodeABG!(1, N2%) - NodeABG!(1, n1%)
V12!(2) = NodeABG!(2, N2%) - NodeABG!(2, n1%)
V12!(3) = NodeABG!(3, N2%) - NodeABG!(3, n1%)
V23!(1) = NodeABG!(1, N3%) - NodeABG!(1, N2%)
V23!(2) = NodeABG!(2, N3%) - NodeABG!(2, N2%)
V23!(3) = NodeABG!(3, N3%) - NodeABG!(3, N2%)
V31!(1) = NodeABG!(1, n1%) - NodeABG!(1, N3%)
V31!(2) = NodeABG!(2, n1%) - NodeABG!(2, N3%)
V31!(3) = NodeABG!(3, n1%) - NodeABG!(3, N3%)
CALL Cross(Ins, V12!(), V23!(), Outs, VCross!())
AreaSign! = VCross!(1) * NodeABG!(1, n1%) + VCross!(2) * NodeABG!(2, n1%) + VCross!(3) * NodeABG!(3, n1%)
IF AreaSign! > 0! THEN
	EMemo%(i%) = 0   'False that there is any problem with element.
ELSE
	EMemo%(i%) = 1   'True (an unconventional form using only one bit)
END IF
END SUB

FUNCTION GETBIT% (NUMBER%, PLACE%)
'RETURNS 0 (False) OR -1 (True) FROM THE 2**PLACE% POSITION, IN NUMBER%
IF (NUMBER% AND CINT(2 ^ PLACE%)) THEN
	GETBIT% = NOT 0   'True
ELSE
	GETBIT% = 0       'False
END IF
END FUNCTION

SUB GetFileName (Text$, LineN%, NewName$)
GetFileN:
LOCATE LineN%, 3: PRINT Text$;
LOCATE LineN%, LEN(Text$) + 3: PRINT "        ";
LOCATE LineN%, LEN(Text$) + 3: INPUT ; "", FileN$
IF LEN(FileN$) < 1 OR LEN(FileN$) > 8 THEN GOTO GetFileN
NewName$ = ""
FOR i% = 1 TO LEN(FileN$)
	T$ = MID$(FileN$, i%, 1)
	T% = ASC(T$)
	OK1% = (T% = 33) OR ((T% >= 35) AND (T% <= 41))  '!#$%&()
	OK2% = (T% = 44) OR (T% = 45)     ',-
	OK3% = (T% >= 48) AND (T% <= 57)  '0123456789
	OK4% = (T% >= 64) AND (T% <= 90)  '@ABC...XYZ
	OK5% = (T% >= 97) AND (T% <= 123)  'abc...xyz{
	OK6% = (T% = 125) OR (T% = 126)    '}~
	OK% = OK1% OR OK2% OR OK3% OR OK4% OR OK5% OR OK6%
	IF OK5% THEN T$ = UCASE$(T$)
	IF OK% THEN NewName$ = NewName$ + T$ ELSE NewName$ = NewName$ + "_"
NEXT i%
LOCATE LineN%, LEN(Text$) + 3: PRINT "        ";
LOCATE LineN%, LEN(Text$) + 3: PRINT NewName$;
WaitTime! = 1!: StartTime! = TIMER: DO: TimeNow! = TIMER: LOOP UNTIL (TimeNow! - StartTime!) > WaitTime!
END SUB

SUB GetNet (Ins, GridFileN$, MXNODE%, MXEL%, MXFEL%, Outs, TITLE$, NUMNOD%, NREALN%, NFAKEN%, N1000%, NodeABG!(), EQCM!(), NMemo%(), NUMEL%, NODES%(), Mus!(), EMemo%(), NFL%, NODEF%(), FMemo%(), FDIP!(), OFFSET!(), FAILED%, TempVec!())
'read in a finite element grid, without checking for topology
'fake nodes are renumbered as just above real ones to save storage
DIM DIPS!(2), IFN%(4), V12!(3), V23!(3), V31!(3), VCross!(3), tMus!(3)
FAILED% = NOT 0  'in case of abend within routine
ASum# = 0#: BSum# = 0#: GSum# = 0#
CLOSE #1
OPEN GridFileN$ FOR INPUT AS #1
ln% = 0: LOCATE 24, 25: PRINT " Reading Line #";
ln% = ln% + 1: LOCATE 24, 41: PRINT ln%;
LINE INPUT #1, TITLE$
ln% = ln% + 1: LOCATE 24, 41: PRINT ln%;
INPUT #1, NUMNOD%, NREALN%, NFAKEN%, N1000%
IF (CLNG(N1000%) + CLNG(NFAKEN%)) > 32767& THEN
	CALL Blanker(16, 24, 2, 79)
	LOCATE 17, 8: PRINT "HIGHEST NODE NUMBERS WILL EXCEED VALUES OF 32,767 SET BY"
	LOCATE 18, 8: PRINT "THE STORAGE OF INTEGERS IN BASIC.  REDUCE N1000% OR USE SMALLER GRID."
	FAILED% = NOT 0: EXIT SUB
END IF
CALL Finish(1)
IF (NUMNOD% <> (NREALN% + NFAKEN%)) THEN
	CALL Blanker(16, 24, 2, 79)
	LOCATE 17, 8: PRINT "INCONSISTENT DATA:";
	LOCATE 18, 8: PRINT "NUMBER OF NODES SHOULD EQUAL TOTAL OF REAL NODES AND FAKE NODES.";
	FAILED% = NOT 0: EXIT SUB
END IF
IF (NUMNOD% > MXNODE%) THEN
	CALL Blanker(16, 24, 2, 79)
	LOCATE 17, 8: PRINT "NUMBER OF NODES IN GRID ("; NUMNOD%; ") WOULD EXCEED"
	LOCATE 18, 8: PRINT "THE LIMIT ("; MXNODE%; ") SET BY MEMORY.";
	FAILED% = NOT 0: EXIT SUB
END IF
IF (NREALN% > N1000%) THEN
	CALL Blanker(16, 24, 2, 79)
	LOCATE 17, 8: PRINT "NUMBER OF REAL NODES IN GRID ("; NREALN%; ") EXCEEDS"
	LOCATE 18, 8: PRINT "PARAMETER N1000% ("; N1000%; ") IN SAME FILE; ERROR.";
	FAILED% = NOT 0: EXIT SUB
END IF
FOR k% = 1 TO NUMNOD%
	ln% = ln% + 1
	LOCATE 24, 41: PRINT ln%;
	INPUT #1, i%, Lon!, Lat!, ELEVI!, QI!, ZMI!, HLI!
	CALL Finish(1)
	IF i% > NREALN% THEN i% = i% - N1000% + NREALN%
	Lon! = Lon! * .017453293#
	Lat! = Lat! * .017453293#
	CosLat! = COS(Lat!)
	NodeABG!(1, i%) = CosLat! * COS(Lon!)
	NodeABG!(2, i%) = CosLat! * SIN(Lon!)
	NodeABG!(3, i%) = SIN(Lat!)
	ASum# = ASum# + NodeABG!(1, i%)
	BSum# = BSum# + NodeABG!(2, i%)
	GSum# = GSum# + NodeABG!(3, i%)
	EQCM!(1, i%) = ELEVI!
	EQCM!(2, i%) = QI!
	EQCM!(3, i%) = ZMI!
	EQCM!(4, i%) = HLI!
	NMemo%(i%) = 0
NEXT k%
ASum# = ASum# / CDBL(NUMNOD%)
BSum# = BSum# / CDBL(NUMNOD%)
GSum# = GSum# / CDBL(NUMNOD%)
CALL Unitize(Mods, ASum#, BSum#, GSum#)
TempVec!(1) = ASum#: TempVec!(2) = BSum#: TempVec!(3) = GSum#
ln% = ln% + 1: LOCATE 24, 41: PRINT ln%;
INPUT #1, NUMEL%
CALL Finish(1)
IF (NUMEL% > MXEL%) THEN
	CALL Blanker(16, 24, 2, 79)
	LOCATE 17, 8: PRINT "NUMBER OF ELEMENTS IN GRID ("; NUMEL%; ") WOULD EXCEED"
	LOCATE 18, 8: PRINT "THE LIMIT ("; MXEL%; ") SET BY MEMORY.";
	FAILED% = NOT 0: EXIT SUB
END IF
FOR k% = 1 TO NUMEL%
	ln% = ln% + 1
	LOCATE 24, 41: PRINT ln%;
	IF (FEP% = 3) THEN
		INPUT #1, i%, IFN%(1), IFN%(2), IFN%(3), tMus!(1), tMus!(2), tMus!(3)
	ELSE
		INPUT #1, i%, IFN%(1), IFN%(2), IFN%(3)
	END IF
	CALL Finish(1)
	IF NUMNOD% > NREALN% THEN
		FOR j% = 1 TO 3
			IF IFN%(j%) > NREALN% THEN IFN%(j%) = IFN%(j%) - N1000% + NREALN%
		NEXT j%
	END IF
	FOR j% = 1 TO 3
		NODES%(j%, i%) = IFN%(j%)
		IF (FEP% = 3) THEN Mus!(j%, i%) = tMus!(j%)
	NEXT j%
NEXT k%
ln% = ln% + 1: LOCATE 24, 41: PRINT ln%;
IF (FEP% > 1) THEN
	NFL% = 0
ELSE
	INPUT #1, NFL%
	CALL Finish(1)
	IF (NFL% > MXFEL%) THEN
		CALL Blanker(16, 24, 2, 79)
		LOCATE 17, 8: PRINT "NUMBER OF FAULTS IN GRID ("; NFL%; ") WOULD EXCEED"
		LOCATE 18, 8: PRINT "THE LIMIT ("; MXFEL%; ") SET BY MEMORY.";
		FAILED% = NOT 0: EXIT SUB
	END IF
	FOR k% = 1 TO NFL%
		ln% = ln% + 1
		LOCATE 24, 41: PRINT ln%;
		INPUT #1, i%, IFN%(1), IFN%(2), IFN%(3), IFN%(4), DIPS!(1), DIPS!(2), OFFST!
		CALL Finish(1)
		IF NUMNOD% > NREALN% THEN
			FOR j% = 1 TO 4
				IF IFN%(j%) > NREALN% THEN IFN%(j%) = IFN%(j%) - N1000% + NREALN%
			NEXT j%
		END IF
		FOR j% = 1 TO 4: NODEF%(j%, i%) = IFN%(j%): NEXT j%
		FOR L% = 1 TO 2
			IF DIPS!(L%) > 0! THEN
				FDIP!(L%, i%) = DIPS!(L%)
			ELSE
				FDIP!(L%, i%) = 180! + DIPS!(L%)
			END IF
		NEXT L%
		OFFSET!(i%) = OFFST!
	NEXT k%
END IF
CLOSE #1
CALL Blanker(24, 24, 2, 79)
LOCATE 24, 12: PRINT "Testing elements for counterclockwise perimeter...";
'and making the 2^0 bit of EMemo% = 0 if area is positive; = 1 if inverted.
FOR i% = 1 TO NUMEL%
	LOCATE 24, 65: PRINT i%;
	CALL Flipped(Ins, i%, NODES%(), NodeABG!(), Outs, EMemo%())
NEXT i%
FAILED% = 0 'False; normal completion
END SUB

FUNCTION Inside% (TempVec!(), NPlanes%, Normals!())
'Determine whether position TempVec!(3) = (Alpha,Beta,Gamma)
'is on the correct (outward normal) side of ALL the NPlanes%
'great circles whose normals are in Normals!(1-3,1-NPoly%).
'Returns value of True or False.
FOR i% = 1 TO NPoly%
	Test! = TempVec!(1) * Normals!(1, i%) + TempVec!(2) * Normals!(2, i%) + TempVec!(3) * Normals!(3, i%)
	IF Test! < 0! THEN
		Inside% = 0  'False
		EXIT FUNCTION
	END IF
NEXT i%
Inside% = NOT 0   'True
END FUNCTION

SUB NEXTto (Ins, i%, j%, NFL%, NODEF%(), NODES%(), NUMEL%, Outs, KFAULT%, jf%, KELE%, je%)
'
'  DETERMINE WHETHER THERE ARE MORE ELEMENTS ADJACENT TO SIDE J% OF
'  TRIANGULAR CONTINUUM ELEMENT I%.
'    J% = 1 MEANS THE SIDE OPPOSITE NODE # NODES%(1,I%).
'    J% = 2 MEANS THE SIDE OPPOSITE NODE # NODES%(2,I%).
'    J% = 3 MEANS THE SIDE OPPOSITE NODE # NODES%(3,I%).
'  IF A FAULT ELEMENT IS ADJACENT, ITS NUMBER IS KFAULT%; OTHERWISE,
'  KFAULT% IS SET TO ZERO.
'  IF ANOTHER TRIANGULAR CONTINUUM ELEMENT IS ADJACENT (EVEN ACROSS
'  FAULT ELEMENT KFAULT%) THEN ITS NUMBER IS KELE%; OTHERWISE, KELE% = 0.
'
'  TWO NODE NUMBERS ALONG THE SIDE OF INTEREST, COUNTERCLOCKWISE:
	  n1% = NODES%((j% MOD 3) + 1, i%)
	  N2% = NODES%(((j% + 1) MOD 3) + 1, i%)
'  CHECK FOR ADJACENT FAULT ELEMENT FIRST:
	  FOUNDF% = 0
	  FOR k% = 1 TO NFL%
		  M1% = NODEF%(1, k%)
		  M2% = NODEF%(2, k%)
		  M3% = NODEF%(3, k%)
		  M4% = NODEF%(4, k%)
		  IF (M1% = N2%) AND (M2% = n1%) THEN
			  FOUNDF% = NOT 0 'True
			  KFAULT% = k%
			  jf% = 1
			  GOTO L11
		  END IF
		  IF (M3% = N2%) AND (M4% = n1%) THEN
			  FOUNDF% = NOT 0 'True
			  KFAULT% = k%
			  jf% = 2
			  GOTO L11
		  END IF
	  NEXT k%
L11:   IF NOT FOUNDF% THEN KFAULT% = 0
'  IF THERE WAS A FAULT, REPLACE 2 NODE NUMBERS THAT WE SEARCH FOR:
	  IF FOUNDF% THEN
		  IF jf% = 1 THEN
			  n1% = M3%
			  N2% = M4%
		  ELSE 'jf% = 2
			  n1% = M1%
			  N2% = M2%
		  END IF
	  END IF
'  SEARCH FOR ADJACENT TRIANGULAR CONTINUUM ELEMENT:
	  KELE% = 0
	  KLow% = i%
	  KHigh% = i%
	 'begin irregular loop, checking nearest neighbors first:
E1:
		  KLow% = KLow% - 1
		  IF (KLow% > 0) THEN
			  FOR L% = 1 TO 3
				  M1% = NODES%((L% MOD 3) + 1, KLow%)
				  M2% = NODES%(((L% + 1) MOD 3) + 1, KLow%)
				  IF (M1% = N2%) AND (M2% = n1%) THEN
					  KELE% = KLow%
					  je% = L%
					  EXIT SUB
				  END IF
			  NEXT L%
		  END IF
		  KHigh% = KHigh% + 1
		  IF (KHigh% <= NUMEL%) THEN
			  FOR L% = 1 TO 3
				  M1% = NODES%((L% MOD 3) + 1, KHigh%)
				  M2% = NODES%(((L% + 1) MOD 3) + 1, KHigh%)
				  IF (M1% = N2%) AND (M2% = n1%) THEN
					  KELE% = KHigh%
					  je% = L%
					  EXIT SUB
				  END IF
			  NEXT L%
		  END IF
	  IF ((KLow% > 1) OR (KHigh% < NUMEL%)) GOTO E1
END SUB

SUB OnAnyE (Ins, n%, n1%, NLast%, NODES%(), Outs, ie%, je%)
'locate any element including the specified node n%
ie% = 0
je% = 0
FOR i% = n1% TO NLast%
	FOR j% = 1 TO 3
		IF NODES%(j%, i%) = n% THEN
			ie% = i%
			je% = j%
			EXIT SUB
		END IF
	NEXT j%
NEXT i%
END SUB

SUB OnAnyF (Ins, n%, n1%, NLast%, NODEF%(), Outs, ie%, je%)
'locate any fault including the specified node n%
ie% = 0
je% = 0
FOR i% = n1% TO NLast%
	FOR j% = 1 TO 4
		IF NODEF%(j%, i%) = n% THEN
			ie% = i%
			je% = j%
			EXIT SUB
		END IF
	NEXT j%
NEXT i%
END SUB

SUB Pixels (Ins, X!, Y!, Outs, col%, row%)
'convert (x,y) coordinates to pixels; prevent overflows.
col! = Scales!(1, 1) * X! + Scales!(1, 2)
row! = Scales!(2, 1) * Y! + Scales!(2, 2)
IF ABS(col!) < 32767 THEN
	col% = CINT(col!)
ELSEIF col! > 0! THEN
	col% = 32767
ELSE
	col% = -32767
END IF
IF ABS(row!) < 32767 THEN
	row% = CINT(row!)
ELSEIF row! > 0! THEN
	row% = 32767
ELSE
	row% = -32767
END IF
END SUB

FUNCTION Principal! (angle!)
'Returns the principal value (-Pi/2 < Prinicpal! <= +Pi/2) associated with any
'angle; vectors pointing left are reversed.
CONST Pi! = 3.141592
CONST HalfPi! = 1.570796
CONST Scant! = -1.570795
T! = angle!
DO WHILE T! > HalfPi!
	T! = T! - Pi!
LOOP
DO WHILE T! <= Scant!
	T! = T! + Pi!
LOOP
IF T! > HalfPi! THEN T! = HalfPi!
Principal! = T!
END FUNCTION

SUB PutNet (GridFileN$, TITLE$, NUMNOD%, NREALN%, NFAKEN%, N1000%, NodeABG!(), EQCM!(), NUMEL%, NODES%(), Mus!(), NFL%, NODEF%(), FDIP!(), OFFSET!())
'write out finite element grid to a file, in format used by
'the SHELLS finite element program; add Mus! after element definition for Restore3.
DIM IFN%(4), DIPS!(2), TempVec!(3)
CLOSE #1
OPEN GridFileN$ FOR OUTPUT AS #1
ln% = 0: LOCATE 24, 25: PRINT "Writing Line #";
ln% = ln% + 1: LOCATE 24, 41: PRINT ln%;
PRINT #1, TITLE$
ln% = ln% + 1: LOCATE 24, 41: PRINT ln%;
PRINT #1, USING "##### ##### ##### ##### F"; NUMNOD%; NREALN%; NFAKEN%; N1000%
FOR k% = 1 TO NREALN%
	ln% = ln% + 1: LOCATE 24, 41: PRINT ln%;
	TempVec!(1) = NodeABG!(1, k%): TempVec!(2) = NodeABG!(2, k%): TempVec!(3) = NodeABG!(3, k%)
	CALL ABG2LonLat(Ins, TempVec!(), Outs, Lon!, Lat!)
	Lon! = Lon! * 57.2957795#: Lat! = Lat! * 57.2957795#
	PRINT #1, USING "##### +###.### +##.### +#.##^^^^ +#.##^^^^ +#.##^^^^ +#.##^^^^"; k%; Lon!; Lat!; EQCM!(1, k%); EQCM!(2, k%); EQCM!(3, k%); EQCM!(4, k%)
NEXT k%
IF NUMNOD% > NREALN% THEN
	FOR k% = (NREALN% + 1) TO NUMNOD%
		i% = N1000% + k% - NREALN%
		ln% = ln% + 1: LOCATE 24, 41: PRINT ln%;
		TempVec!(1) = NodeABG!(1, k%): TempVec!(2) = NodeABG!(2, k%): TempVec!(3) = NodeABG!(3, k%)
		CALL ABG2LonLat(Ins, TempVec!(), Outs, Lon!, Lat!)
		Lon! = Lon! * 57.2957795#: Lat! = Lat! * 57.2957795#
		PRINT #1, USING "##### +###.### +##.### +#.##^^^^ +#.##^^^^ +#.##^^^^ +#.##^^^^"; i%; Lon!; Lat!; EQCM!(1, k%); EQCM!(2, k%); EQCM!(3, k%); EQCM!(4, k%)
	NEXT k%
END IF
ln% = ln% + 1: LOCATE 24, 41: PRINT ln%;
PRINT #1, USING "#####"; NUMEL%
FOR k% = 1 TO NUMEL%
	FOR j% = 1 TO 3: IFN%(j%) = NODES%(j%, k%): NEXT j%
	IF NUMNOD% > NREALN% THEN
		FOR j% = 1 TO 3
			IF IFN%(j%) > NREALN% THEN IFN%(j%) = N1000% + IFN%(j%) - NREALN%
		NEXT j%
	END IF
	ln% = ln% + 1: LOCATE 24, 41: PRINT ln%;
	IF (FEP% = 3) THEN
		PRINT #1, USING "##### ##### ##### ##### #.#^^^^ ####.# #.#^^^^"; k%; IFN%(1); IFN%(2); IFN%(3); Mus!(1, k%); Mus!(2, k%); Mus!(3, k%)
	ELSE
		PRINT #1, USING "##### ##### ##### #####"; k%; IFN%(1); IFN%(2); IFN%(3)
	END IF
NEXT k%
ln% = ln% + 1: LOCATE 24, 41: PRINT ln%;
PRINT #1, USING "#####"; NFL%
FOR i% = 1 TO NFL%
	FOR j% = 1 TO 4: IFN%(j%) = NODEF%(j%, i%): NEXT j%
	IF NUMNOD% > NREALN% THEN
		FOR j% = 1 TO 4
			IF IFN%(j%) > NREALN% THEN IFN%(j%) = N1000% + IFN%(j%) - NREALN%
		NEXT j%
	END IF
	FOR j% = 1 TO 2
		IF FDIP!(j%, i%) <= 90.01 THEN
			DIPS!(j%) = FDIP!(j%, i%)
		ELSE
			DIPS!(j%) = FDIP!(j%, i%) - 180!
		END IF
	NEXT j%
	ln% = ln% + 1: LOCATE 24, 41: PRINT ln%;
	PRINT #1, USING "##### ##### ##### ##### ##### +##.# +##.# +#.##^^^^"; i%; IFN%(1); IFN%(2); IFN%(3); IFN%(4); DIPS!(1); DIPS!(2); OFFSET!(i%)
NEXT i%
CLOSE #1
END SUB

SUB RaiseOne (Ins, VCol%(), VRow%())
'Colors in an element, by increasing the color-code of every pixel one unit.
'Note: For routine to succeed, vertices must be input in counterclockwise
'order.
DIM X!(4), Y!(4), Use%(3, 2)
FALSE% = 0: TRUE% = NOT 0
Row1% = 32000: Row2% = -32000
FOR j% = 1 TO 3
	IF VRow%(j%) > Row2% THEN Row2% = VRow%(j%)
	IF VRow%(j%) < Row1% THEN Row1% = VRow%(j%)
NEXT j%
IF Row1% < 0 THEN Row1% = 0
IF Row1% > HiRow% THEN Row1% = HiRow%
IF Row2% < 0 THEN Row2% = 0
IF Row2% > HiRow% THEN Row2% = HiRow%
IF Row2% <= Row1% THEN EXIT SUB
IF VRow%(1) <> VRow%(2) THEN
	IF VRow%(2) > VRow%(1) THEN
		Use%(1, 1) = TRUE%
		Use%(1, 2) = FALSE%
	ELSE
		Use%(1, 1) = FALSE%
		Use%(1, 2) = TRUE%
	END IF
ELSE
	Use%(1, 1) = FALSE%
	Use%(1, 2) = FALSE%
END IF
IF VRow%(2) <> VRow%(3) THEN
	IF VRow%(3) > VRow%(2) THEN
		Use%(2, 1) = TRUE%
		Use%(2, 2) = FALSE%
	ELSE
		Use%(2, 1) = FALSE%
		Use%(2, 2) = TRUE%
	END IF
ELSE
	Use%(2, 1) = FALSE%
	Use%(2, 2) = FALSE%
END IF
IF VRow%(3) <> VRow%(1) THEN
	IF VRow%(1) > VRow%(3) THEN
		Use%(3, 1) = TRUE%
		Use%(3, 2) = FALSE%
	ELSE
		Use%(3, 1) = FALSE%
		Use%(3, 2) = TRUE%
	END IF
ELSE
	Use%(3, 1) = FALSE%
	Use%(3, 2) = FALSE%
END IF
X!(1) = CSNG(VCol%(1)): Y!(1) = CSNG(VRow%(1))
X!(2) = CSNG(VCol%(2)): Y!(2) = CSNG(VRow%(2))
X!(3) = CSNG(VCol%(3)): Y!(3) = CSNG(VRow%(3))
X!(4) = X!(1): Y!(4) = Y!(1)
FOR row% = Row1% TO (Row2% - 1)
	r! = CSNG(row%)
	Col1% = 0: Col2% = HiCol%
	FOR i% = 1 TO 3
		j% = i% + 1
		IF Use%(i%, 1) THEN
			Edge% = 1 + CINT(X!(i%) + (X!(j%) - X!(i%)) * (r! - Y!(i%)) / (Y!(j%) - Y!(i%)))
			IF Edge% > Col1% THEN Col1% = Edge%
		END IF
		IF Use%(i%, 2) THEN
			Edge% = CINT(X!(i%) + (X!(j%) - X!(i%)) * (r! - Y!(i%)) / (Y!(j%) - Y!(i%)))
			IF Edge% < Col2% THEN Col2% = Edge%
		END IF
	NEXT i%
	FOR col% = Col1% TO Col2%
		Hue% = POINT(col%, row%)
		Hue% = Hue% + 1
		PSET (col%, row%), Hue%
	NEXT col%
NEXT row%
END SUB

SUB RunMenu (Ins, Background%, BestMode%, CD$, a$, CW1%, CW2%, Foreground%, HCol%(), HotCol%, HotColor%, HotLine%, HRow%, Lines%, Mask%(), MXEL%, MXFEL%, MXNODE%, NFL%, NUMEL%, NUMNOD%, Places%(), RAM&, RW1%, RW2%, Trial$, Outs, Current$)
'operate the menu in an "endless" loop (until Enter is pressed)
IF LEN(a$) > 0 THEN
	IF ASC(a$) = 13 THEN 'ignore ENTER or an infinite loop occurs
		a$ = ""
	ELSE
		GOTO GetMoving
	END IF
END IF
DO
	'blank out the dialog box (unless more commands are coming in)
	a$ = INKEY$: IF LEN(a$) > 0 THEN GOTO GetMoving
	'place details in dialog box (unless more commands are coming in)
	IF Colored% THEN
		IF BestMode% <= 10 THEN
			COLOR Foreground%, Background%
			LINE (CW1%, RW1%)-(CW2%, RW2%), Brown%, BF
		ELSE
			COLOR Foreground%
			CALL Blanker(16, 24, 2, 79)
		END IF
	ELSE
		CALL Blanker(16, 24, 2, 79)
	END IF
	a$ = INKEY$: IF LEN(a$) > 0 THEN GOTO GetMoving
	IF Trial$ = "D" THEN
		LOCATE 16, 10
		PRINT "Check or change the current directory (CD of DOS),"
		LOCATE 17, 10
		PRINT "which will be the default directory when a grid file"
		LOCATE 18, 10
		PRINT "is Loaded or Saved, or when a Basemap file is displayed."
		IF LEFT$(CD$, 1) = "A" OR LEFT$(CD$, 1) = "B" THEN
			LOCATE 20, 3
			PRINT "* Working from a floppy disk drive means that Loading, Saving, and drawing";
			LOCATE 21, 5
			PRINT "the Basemap will be VERY SLOW.  You may wish to eXit this program,";
			LOCATE 22, 5
			PRINT "move your .DIG and .FEG files to a hard disk, and restart.";
		END IF
	ELSEIF Trial$ = "B" THEN
		LOCATE 16, 10
		PRINT "Designate a basemap file to be displayed along with"
		LOCATE 17, 10
		PRINT "the finite element grid, for location reference."
		LOCATE 18, 10
		PRINT "Basemap files are created with program DIGITISE, have extension"
		LOCATE 19, 10
		PRINT ".DIG, and consist of line segments and polylines only"
		LOCATE 20, 10
		PRINT "(no text, color, or shading).  Basemap files could include"
		LOCATE 21, 10
		PRINT "faults, surface geology, state lines, coastlines, latitude-"
		LOCATE 22, 10
		PRINT "longitude grids, et cetera.";
	ELSEIF Trial$ = "C" THEN
		IF Colored% THEN COLOR HotColor%
		LOCATE 17, 35
		PRINT "CAUTION !!!"
		LOCATE 18, 8
		PRINT "This command erases any finite element grid currently in memory."
		LOCATE 19, 26
		PRINT "It does NOT save it first!"
		LOCATE 20, 12
		PRINT "Press Enter ONLY if you are certain you wish to clear."
		LOCATE 21, 17
		PRINT "(Otherwise, select Save, or another command.)"
	ELSEIF Trial$ = "L" THEN
		LOCATE 16, 10
		PRINT "Load a finite element grid file into memory for display"
		LOCATE 17, 10
		PRINT "and/or editing.  (Note: A Load command will not execute when"
		LOCATE 18, 10
		PRINT "another grid is already in memory; you must Clear it first.)"
	ELSEIF Trial$ = "S" THEN
		LOCATE 16, 10
		PRINT "Save the finite element grid currently in memory"
		LOCATE 17, 10
		PRINT "as a file on the hard disk."
	ELSEIF Trial$ = "X" THEN
		LOCATE 16, 10
		PRINT "Exit and return to DOS."
		LOCATE 18, 10
		PRINT "(Note: This command will not execute while there is a finite"
		LOCATE 19, 17
		PRINT " element grid in memory; you must Clear it first.)"
	ELSEIF Trial$ = "W" THEN
		LOCATE 16, 10
		PRINT "Move the viewing window further or more precisely than is"
		LOCATE 17, 10
		PRINT "possible with the Zoom command, by specifying coordinates."
	ELSEIF Trial$ = "O" THEN
		LOCATE 16, 10
		PRINT "To create a second polar coordinate system, for convenience:"
		LOCATE 17, 10
		PRINT "Click the left mouse button on the desired Northern pole."
		LOCATE 18, 10
		PRINT "Click the right mouse button on another point to define the"
		LOCATE 19, 10
		PRINT "prime meridian.  The new (Lat',Lon') coordinates will"
		LOCATE 20, 10
		PRINT "be reported at lower right, in addition to old coordinates"
		LOCATE 21, 10
		PRINT "at lower left.  Grid and basemap files are unchanged.  This"
		LOCATE 22, 10
		PRINT "command can be used to measure distance between two points"
		LOCATE 23, 10
		PRINT "which will be the complement of Lat' after the second click.";
	ELSEIF Trial$ = "Z" THEN
		LOCATE 16, 10
		PRINT "After you supply a new window width below, you then click"
		LOCATE 17, 10
		PRINT "either mouse button on the point which is to be the center"
		LOCATE 18, 10
		PRINT "of the new window.  (If you leave the window width unchanged,"
		LOCATE 19, 10
		PRINT "this command can also be used for small window movements.)"
	ELSEIF Trial$ = "K" THEN
		LOCATE 16, 10
		PRINT "Choose the colors to be used in displaying the nodes,";
		LOCATE 17, 10
		PRINT "   elements, faults, outer circle, and basemap."
		LOCATE 19, 10
		PRINT "One choice is black, which will suppress plotting if desired."
	ELSEIF Trial$ = "R" THEN
		LOCATE 16, 10
		PRINT "Redraw the contents of the graphics window,"
		LOCATE 17, 10
		PRINT "to eliminate unwanted messages, secondary coordinates,"
		LOCATE 18, 10
		PRINT "or debris resulting from editing."
		LOCATE 20, 10
		PRINT "Also redraws this menu (when Enter is pressed for the second time)."
	ELSEIF Trial$ = "G" THEN
		LOCATE 16, 10
		PRINT "Covers the whole sphere with a regular grid of elements"
		LOCATE 17, 10
		PRINT "obtained by iterative subdivision of the icosahedron."
		LOCATE 18, 10
		PRINT "You will be allowed to choose the level of refinement."
		LOCATE 19, 10
		PRINT "This is one way to begin building a global grid; local"
		LOCATE 20, 10
		PRINT "editing can be performed with commands Tile, Adjust, Nodes,"
		LOCATE 21, 10
		PRINT "Elements, Faults, and Inclination after Global grid."
		LOCATE 22, 9
		PRINT "(This command is not permitted with a finite element grid in "
		LOCATE 23, 10
		PRINT "memory; you must Clear it first.)"
	ELSEIF Trial$ = "T" THEN
		LOCATE 16, 10
		PRINT "Define a convex polygon by clicking the left mouse button on"
		LOCATE 17, 10
		PRINT "the vertices in the counterclockwise direction; end the polygon"
		LOCATE 18, 10
		PRINT "by clicking the right mouse button near the initial vertex."
		LOCATE 19, 10
		PRINT "A regular grid of elements will be created within the"
		LOCATE 20, 10
		PRINT "region specified, such that no element crosses the boundary."
		LOCATE 21, 10
		PRINT "You will be prompted for the element size."
	ELSEIF Trial$ = "A" THEN
		LOCATE 16, 10
		PRINT "Adjust node locations by placing cursor over a node, holding"
		LOCATE 17, 10
		PRINT "either button down, and dragging to the desired position."
		LOCATE 18, 10
		PRINT "If an X symbol appears in any of the neighboring elements,"
		LOCATE 19, 10
		PRINT "then the node displacement was too extreme; reduce it."
	ELSEIF Trial$ = "N" THEN
		LOCATE 16, 10
		PRINT "Add a node by clicking the left mouse button."
		LOCATE 18, 10
		PRINT "Remove a node by clicking the right mouse button."
		LOCATE 19, 10
		PRINT "(You will not be permitted to remove any node which is"
		LOCATE 20, 10
		PRINT "part of any element or fault; you must heal any Fault(s)"
		LOCATE 21, 10
		PRINT "and delete any Element(s) first.)"
	ELSEIF Trial$ = "E" THEN
		LOCATE 16, 10
		PRINT "Add a triangular element by clicking the left mouse button"
		LOCATE 17, 10
		PRINT "on any 3 nodes in counterclockwise order."
		LOCATE 19, 10
		PRINT "Remove an element by clicking the right mouse button"
		LOCATE 20, 10
		PRINT "on the little triangle symbol in its center.  (If none, use Kolors.)"
		LOCATE 22, 9
		PRINT "(Note: You will not be permitted to disturb any node which"
		LOCATE 23, 10
		PRINT "is already involved in a fault.  You must heal the fault"
		LOCATE 24, 10
		PRINT "using the F command, redefine all Elements, and recut the Fault.)";
	ELSEIF Trial$ = "F" THEN
		LOCATE 16, 6
		PRINT "Add a fault segment by clicking the left mouse button"
		LOCATE 17, 6
		PRINT "on any boundary between two elements, at its midpoint."
		LOCATE 19, 6
		PRINT "Heal (eliminate) any fault by clicking the right mouse button."
		LOCATE 21, 6
		PRINT "NOTE: An isolated fault in the middle of the grid has no meaning."
		LOCATE 22, 6
		PRINT "Connect faults to a grid boundary (if any) or to another fault."
		LOCATE 23, 6
		PRINT "Isolated faults will dissapear when you exit this command."
	ELSEIF Trial$ = "I" THEN
		LOCATE 16, 7
		PRINT "By default, faults are vertical (strike-slip) when created."
		LOCATE 17, 7
		PRINT "To change, specify dip at 2 points per fault (at each end)"
		LOCATE 18, 7
		PRINT "by clicking left mouse button for shallow dip and right mouse"
		LOCATE 19, 7
		PRINT "button for steeper dip, on the side to which the fault dips."
		LOCATE 20, 7
		PRINT "Clicking both buttons at once returns that end to vertical dip."
		LOCATE 22, 7
		PRINT "Note: Dips < 27.5 indicate subduction zones; in SHELLS these"
		LOCATE 23, 7
		PRINT "faults will have an upper limit on their traction; in global"
		LOCATE 24, 7
		PRINT "models they will require boundary conditions.";
	ELSEIF Trial$ = "H" THEN
		LOCATE 16, 7
		PRINT "Heading (azimuth) of a fault may be adjusted by clicking on";
		LOCATE 17, 7
		PRINT "its center, holding the mouse button down, and dragging the mouse";
		LOCATE 18, 7
		PRINT "away along the desired heading.  The heading will be displayed";
		LOCATE 19, 7
		PRINT "at the bottom of the screen, in degrees clockwise from North.";
		LOCATE 20, 6
		PRINT "(If you do not drag the mouse, the old heading will not be changed.)";
	ELSEIF Trial$ = "P" THEN
		LOCATE 16, 4
		PRINT "This command tests the topology of the grid, specifically:"
		LOCATE 17, 5
		PRINT "(a) whether its boundary (if any) is simply connected, and"
		LOCATE 18, 5
		PRINT "(b) whether there are internal holes or overlap areas."
		LOCATE 19, 4
		PRINT "The actions taken are:"
		LOCATE 20, 5
		PRINT "(1) Locate a free element side (or determine that sphere is full);"
		LOCATE 21, 5
		PRINT "(2) Trace around the boundary (if any) back to the initial point;"
		LOCATE 22, 5
		PRINT "(3) Compute the area of the sphere within this perimeter (if any);"
		LOCATE 23, 5
		PRINT "(4) Compute the sum of the areas of the elements;"
		LOCATE 24, 5
		PRINT "(5) Report results of (3) and (4) for comparison.";
	ELSEIF Trial$ = "V" THEN
		LOCATE 16, 5
		PRINT "This graphical test is designed to show any areas where elements"
		LOCATE 17, 5
		PRINT "may fail to meet (gaps) or may cover the same area twice (overlaps)."
		LOCATE 18, 5
		PRINT "After this command is selected, the insides of the elements will be"
		LOCATE 19, 5
		PRINT "colored in all succeeding views (until you Redraw or edit the grid)."
		LOCATE 20, 5
		PRINT "A gap or overlap should show up as a different color."
		LOCATE 21, 10
		PRINT "*** THIS COMMAND IS ONLY AVAILABLE ON COLOR MONITORS! ***";
		LOCATE 22, 5
		PRINT "In this view, you can click on any element, and the color will change";
		LOCATE 23, 5
		PRINT "in all elements of the same plate (but color change will not bleed across";
		LOCATE 24, 5
		PRINT "faults.)";
	ELSEIF Trial$ = "M" THEN
		KRAM% = RAM& / 1024&
		LOCATE 17, 16
		PRINT "Memory available for finite element grid: "; KRAM%; " K"
		LOCATE 19, 27
		PRINT "Number loaded   Limit     % of Limit"
		Npc% = CINT((100& * NUMNOD%) / MXNODE%)
		Epc% = CINT((100& * NUMEL%) / MXEL%)
		IF (FEP% > 1) THEN
			Fpc% = 100
		ELSE
			Fpc% = CINT((100& * NFL%) / MXFEL%)
		END IF
		LOCATE 20, 20: PRINT USING "NODES:    #####        #####        ###%"; NUMNOD%; MXNODE%; Npc%
		LOCATE 21, 20: PRINT USING "ELEMENTS: #####        #####        ###%"; NUMEL%; MXEL%; Epc%
		LOCATE 22, 20: PRINT USING "FAULTS:   #####        #####        ###%"; NFL%; MXFEL%; Fpc%
	ELSEIF Trial$ = "Q" THEN
		LOCATE 16, 7
		PRINT "This command causes one of the four nodal variables";
		LOCATE 17, 7
		PRINT "(elevation for Shells or mu_ for Restore v.2; heat-flow;";
		LOCATE 18, 7
		PRINT "crustal thickness; or mantle lithosphere thickness)";
		LOCATE 19, 7
		PRINT "to be contoured in color (on color monitors) in all views, until";
		LOCATE 20, 7
		PRINT "you select Redraw.  By clicking on any node, you may check its";
		LOCATE 21, 7
		PRINT "value and (optionally) enter a new value.  If you don't want to";
		LOCATE 22, 7
		PRINT "change the value, just press Enter.  To repeat the last value";
		LOCATE 23, 7
		PRINT "typed, press ' [single quote] and Enter.";
	ELSEIF Trial$ = "U" THEN
		LOCATE 16, 7
		PRINT "This command causes one of the three element variables";
		LOCATE 17, 7
		PRINT "       recent mu_,   transition age,   ancient mu_";
		LOCATE 18, 7
		PRINT "   (used ONLY by programs Restore v.3 and RetroMap v.3)";
		LOCATE 19, 7
		PRINT "to be displayed in color (on color monitors) in all views, until";
		LOCATE 20, 7
		PRINT "you select Redraw.  By clicking on any element, you may check its";
		LOCATE 21, 7
		PRINT "value and (optionally) enter a new value.  If you don't want to";
		LOCATE 22, 7
		PRINT "change the value, just press Enter.  To repeat the last value";
		LOCATE 23, 7
		PRINT "typed, press ' [single quote] and Enter.";
	END IF
	CALL WaitForKey(a$)
GetMoving:
	NewHotCol% = HotCol%
	NewHotLine% = HotLine%
	IF LEN(a$) = 1 THEN  'letter command, or Esc, or Return
		a$ = UCASE$(a$)
		m% = ASC(a$)
		n% = m% - 64
		IF (n% > 0) AND (n% < 27) THEN
			IF Places%(n%, 1) > 0 THEN
				Trial$ = a$
				NewHotCol% = Places%(n%, 1)
				NewHotLine% = Places%(n%, 2)
			END IF
		ELSEIF m% = 27 THEN 'Esc key
			Trial$ = Current$
			n% = ASC(Current$) - 64
			NewHotCol% = Places%(n%, 1)
			NewHotLine% = Places%(n%, 2)
		ELSEIF m% = 13 THEN 'Enter key; ONLY EXIT FROM LOOP!
			Current$ = Trial$
			a$ = "" 'BLANK OUT COMMAND KEY TO PREVENT INFINITE LOOP!
			EXIT SUB
		ELSEIF m% = 50 THEN 'digit 2 is pseudonym for letter O.
			Trial$ = "O"
			NewHotCol% = Places%(15, 1)
			NewHotLine% = Places%(15, 2)
		END IF
	ELSE 'extended character = possible cursor command
		m% = ASC(RIGHT$(a$, 1))
		IF m% = 72 THEN      'Up cursor
			DO: NewHotLine% = ((NewHotLine% + 12) MOD 14) + 1
			LOOP UNTIL Letters(NewHotCol%, NewHotLine%) <> " "
		ELSEIF m% = 80 THEN  'Down cursor
			DO: NewHotLine% = ((NewHotLine% + 14) MOD 14) + 1
			LOOP UNTIL Letters(NewHotCol%, NewHotLine%) <> " "
		ELSEIF m% = 75 OR m% = 77 THEN 'Left/Right cursors
			NewHotCol% = HotCol% MOD 2 + 1
			IF Letters(NewHotCol%, NewHotLine%) = " " THEN
				NewHotCol% = HotCol%
				BEEP
			END IF
		END IF
	END IF
	'place highlight on tentative command choice
	PUT (HCol%(HotCol%), HRow%), Mask%, XOR 'erase old highlight
	HotCol% = NewHotCol%
	HotLine% = NewHotLine%
	HRow% = (HotLine% - 1) * (HiRow% + 1) / Lines%
	PUT (HCol%(HotCol%), HRow%), Mask%, XOR 'place new highlight
	Trial$ = Letters(HotCol%, HotLine%)
LOOP
END SUB

SUB Scaler
'Uses WindowHeight!, HiCol%, and HiRow% (from COMMON SHARED) to compute
'R2C!, Scales!(), and UnScale!() (in COMMON SHARED) which
'convert back and forth from (x,y) to (col,row).
R2C! = (.5 * WindowHeight) ^ 2 + (.666667 * WindowHeight!) ^ 2
Scales!(1, 1) = .75 * CSNG(HiCol%) / WindowHeight!
Scales!(1, 2) = .5 * CSNG(HiCol%)
Scales!(2, 1) = -CSNG(HiRow%) / WindowHeight!
Scales!(2, 2) = .5 * CSNG(HiRow%)
UnScale!(1, 1) = 1.33333 * WindowHeight! / CSNG(HiCol%)
UnScale!(1, 2) = -.66667 * WindowHeight!
UnScale!(2, 1) = -WindowHeight! / CSNG(HiRow%)
UnScale!(2, 2) = .5 * WindowHeight!
END SUB

SUB SETBIT (NUMBER%, PLACE%, ONOFF%)
'IN NUMBER%, SETS THE 2**PLACE% BIT TO VALUE ONOFF% (O OR 1)
IF ONOFF% THEN
	NUMBER% = NUMBER% OR CINT(2 ^ PLACE%)
ELSE
	NUMBER% = NUMBER% AND (NOT CINT(2 ^ PLACE%))
END IF
END SUB

SUB SetMenu (Ins, Background%, BestMode%, Current$, Foreground%, HCol%(), HotCol%, HotLine%, HighPage%, Lines%, Mask%(), Mods, NeedMenu%, Outs, HRow%, PageNow%, Places%(), Trial$)
'draws the basic menu of commands
Trial$ = Current$
SCREEN BestMode%, , HighPage%, HighPage%
PageNow% = HighPage%
IF Colored% THEN
	IF BestMode% <= 10 THEN
		COLOR Foreground%, Background%
	ELSE
		COLOR Foreground%
	END IF
END IF
NeedMenu% = 0
CLS
LOCATE 1, 10
PRINT "CHOOSE COMMAND BY CURSOR OR CAPITAL LETTER, THEN PRESS Enter:"
LINE (0, ((HiRow% + 1) / Lines%) - 1)-(HiCol%, ((HiRow% + 1) / Lines%) - 1), Foreground%
LOCATE 2, 9: PRINT "FILE COMMANDS:"
'LOCATE 2, 45: PRINT "GRID COMMANDS:" (now converted to a vertical label)
LOCATE 9, 9: PRINT "WINDOW COMMANDS:"
'initialize highlighting, and color the group labels
HCol%(1) = HiCol% * .1
HCol%(2) = HiCol% * .55
BoxLine% = 2
HRow% = (BoxLine% - 1) * (HiRow% + 1) / Lines%
PUT (HCol%(1), HRow%), Mask%, AND
PUT (HCol%(2), HRow%), Mask%, AND
BoxLine% = 9
HRow% = (BoxLine% - 1) * (HiRow% + 1) / Lines%
PUT (HCol%(1), HRow%), Mask%, AND
'list all commands, and create cross-reference tables
FOR i% = 1 TO 26: Places%(i%, 1) = 0: Places%(i%, 2) = 0: NEXT i%
FOR i% = 1 TO 14: Letters(1, i%) = " ": Letters(2, i%) = " ": NEXT i%
LOCATE 3, 10: PRINT "Directory (path) for files"
Letters(1, 3) = "D"
Places%(4, 1) = 1
Places%(4, 2) = 3
LOCATE 4, 10: PRINT "Basemap file"
Letters(1, 4) = "B"
Places%(2, 1) = 1
Places%(2, 2) = 4
LOCATE 5, 10: PRINT "Clear grid file"
Letters(1, 5) = "C"
Places%(3, 1) = 1
Places%(3, 2) = 5
LOCATE 6, 10: PRINT "Load  grid file"
Letters(1, 6) = "L"
Places%(12, 1) = 1
Places%(12, 2) = 6
LOCATE 7, 10: PRINT "Save  grid file"
Letters(1, 7) = "S"
Places%(19, 1) = 1
Places%(19, 2) = 7
LOCATE 8, 10: PRINT "eXit to DOS"
Letters(1, 8) = "X"
Places%(24, 1) = 1
Places%(24, 2) = 8
LOCATE 10, 10: PRINT "Window position"
Letters(1, 10) = "W"
Places%(23, 1) = 1
Places%(23, 2) = 10
LOCATE 11, 10: PRINT "2nd Origin (re)set"
Letters(1, 11) = "O"
Places%(15, 1) = 1
Places%(15, 2) = 11
LOCATE 12, 10: PRINT "Zoom in or out"
Letters(1, 12) = "Z"
Places%(26, 1) = 1
Places%(26, 2) = 12
LOCATE 13, 10: PRINT "Kolors"
Letters(1, 13) = "K"
Places%(11, 1) = 1
Places%(11, 2) = 13
LOCATE 14, 10: PRINT "Redraw window & menu"
Letters(1, 14) = "R"
Places%(18, 1) = 1
Places%(18, 2) = 14
LOCATE 2, 43: PRINT "G  Global grid of elements"
Letters(2, 2) = "G"
Places%(7, 1) = 2
Places%(7, 2) = 2
LOCATE 3, 43: PRINT "r  Tile a region with elements"
Letters(2, 3) = "T"
Places%(20, 1) = 2
Places%(20, 2) = 3
LOCATE 4, 43: PRINT "i  Adjust node position(s)"
Letters(2, 4) = "A"
Places%(1, 1) = 2
Places%(1, 2) = 4
LOCATE 5, 43: PRINT "d  Node(s): add or delete"
Letters(2, 5) = "N"
Places%(14, 1) = 2
Places%(14, 2) = 5
LOCATE 6, 43: PRINT "   Element(s): add or delete"
Letters(2, 6) = "E"
Places%(5, 1) = 2
Places%(5, 2) = 6
IF (FEP% > 1) THEN
	LOCATE 7, 43: PRINT "c"
	LOCATE 8, 43: PRINT "o"
	LOCATE 9, 43: PRINT "m"
ELSE
	LOCATE 7, 43: PRINT "c  Fault(s): cut or heal"
	Letters(2, 7) = "F"
	Places%(6, 1) = 2
	Places%(6, 2) = 7
	LOCATE 8, 43: PRINT "o  Inclination (dip) of faults"
	Letters(2, 8) = "I"
	Places%(9, 1) = 2
	Places%(9, 2) = 8
	LOCATE 9, 43: PRINT "m  Heading (azimuth) of faults"
	Letters(2, 9) = "H"
	Places%(8, 1) = 2
	Places%(8, 2) = 9
END IF
LOCATE 10, 43: PRINT "m  Perimeter and area tests"
Letters(2, 10) = "P"
Places%(16, 1) = 2
Places%(16, 2) = 10
LOCATE 11, 43: PRINT "a  View gaps/overlaps"
Letters(2, 11) = "V"
Places%(22, 1) = 2
Places%(22, 2) = 11
LOCATE 12, 43: PRINT "n  Memory check"
Letters(2, 12) = "M"
Places%(13, 1) = 2
Places%(13, 2) = 12
LOCATE 13, 43: PRINT "d  elevation/Q/crust/mantle"
Letters(2, 13) = "Q"
Places%(17, 1) = 2
Places%(17, 2) = 13
IF (FEP% = 3) THEN
	LOCATE 14, 43: PRINT "s  mU_ (view/edit element data)"
	Letters(2, 14) = "U"
	Places%(21, 1) = 2
	Places%(21, 2) = 14
ELSE
	LOCATE 14, 43: PRINT "s"
END IF
'
'frame the dialog box
LOCATE 15, 1: PRINT CHR$(201);
FOR i% = 2 TO 79: PRINT CHR$(205); : NEXT i%
PRINT CHR$(187);
FOR i% = 16 TO 24
LOCATE i%, 1: PRINT CHR$(186);
LOCATE i%, 80: PRINT CHR$(186);
NEXT i%
LOCATE 25, 1: PRINT CHR$(200);
FOR i% = 2 TO 79: PRINT CHR$(205); : NEXT i%
PRINT CHR$(188);
'initialize highlight
HRow% = (HotLine% - 1) * (HiRow% + 1) / Lines%
PUT (HCol%(HotCol%), HRow%), Mask%, XOR
END SUB

SUB Sterradians (Ins, V1!(), V2!(), V3!(), Outs, SolidAngle!)
'If V1,V2,V3 are unit vectors (in the alpha/beta/gamma Cartesian system)
'from the center of the sphere to the three vertices of a spherical
'triangle, then this sub returns the solid angle of the triangle in
'sterradians.  Note that the calculations are all double-precision, but
'the result is only single-precision.  This is due to the inherent ill-
'conditioning of the formula in the last line, where very similar
'numbers are differenced.
'If the 3 corners are given in counterclockwise order
'(as seen from outside the planet),
'the area will be positive; else, negative.
'This sign feature is needed for computation of model area
'by the perimeter method.
DIM VCross#(3), RA#(3), RB#(3), RC#(3), CmB#(3), BmA#(3), center#(3)
'--------------------------------
'Convert to double precision:
RA#(1) = CDBL(V1!(1)): RA#(2) = CDBL(V1!(2)): RA#(3) = CDBL(V1!(3))
RB#(1) = CDBL(V2!(1)): RB#(2) = CDBL(V2!(2)): RB#(3) = CDBL(V2!(3))
RC#(1) = CDBL(V3!(1)): RC#(2) = CDBL(V3!(2)): RC#(3) = CDBL(V3!(3))
'---------------------------------
'Make sure that these are truly unit vectors, after conversion:
lA# = SQR(RA#(1) * RA#(1) + RA#(2) * RA#(2) + RA#(3) * RA#(3))
RA#(1) = RA#(1) / lA#: RA#(2) = RA#(2) / lA#: RA#(3) = RA#(3) / lA#
lB# = SQR(RB#(1) * RB#(1) + RB#(2) * RB#(2) + RB#(3) * RB#(3))
RB#(1) = RB#(1) / lB#: RB#(2) = RB#(2) / lB#: RB#(3) = RB#(3) / lB#
lC# = SQR(RC#(1) * RC#(1) + RC#(2) * RC#(2) + RC#(3) * RC#(3))
RC#(1) = RC#(1) / lC#: RC#(2) = RC#(2) / lC#: RC#(3) = RC#(3) / lC#
'--------------------------------
'Compute lengths of triangle sides:
aDot# = RB#(1) * RC#(1) + RB#(2) * RC#(2) + RB#(3) * RC#(3)
VCross#(1) = RB#(2) * RC#(3) - RB#(3) * RC#(2)
VCross#(2) = RB#(3) * RC#(1) - RB#(1) * RC#(3)
VCross#(3) = RB#(1) * RC#(2) - RB#(2) * RC#(1)
aCross# = SQR(VCross#(1) * VCross#(1) + VCross#(2) * VCross#(2) + VCross#(3) * VCross#(3))
IF aCross# = 0# THEN
	SolidAngle! = 0!
	EXIT SUB
END IF
aSide# = ATAN2D#(aCross#, aDot#)
'
bDot# = RC#(1) * RA#(1) + RC#(2) * RA#(2) + RC#(3) * RA#(3)
VCross#(1) = RC#(2) * RA#(3) - RC#(3) * RA#(2)
VCross#(2) = RC#(3) * RA#(1) - RC#(1) * RA#(3)
VCross#(3) = RC#(1) * RA#(2) - RC#(2) * RA#(1)
bCross# = SQR(VCross#(1) * VCross#(1) + VCross#(2) * VCross#(2) + VCross#(3) * VCross#(3))
IF bCross# = 0# THEN
	SolidAngle! = 0!
	EXIT SUB
END IF
bSide# = ATAN2D#(bCross#, bDot#)
'
cDot# = RA#(1) * RB#(1) + RA#(2) * RB#(2) + RA#(3) * RB#(3)
VCross#(1) = RA#(2) * RB#(3) - RA#(3) * RB#(2)
VCross#(2) = RA#(3) * RB#(1) - RA#(1) * RB#(3)
VCross#(3) = RA#(1) * RB#(2) - RA#(2) * RB#(1)
cCross# = SQR(VCross#(1) * VCross#(1) + VCross#(2) * VCross#(2) + VCross#(3) * VCross#(3))
IF cCross# = 0# THEN
	SolidAngle! = 0!
	EXIT SUB
END IF
cSide# = ATAN2D#(cCross#, cDot#)
'--------------------------------
'determine half-angles, and then double them:
s# = (aSide# + bSide# + cSide#) / 2#
'
Denominator# = SIN(bSide#) * SIN(cSide#)
SinHalf# = SQR((SIN(s# - bSide#) * SIN(s# - cSide#)) / Denominator#)
CosHalf# = SQR((SIN(s#) * SIN(s# - aSide#)) / Denominator#)
Half# = ATAN2D#(SinHalf#, CosHalf#)
a# = 2# * Half#
'
Denominator# = SIN(cSide#) * SIN(aSide#)
SinHalf# = SQR((SIN(s# - cSide#) * SIN(s# - aSide#)) / Denominator#)
CosHalf# = SQR((SIN(s#) * SIN(s# - bSide#)) / Denominator#)
Half# = ATAN2D#(SinHalf#, CosHalf#)
b# = 2# * Half#
'
Denominator# = SIN(aSide#) * SIN(bSide#)
SinHalf# = SQR((SIN(s# - aSide#) * SIN(s# - bSide#)) / Denominator#)
CosHalf# = SQR((SIN(s#) * SIN(s# - cSide#)) / Denominator#)
Half# = ATAN2D#(SinHalf#, CosHalf#)
c# = 2# * Half#
'--------------------------------
'magnitude of solid angle, in sterradians:
SolidAngle! = CSNG((a# + b# + c#) - 3.141592654#)
'
'Consider possible change of sign:
'(a) first, get outward normal vector by cross-product of two sides:
BmA#(1) = RB#(1) - RA#(1): BmA#(2) = RB#(2) - RA#(2): BmA#(3) = RB#(3) - RA#(3)
CmB#(1) = RC#(1) - RB#(1): CmB#(2) = RC#(2) - RB#(2): CmB#(3) = RC#(3) - RB#(3)
VCross#(1) = BmA#(2) * CmB#(3) - BmA#(3) * CmB#(2)
VCross#(2) = BmA#(3) * CmB#(1) - BmA#(1) * CmB#(3)
VCross#(3) = BmA#(1) * CmB#(2) - BmA#(2) * CmB#(1)
'(b) second, get rough center position by averaging nodes:
center#(1) = RA#(1) + RB#(1) + RC#(1)
center#(2) = RA#(2) + RB#(2) + RC#(2)
center#(3) = RA#(3) + RB#(3) + RC#(3)
'(c) compare their general directions:
dot# = center#(1) * VCross#(1) + center#(2) * VCross#(2) + center#(3) * VCross#(3)
IF (dot# < 0#) THEN SolidAngle! = -SolidAngle!
END SUB

SUB UnitVec (Mods, TempVec!())
'Converts any vector to a unit vector.
R2! = TempVec!(1) * TempVec!(1) + TempVec!(2) * TempVec!(2) + TempVec!(3) * TempVec!(3)
IF R2! > 0! THEN
	r! = SQR(R2!)
	TempVec!(1) = TempVec!(1) / r!
	TempVec!(2) = TempVec!(2) / r!
	TempVec!(3) = TempVec!(3) / r!
ELSE
	TempVec!(1) = 1!
	TempVec!(2) = 0!
	TempVec!(3) = 0!
END IF
END SUB

SUB WaitForKey (a$)
'does nothing until a key is pressed; then returns the key as a$
a$ = ""
DO WHILE LEN(a$) = 0
	a$ = INKEY$
LOOP
END SUB

SUB XandY (Ins, col%, row%, Outs, X!, Y!)
'convert pixel coordinates to (x,y) using predefined transformation
X! = UnScale!(1, 1) * col% + UnScale!(1, 2)
Y! = UnScale!(2, 1) * row% + UnScale!(2, 2)
END SUB

SUB XORLine (X1%, Y1%, X2%, Y2%)
'Draws a line, like LINE, but XOR's each pixel;
'thus, two successive calls erase the line,
'leaving the drawing completely unchanged.
'Note: HiColor% is in COMMON SHARED.
DX% = X2% - X1%
DY% = Y2% - Y1%
IF DX% = 0 AND DY% = 0 THEN
	EXIT SUB
ELSEIF ABS(DX%) > ABS(DY%) THEN 'quasihorizontal line
	Slope! = CSNG(DY%) / CSNG(DX%)
	IF DX% > 0 THEN
		TheStep% = 1
	ELSE
		TheStep% = -1
	END IF
	FOR X% = X1% TO X2% STEP TheStep%
		Y% = CINT(CSNG(X% - X1%) * Slope! + Y1%)
		Hue% = POINT(X%, Y%)
		Hue% = HiColor% - Hue%
		PSET (X%, Y%), Hue%
	NEXT X%
ELSE 'quasivertical line
	Slope! = CSNG(DX%) / CSNG(DY%)
	IF DY% > 0 THEN
		TheStep% = 1
	ELSE
		TheStep% = -1
	END IF
	FOR Y% = Y1% TO Y2% STEP TheStep%
		X% = CINT(CSNG(Y% - Y1%) * Slope! + X1%)
		Hue% = POINT(X%, Y%)
		Hue% = HiColor% - Hue%
		PSET (X%, Y%), Hue%
	NEXT Y%
END IF
END SUB

SUB XY2ABG (Ins, X!, Y!, Outs, Outside%, TempVec!())
'Converts (X,Y) of orthographic projection plane to (alpha,beta,gamma)
'orthogonal Cartesian coordinates.
R2! = X! * X! + Y! * Y!
IF R2! >= 1! THEN
	Outside% = NOT 0 'true
	TempVec!(1) = 0!
	TempVec!(2) = 0!
	TempVec!(3) = 0!
ELSE
	Outside% = 0 'false
	OutOfPlane! = SQR(1! - R2!)
	TempVec!(1) = X! * WinRight!(1) + Y! * WinUp!(1) + OutOfPlane! * WinOut!(1)
	TempVec!(2) = X! * WinRight!(2) + Y! * WinUp!(2) + OutOfPlane! * WinOut!(2)
	TempVec!(3) = X! * WinRight!(3) + Y! * WinUp!(3) + OutOfPlane! * WinOut!(3)
END IF
END SUB

SUB XY2LonLat (Ins, Which%, X!, Y!, Outs, Outside%, Lon!, Lat!)
'Converts (x,y) in projection plane to (Lat,Lon) if x^2+y^2 < 1;
'   else declares Outside%.
'The switch Which%=1,2 allows use of either the primary or secondary
'   coordinate systems.
DIM TempVec!(3), TempV2!(3)
R2! = X! * X! + Y! * Y!
IF R2! >= 1! THEN
	Outside% = NOT 0 'true
	Lon! = 0!
	Lat! = 0!
ELSE
	Outside% = 0 'false
	OutOfPlane! = SQR(1! - R2!)
	TempVec!(1) = X! * WinRight!(1) + Y! * WinUp!(1) + OutOfPlane! * WinOut!(1)
	TempVec!(2) = X! * WinRight!(2) + Y! * WinUp!(2) + OutOfPlane! * WinOut!(2)
	TempVec!(3) = X! * WinRight!(3) + Y! * WinUp!(3) + OutOfPlane! * WinOut!(3)
	IF Which% = 1 THEN
		CALL ABG2LonLat(Ins, TempVec!(), Outs, Lon!, Lat!)
	ELSE
		FOR i% = 1 TO 3
			TempV2!(i%) = 0!
			FOR j% = 1 TO 3
				TempV2!(i%) = TempV2!(i%) + TempVec!(j%) * Cart2!(j%, i%)
			NEXT j%
		NEXT i%
		CALL ABG2LonLat(Ins, TempV2!(), Outs, Lon!, Lat!)
	END IF
END IF
END SUB