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
 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.