C* */ /* THREEDS.F77   15-Feb-86   G. Doughty				     */
C*****************************************************************************/
C									     */
C Subroutines to plot function of two variables given as a 100 by 100	     */
C array of real values.  Current version drives a Televideo 921 terminal     */
C with the graphics option.						     */
C									     */
C This code was derived from a program I obtained from a graduate chemistry  */
C student at the University of North Carolina at Chapel Hill in 1973.	     */
C (I wish I could give his full name but I only remember that his last name  */
C was Hoskins and that he was working for Dr. Peterson at the time.)  I've   */
C only modified the original code as necessary for compilation with	     */
C Digital Research's FORTRAN compiler (as well as indenting it to make it    */
C easier to read).  I've also added routines to replace the original plot    */
C functions (PLOT, PICSIZ, ORIGIN) available on the university's computation */
C center computer.  Finally, I've added comments where I've been able to     */
C fathom the operation of the routines (The original program is largely      */
C devoid of program... but then many of the programs I wrote at UNC are      */
C similarly unadorned).							     */
C									     */
C*****************************************************************************/

	SUBROUTINE THREED(A, N, M, K)
C*****************************************************************************/
C									     */
C A is a 100 by 100 real array containing the function values at uniform     */
C grid points.  This array need not be completely filled, however.  N and    */
C M define the limits of the first and second (respectively) array index.    */
C Thus the function occupies A(1 - N, 1 - M).  K is a mystery as the momemt. */
C The example main-line that drives this always passes the value of 3 for K. */
C Unfortunately the other arguments are passed in COMMON.  They are ANGA,    */
C ANGB, and HV.  ANGB appears to tilt the 3-D plot and ANGA appears to ro-   */
C tate it (I'm not sure about the order of these operations on the image).   */
C HV seems to be a height parameter which I have been setting to 5.   I      */
C thinks it controls how much perspective affect one gets.  Beware that      */
C greater or lesser values may cause the image not to fit within 640 x 240.  */
C The program works best if the real values are scalled to the range from    */
C -4.0 to +4.0.								     */
C									     */
C*****************************************************************************/

	COMMON ANGA, ANGB, HV, D
	COMMON DL, DM, DN, CX, CY, CZ, QX, QY, QZ, DD
	DIMENSION H(10), V(10), X(2), Y(2), Z(2), XP(8), A(100,100)
	CALL	PICSIZ(15.,10.)
	DL = -COS(ANGA) * COS(ANGB)
	DM = -SIN(ANGA) * COS(ANGB)
	DN = -SIN(ANGB)
	IF(ABS(DN).NE.1.0) GO TO 10
	    WRITE(6,20)
20	    FORMAT(1X,'******* You are attempting to look straight ',
+	    'down (or up) at the surface')
	    GO TO 2150
C
10	CONTINUE
	DD = 1.0/SQRT(1.0 - DN * DN)
	X(1) = 0
	X(2) = N
	Y(1) = 0
	Y(2) = M
	T = N
	IF(N .LT. M) T = M
	D = M * M + N * N + T * T
	D = SQRT(D)
	SCALE = 10.0 * D
	CX = -DL * SCALE
	CY = -DM * SCALE
	CZ = -DN * SCALE
	QX = CX + D * DL
	QY = CY + D * DM
	QZ = CZ + D * DN
2060	CALL MINMAX(A, N, M, Z)
	DOLR = T / (Z(2) - Z(1))
	DO 30 I = 1, N
	    DO 30 J = 1, M
		A(I,J) = (A(I,J) - Z(1)) * DOLR
30	CONTINUE
	Z(1) = 0.0
	Z(2) = T
2080	CALL CUBE(X, Y, Z, XP, H, V)
	DO 2130 I = 1, 8
	    H(I) = ((XP(I) - QX) * DM - (H(I) - QY)*DL) * DD
	    V(I) = (V(I) - QZ) * DD
2100	CALL MINMAX(H, 8, 1, H(9))
2120	CALL MINMAX(V, 8, 1, V(9))
2130	CONTINUE
	IF(ANGB .GT. 0.0) GOTO 2140
	    T = V(9)
	    V(9) = V(10)
	    V(10) = T
2140	CALL ORTHO (X, Y, A, N, M, H, V, K)
C	CALL PLOT(HV + 2.0, 0.0, -3)
	XXXX = HV + 2.
	CALL ORIGIN(XXXX, 0., 1)
2150	CONTINUE
	END
C
	SUBROUTINE MINMAX(A, N, M, Z)
	DIMENSION Z(1), A(100, 1)
C
1050	U = A(1, 1)
1060	V = U
C
1080	DO 1190 J = 1, M
1100	    DO 1180 I = 1, N
1120		IF(U - A(I,J)) 1150, 1150, 1130
1130		    U = A(I, J)
1150		IF(A(I, J) - V) 1180, 1180, 1160
1160		    V = A(I, J)
1180	    CONTINUE
1190	CONTINUE
1210	Z(1) = U
1220	Z(2) = V
1230	END
C
	SUBROUTINE CUBE(X, Y, Z, XP, H, V)
	DIMENSION X(1), Y(1), Z(1), H(1), V(1), XP(1)
C
50	L = 0
70	DO 180 I = 1, 2
90	    DO 170 J = 1, 2
110		DO 160 K = 1, 2
130		    L = L + 1
140		    CALL ROTATE(X(I), Y(J), Z(K), XP(L), H(L), V(L))
160		CONTINUE
170	    CONTINUE
180	CONTINUE
	END
C
	SUBROUTINE ROTATE(X, Y, Z, XP, YP, ZP)
	COMMON ANGA, ANGB, HV, D
	COMMON DL, DM, DN, CX, CY, CZ, QX, QY, QZ, DD
	DK = D/((X - CX) * DL + (Y - CY) * DM + (Z - CZ) * DN)
	XP = CX + DK * (X - CX)
	YP = CY + DK * (Y - CY)
	ZP = CZ + DK * (Z - CZ)
	END
C
	SUBROUTINE ORTHO(X, Y, A, N, M, H, V, K)
	COMMON ANGA, ANGB, HV, D
	COMMON DL, DM, DN, CX, CY, CZ, QX, QY, QZ, DD
C
	DIMENSION X(1), Y(1), H(1), V(1), A(100,100)
	INTEGER UP, DOWN, PEN, P, Q, P1, PO
C
C
C
	END = 1.0 / 16.0
C
C Can use 1/32 or 1/64 for finer interpolation
C
	UP = 1
	DOWN = 0
	SH = HV / (H(10) - H(9))
	SV = HV / (V(10) - V(9))
	MM = M
	NN = N
80	IF(K - 1) 100, 120, 100
100	    IF(K - 3) 1110, 120, 1110
120		L = 0
		LD = 1
		DDD = 0.5 * LD
140		DO 1060 J = 1, M
		    Q = 0
		    YJ = J
160		    DO 1030 I = 1, NN
			L = L + LD
			XI = L
			CALL CANSEE(A, XI, YJ, N, M, P)
			PEN = UP
			IF(P) 510, 520, 530
510			CONTINUE
C			IF(Q) 540, 550, 170
			IF(Q) 540, 550, 540
520			CONTINUE
			IF(Q) 610, 1020, 610
530			CONTINUE
C			IF(Q) 170, 550, 540
			IF(Q) 540, 550, 540
540			CONTINUE
			PEN = DOWN
			GO TO 170
550			CONTINUE
			IF(I .EQ. 1) GO TO 170
			DI = DDD
			TO = L - LD
			T = TO + DI
			P1 = Q
560			IF(ABS(DI) .LT. END) GO TO 570
			    CALL CANSEE(A, T, YJ, N, M, PO)
			    DI = DI * 0.5
			    IF(PO .EQ. 0) GO TO 565
				TO = T
				P1 = PO
				T = T - DI
				GO TO 560
565			    T = T + DI
			    GOTO 560
570			CONTINUE
			T = TO
			IF(P1 * P) 170, 170, 580
580			CONTINUE
590			CONTINUE
			ZP = A(L-LD,J) + (T-L+LD)*(A(L,J)-A(L-LD,J))/LD
			CALL ROTATE(T,YJ,ZP,XP,HH,VV)
			HH = ((XP - QX)*DM - (HH - QY) * DL) * DD
			VV = (VV - QZ) * DD
			HH = (HH - H(9)) * SH
			VV = (VV - V(9)) * SV
			CALL PLOT(HH, VV, PEN)
C600			PEN = 5 - PEN
600			PEN = 1 - PEN
			GOTO 170
C
610			CONTINUE
			PEN = DOWN
			DI = DDD
			TO = L - LD
			T = TO + DI
			P1 = Q
620			IF(ABS(DI) .LT. END) GO TO 630
			    CALL CANSEE(A, T, YJ, N, M, PO)
			    DI = DI * 0.5
			    IF(PO .EQ. 0) GO TO 625
				TO = T
				P1 = PO
				T = T + DI
				GO TO 620
625			    T = T - DI
			    GOTO 620
630			CONTINUE
			T = TO
			IF(P1 * Q) 600, 600, 590
170			CALL ROTATE(XI, YJ, A(L, J), XP, HH, VV)
			VV = (VV - QZ) * DD
			HH = ((XP-QX)*DM - (HH - QY)*DL) * DD
190			HH = (HH - H(9)) * SH
			VV = (VV - V(9)) * SV
			CALL PLOT(HH, VV, PEN)
1020			Q = P
1030		    CONTINUE
C
C
		    L = L + LD
		    LD = -LD
		    DDD = -DDD
C
1060		CONTINUE
C
C
1090		IF(K - 3) 2060, 1110, 2060
C
1110		CONTINUE
C
		L = 0
		LD = 1
		DDD = 0.5 * LD
1140		DO 2040 I = 1, N
		    XI = I
		    Q = 0
1160		    DO 2020 J = 1, MM
			L = L + LD
			YJ = L
			CALL CANSEE(A, XI, YJ, N, M, P)
			PEN = UP
			IF(P) 1510, 1520, 1530
1510			CONTINUE
C			IF(Q) 1540, 1550, 1170
			IF(Q) 1540, 1550, 1540
1520			CONTINUE
			IF(Q) 1610, 2010, 1610
1530			CONTINUE
C			IF(Q) 1170, 1550, 1540
			IF(Q) 1540, 1550, 1540
1540			CONTINUE
			PEN = DOWN
			GO TO 1170
1550			CONTINUE
			IF(J .EQ. 1) GO TO 1170
			DI = DDD
			TO = L - LD
			T = TO + DI
			P1 = Q
1560			IF(ABS(DI) .LT. END) GO TO 1570
			CALL CANSEE(A, XI, T, N, M, PO)
			DI = DI * 0.5
			IF(PO .EQ. 0) GO TO 1565
			TO = T
			P1 = PO
			T = T - DI
			GO TO 1560
1565			T = T + DI
			GO TO 1560
1570			CONTINUE
			T = TO
			IF(P1 * P) 1170, 1170, 1580
1580			CONTINUE
1590			CONTINUE
			ZP = A(I,L-LD) +(T-L+LD)*(A(I,L) -A(I,L-LD))/LD
			CALL ROTATE(XI, T, ZP, XP, HH, VV)
			HH = ((XP - QX) * DM - (HH - QY) * DL) * DD
			VV = (VV - QZ) * DD
			HH = (HH - H(9)) * SH
			VV = (VV - V(9)) * SV
			CALL PLOT(HH, VV, PEN)
C1600			PEN = 5 - PEN
1600			PEN = 1 - PEN
			GO TO 1170
1610			CONTINUE
			PEN = DOWN
			DI = DDD
			TO = L - LD
			T = TO + DI
			P1 = Q
1620			IF(ABS(DI) .LT. END) GO TO 1630
			    CALL CANSEE(A, XI, T, N, M, PO)
			    DI = DI * 0.5
			    IF(PO .EQ. 0) GO TO 1625
				TO = T
				P1 = PO
				T = T + DI
				GO TO 1620
1625			    T = T - DI
			    GO TO 1620
1630			CONTINUE
			T = TO
			IF(P1 * Q) 1600, 1600, 1590
1170			CALL ROTATE(XI, YJ, A(I,L), XP, HH, VV)
			HH = ((XP-QX)*DM - (HH-QY)*DL) * DD
			VV = (VV - QZ) * DD
1180			HH = (HH - H(9)) * SH
1190			VV = (VV - V(9)) * SV
			CALL PLOT(HH, VV, PEN)
2010			Q = P
2020		    CONTINUE
		    L = L + LD
		    LD = -LD
		    DDD = -DDD
2040		CONTINUE
2060	    CONTINUE
2130	CONTINUE
	END
C
	SUBROUTINE CANSEE(Z, XI, YJ, M, N, P)
	COMMON ANGA, ANGB, HV, D
	COMMON DL, DM, DN, CX, CY, CZ, QX, QY, QZ, DD
	INTEGER CUM, CNT, P, POLD
	REAL I, J, II, JJ
	DIMENSION Z(100, 100)
	IR = XI
	JC = YJ
	ZB = Z(IR, JC)
	IF(XI .EQ. IR) GO TO 2
	    ZB = Z(IR,JC) + (XI-IR) * (Z(IR+1,JC) - Z(IR,JC))
	GO TO 4
2	IF(YJ .EQ. JC) GO TO 4
	    ZB = Z(IR,JC) + (YJ-JC) * (Z(IR,JC+1) - Z(IR,JC))
4	CONTINUE
	XEND = 0.0
	DX = 0.0
	YMULT = 0.0
	ZMULT = 0.0
	IF(XI .EQ. CX) GO TO 10
	    YMULT = (YJ - CY) / (XI - CX)
	    ZMULT = (ZB - CZ) / (XI - CX)
	    DX = 1.0
	    XEND = M + 1
	    IF(XI .LT. CX) GO TO 10
		DX = -1.0
		XEND = 0.0
10	CONTINUE
	YEND = 0.0
	DY = 0.0
	XMULT = 0.0
	IF(YJ .EQ. CY) GO TO 20
	    XMULT = (XI - CX) / (YJ - CY)
	    IF(ZMULT .EQ. 0.0) ZMULT = (ZB - CZ) / (YJ - CY)
	    DY = 1.0
	    YEND = N + 1
	    IF(YJ .LT. CY) GO TO 20
		DY = -1.0
		YEND = 0.0
20	CONTINUE
	CUM = 0
	CNT = 0
	P = 0
	XB = XI
	YB = YJ
30	CONTINUE
	II = AINT(XB)
	JJ = AINT(YB)
	XSTEP = DX
	YSTEP = DY
	IF(XB .EQ. II) GO TO 40
	    IF(DX. LT. 0.0) XSTEP = 0.0
	    GOTO 45
40	IF(YB .EQ. JJ) GO TO 45
	    IF(DY .LT. 0.0) YSTEP = 0.0
45	CONTINUE
	I = II + XSTEP
	J = JJ + YSTEP
	IF(I .EQ. XEND) GO TO 81
	IF(J .EQ. YEND) GO TO 81
	    XB = CX + XMULT * (J - CY)
	    YB = CY + YMULT * (I - CX)
	    IF(DX .LT. 0.0) GO TO 55
		IF(XB .LT. I) GO TO 60
50		    XB = I
		    GO TO 65
55	    IF(XB .LT. I) GO TO 50
60		YB = J
65	    CONTINUE
	    ZB = CZ + ZMULT * (XB - CX)
	    IR = I
	    JC = J
	    IF(YB .NE. J) GO TO 70
		IDX = I - DX
		ZS = Z(IR, JC) - DX*(XB-I)*(Z(IDX,JC) - Z(IR,JC))
	    GO TO 75
70		JDY = J - DY
		ZS = Z(IR,JC) - DY*(YB-J)*(Z(IR,JDY) - Z(IR,JC))
75	    CONTINUE
	    SGN = 1
	    IF(ZB .LT. ZS) SGN = -1
	    CUM = CUM + SGN
	    CNT = CNT + 1
	    IF(IABS(CUM) .EQ. CNT) GO TO 30
	    GO TO 90
81	CONTINUE
	P = 1
	IF(CUM) 84, 86, 90
84	    P = -1
	GO TO 90
86	CONTINUE
	    IF(ZB .LE. CZ) GO TO 90
		P = -1
90	CONTINUE
	END
C
C DUMMY ROUTINES TO REPLACE PLOT FUNCTIONS
C
	SUBROUTINE PICSIZ(X, Y)
	REAL X, Y
	CHARACTER GFPL*8
	DATA GFPL/' * mcgW1'/
	IF(X.EQ.0.0) GOTO 3001
	IF(Y.NE.0.0) GOTO 3005
3001    WRITE(6, 3002)
3002    FORMAT(1X,'D')
	GOTO 3011
C
3005	WRITE(6,5000)X,Y
5000	FORMAT(1X,'PICSIZ: X DIMENSION = ',F10.5,'  Y DIMENSION = ',
+	F10.5)
	GFPL(1:1) = CHAR(27)
	GFPL(3:3) = CHAR(27)
	WRITE(6, 3010)GFPL
3010	FORMAT(1X,A)
3011	CONTINUE
	END
C
	SUBROUTINE PLOT(X, Y, PEN)
	COMMON ANGA, ANGB, HV, D
	COMMON DL, DM, DN, CX, CY, CZ, QX, QY, QZ, DD
	REAL X, Y
	REAL YSCALED
	INTEGER PEN
	INTEGER IX, IY, IXO, IYO
C	WRITE(6,5001)PEN, X, Y
	IX = IFIX(X * 639./HV)
	YSCALED = 50.*Y
	IF(YSCALED.LT.0.0) YSCALED = 0.0
	IF(YSCALED.GT.239.) YSCALED = 239.
	IY = IFIX(YSCALED)
	IF(PEN.LT.0) GOTO 3111
	IF(PEN.GT.1) GOTO 3111
	IF(PEN.NE.0) GOTO 3110
	    WRITE(6, 3100) IXO, IYO, IX, IY
3100	    FORMAT(1X,'M',I3.3,',',I3.3,';L',I3.3,',',I3.3,';')
3110	IXO = IX
	IYO = IY
	RETURN
C
3111	CONTINUE
	WRITE(6,5005)
5005	FORMAT(1X,'D')
	WRITE(6,5001)PEN, X, Y
5001	FORMAT(1X,'PLOT : Pen = ',I4,'  at location (',F10.5,F10.5,')')
	WRITE(6,5004)CHAR(27)
5004	FORMAT(1X,A1,'m')
	END
C
	SUBROUTINE ORIGIN(X,Y,PEN)
	REAL X, Y
	INTEGER PEN
C	WRITE(6,5003)
C	WRITE(6,5002)PEN, X, Y
C5002	FORMAT(1X,'ORIGIN : 3rd arg = ',I4,' 1st arg = ',F10.5,
C+	' 2nd arg = ',F10.5)
C5003	FORMAT(1X,'D')
	END
