home *** CD-ROM | disk | FTP | other *** search
- 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.
-