home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Simtel MSDOS 1992 September
/
Simtel20_Sept92.cdr
/
msdos
/
plot
/
3d_plot.for
< prev
next >
Wrap
Text File
|
1992-09-08
|
9KB
|
269 lines
Date: Tue, 29 Dec 87 13:50:53 EST
From: Larry Granroth 319/335-1960 <"IOWASP::GRANROTH"@nssdca.GSFC.NASA.GOV>
Subject: fortran code for 3d surface plotting
Info-IBMPC Digest Volume 6, Issue 72, contained a request from
<SIMPSON_P%MERCURY.ceo.dg.com@adam.dg.com> for a public domain 3-D
plotting program. I wrote a quick-and-dirty program which uses
Tektronix PLOT10/TCS calls and ran on a VAX several years ago. The
algorithm effectively plots "slices" through a surface, so it isn't
perfect, but it is faster than more sophisticated "web" algorithms.
It should be fairly easy to convert to any FORTRAN system and
substitute different plot drivers, so if anyone is interested, here is
the source code:
PROGRAM TEK3D
C----------------------------------------------------------------------
C
C DEMONSTRATION PROGRAM TO DRIVE TEK 4014 TYPE TERMINAL
C MUST BE LINKED WITH PLOT10/TCS TYPE LIBRARY.
C DRAWS PERSPECTIVE SURFACE OF SINC(R) ON 40X40 GRID
C WITH R=SQRT(X*X+Y*Y)/2.
C
C----------------------------------------------------------------------
DIMENSION Z(40,40)
COMMON /BLK3D/ IXCTR,IZCTR,IWIDE,IDEEP,IHIGH,VDIST
IXCTR=512
IZCTR=390
IWIDE=512
IDEEP=512
IHIGH=480
VDIST=2048.
C
H 67 528 528 0 0 528 0 1 1 0H DO 20 J=1,40
Y=FLOAT(J)-20.5
YY=Y*Y
DO 10 I=1,40
X=FLOAT(I)-20.5
R=SQRT(X*X+YY)/2.
Z(I,J)=SIN(R)/R
10 CONTINUE
20 CONTINUE
C
99 TYPE '('' ENTER AZ,EL, AND IFLAG: '',$)'
ACCEPT *,AZ,EL,IFLAG
IF (IFLAG.LT.0) GOTO 9999
CALL PLT3D(Z,40,40,-.8,.8,AZ,EL,IFLAG)
PAUSE
GOTO 99
C
9999 STOP
END
C----------------------------------------------------------------------
C
C RULED SURFACE VERSION OF PLT3D WRITTEN BY LARRY GRANROTH
C
C TEKTRONIX 4014 VERSION USES PLOT10/TCS ROUTINES
C
C VAX FORTRAN VERSION
C
C PARAMETERS:
C DATA 2-DIM REAL DATA ARRAY
C NX,NY X AND Y INTEGER DIMENSIONS OF DATA ARRAY
C DMIN DATA VALUE TO BE SCALED TO BOTTOM OF 3-D PLOT CUBE
C DMAX VALUE TO BE SCALED TO TOP OF PLOT CUBE
C AZ ANGLE (DEG) TO ROTATE PLOT CUBE CCW
C EL ANGLE TO TILT PLOT CUBE TOWARD VIEWER
C IFLAG 1=X-SCAN ONLY, 2=Y-SCAN ONLY, OTHER=BOTH
C
C COMMON BLOCK: /BLK3D/
C IXCTR HORIZONTAL PIXEL COORDINATE TO CENTER PLOT CUBE
C IZCTR VERTICLE COORDINATE (16 BIT INTEGERS!)
C IWIDE PLOT CUBE WIDTH IN PIXELS
C IDEEP PLOT CUBE DEPTH
C IHIGH PLOT CUBE HIGHT
C VDIST VIEWING DISTANCE IN PIXELS (FLOATING POINT)
C
C MODIFICATIONS:
C 1-28-82 X OR Y SCAN OPTION AND TOTAL 360 AZ -LJG
C
C 8-26-82 PERSPECTIVE FORSHORTENING ADDED
C VDIST ADDED TO COMMON BLOCK -LJG
C
C This source code is placed in the public domain by
C Larry Granroth as of 29 December 1987. I claim no
C responsibility for any problems associated with this
C code.
C
C----------------------------------------------------------------------
SUBROUTINE PLT3D(DATA,NX,NY,DMIN,DMAX,AZ,EL,IFLAG)
INTEGER MAXHID(1024), MINHID(1024)
LOGICAL PEN,NOHID,ABOV,BELO
DIMENSION DATA(NX,NY)
COMMON /BLK3D/ IXCTR,IZCTR,IWIDE,IDEEP,IHIGH,VDIST
COMMON /ROTRBK/ SINAZ,COSAZ,SINEL,COSEL
C
CALL INITT (1200)
C
XSCL=FLOAT(IWIDE)/FLOAT(NX)
H 67 528 528 0 0 528 0 1 1 0H YSCL=FLOAT(IDEEP)/FLOAT(NY)
IF (DMIN.EQ.DMAX) DMAX=DMIN+1.
ZSCL=FLOAT(IHIGH)/(DMAX-DMIN)
AZ=AZ*.01745329
EL=EL*.01745329
COSAZ=COS(AZ)
SINAZ=SIN(AZ)
COSEL=COS(EL)
SINEL=SIN(EL)
X0=(FLOAT(NX)-1.)/2.*XSCL
Y0=(FLOAT(NY)-1.)/2.*YSCL
D0=-(DMIN+(DMAX-DMIN)/2.)
IF (IFLAG.EQ.2) GOTO 5
IF (COSAZ.LT.0.) GOTO 3
LOY1=1
LOY2=NY
IEX=1
LIX1=2
LIX2=NX
YST1=-Y0
XSC1=-X0
YSTI=YSCL
XSCI=XSCL
GO TO 4
3 LOY1=-NY
LOY2=-1
IEX=NX
LIX1=1-NX
LIX2=-1
YST1=Y0
XSC1=X0
YSTI=-YSCL
XSCI=-XSCL
4 IF (IFLAG.EQ.1) GOTO 9
5 IF (SINAZ.LT.0.) GOTO 6
LOX1=1
LOX2=NX
IEY=NY
LIY1=1-NY
LIY2=-1
XST1=-X0
YSC1=Y0
XSTI=XSCL
YSCI=-YSCL
GO TO 8
6 LOX1=-NX
LOX2=-1
IEY=1
LIY1=2
LIY2=NY
XST1=X0
YSC1=-Y0
XSTI=-XSCL
YSCI=YSCL
8 IF (IFLAG.EQ.2) GOTO 100
9 DO 10 I=1,1024
MINHID(I)=959
10 MAXHID(I)=-960
Y1=YST1-YSTI
DO 50 LJ=LOY1,LOY2
J=IABS(LJ)
X1=XSC1
Y1=Y1+YSTI
Z1=(DATA(IEX,J)+D0)*ZSCL
CALL ROTR(X1,Y1,Z1,IX2,IZ2)
PEN=.FALSE.
H 67 528 528 0 0 528 0 1 1 0H IF (IZ2.GT.MAXHID(IX2).OR.IZ2.LT.MINHID(IX2)) PEN=.TRUE.
CALL MOVABS (IX2,IZ2) ! MOVE PEN TO BEGINNING OF SCAN LINE
DO 40 LI=LIX1,LIX2
I=IABS(LI)
X1=X1+XSCI
IX1=IX2
IZ1=IZ2
Z1=(DATA(I,J)+D0)*ZSCL
CALL ROTR(X1,Y1,Z1,IX2,IZ2)
NOHID=.TRUE.
IF (IX2.EQ.IX1) GOTO 40
SLOPE=FLOAT(IZ2-IZ1)/FLOAT(IX2-IX1)
XX=0.
IZ=IZ1
IX=IX1+1
DO 30 K=IX,IX2
XX=XX+1.
LAST=IZ1
IZ1=INT(SLOPE*XX+.5)+IZ
ABOV=(IZ1.GT.MAXHID(K))
BELO=(IZ1.LT.MINHID(K))
IF (ABOV.OR.BELO) GOTO 25
NOHID=.FALSE.
IF (PEN) CALL DRWABS(K-1,LAST) ! DRAW TO LAST UNHIDDEN POINT
PEN=.FALSE.
GOTO 30
25 IF (ABOV) MAXHID(K)=IZ1
IF (BELO) MINHID(K)=IZ1
IF (PEN) GOTO 30
CALL MOVABS(K,IZ1) ! MOVE PEN
PEN=.TRUE.
30 CONTINUE
40 IF (NOHID) CALL DRWABS(IX2,IZ2)
50 CONTINUE
IF (IFLAG.EQ.1) GOTO 150
100 DO 110 I=1,1024
MINHID(I)=959
110 MAXHID(I)=-960
X1=XST1-XSTI
DO 150 LJ=LOX1,LOX2
J=IABS(LJ)
Y1=YSC1
X1=X1+XSTI
Z1=(DATA(J,IEY)+D0)*ZSCL
CALL ROTR(X1,Y1,Z1,IX2,IZ2)
PEN=.FALSE.
IF (IZ2.GT.MAXHID(IX2).OR.IZ2.LT.MINHID(IX2)) PEN=.TRUE.
CALL MOVABS(IX2,IZ2) ! MOVE PEN TO BEGINNING OF SCAN LINE
DO 140 LI=LIY1,LIY2
I=IABS(LI)
Y1=Y1+YSCI
IX1=IX2
IZ1=IZ2
Z1=(DATA(J,I)+D0)*ZSCL
CALL ROTR(X1,Y1,Z1,IX2,IZ2)
NOHID=.TRUE.
IF (IX2.EQ.IX1) GOTO 140
SLOPE=FLOAT(IZ2-IZ1)/FLOAT(IX2-IX1)
XX=0.
IZ=IZ1
IX=IX1+1
DO 130 K=IX,IX2
XX=XX+1.
LAST=IZ1
IZ1=INT(SLOPE*XX+.5)+IZ
ABOV=(IZ1.GT.MAXHID(K))
H 67 528 528 0 0 528 0 1 1 0H BELO=(IZ1.LT.MINHID(K))
IF (ABOV.OR.BELO) GOTO 125
NOHID=.FALSE.
IF (PEN) CALL DRWABS(K-1,LAST) ! DRAW TO LAST UNHIDDEN POINT
PEN=.FALSE.
GOTO 130
125 IF (ABOV) MAXHID(K)=IZ1
IF (BELO) MINHID(K)=IZ1
IF (PEN) GOTO 130
CALL MOVABS(K,IZ1) ! MOVE PEN
PEN=.TRUE.
130 CONTINUE
140 IF (NOHID) CALL DRWABS(IX2,IZ2)
150 CONTINUE
CALL FINITT(0,780)
RETURN
END
C----------------------------------------------------------------------
SUBROUTINE ROTR(X1,Y1,Z1,IP,KP)
COMMON /BLK3D/ IXCTR,IZCTR,IWIDE,IDEEP,IHIGH,VDIST
COMMON /ROTRBK/ SINAZ,COSAZ,SINEL,COSEL
XROT =X1*COSAZ-Y1*SINAZ
YROT =X1*SINAZ+Y1*COSAZ
ZTILT=YROT*SINEL+Z1*COSEL
PERSP=VDIST/(VDIST+(YROT*COSEL-Z1*SINEL))
IP =INT(XROT *PERSP+0.5)+IXCTR
KP =INT(ZTILT*PERSP+0.5)+IZCTR
RETURN
END
Various ways to reach me include:
SPAN IOWASP::GRANROTH
ARPAnet GRANROTH%IOWASP.SPAN@NSSDC.GSFC.NASA.GOV
or @STAR.STANFORD.EDU
or @VLSI.JPL.NASA.GOV
or @SDS
BITNET CCOLJGPG@UIAMVS
TELEMAIL LGRANROTH/J.P.L.