aboutsummaryrefslogtreecommitdiff
path: root/plotlib/plt_3D.f
diff options
context:
space:
mode:
Diffstat (limited to 'plotlib/plt_3D.f')
-rw-r--r--plotlib/plt_3D.f405
1 files changed, 405 insertions, 0 deletions
diff --git a/plotlib/plt_3D.f b/plotlib/plt_3D.f
new file mode 100644
index 0000000..0be4f82
--- /dev/null
+++ b/plotlib/plt_3D.f
@@ -0,0 +1,405 @@
+C***********************************************************************
+C Module: plt_3D.f
+C
+C Copyright (C) 1996 Harold Youngren, Mark Drela
+C
+C This library is free software; you can redistribute it and/or
+C modify it under the terms of the GNU Library General Public
+C License as published by the Free Software Foundation; either
+C version 2 of the License, or (at your option) any later version.
+C
+C This library is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+C Library General Public License for more details.
+C
+C You should have received a copy of the GNU Library General Public
+C License along with this library; if not, write to the Free
+C Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+C
+C Report problems to: guppy@maine.com
+C or drela@mit.edu
+C***********************************************************************
+C
+C***********************************************************************
+C --- Xplot11 3D routines
+C
+C Version 4.46 11/28/01
+C
+C Note: These are routine(s) that provide some means of displaying
+C 3D objects in conjunction with the usual XPlot11 routines.
+C They are by no means complete but can serve as a starting
+C point for doing simple 3D graphics.
+C
+C***********************************************************************
+
+
+ subroutine VIEW(X,Y,Z,N,XP,YP,XOB,YOB,ZOB,ROBINV,XUP,YUP,ZUP)
+ DIMENSION X(N), Y(N), Z(N)
+ DIMENSION XP(N), YP(N)
+C........................................................................
+C
+C Projects one or more points in 3-D Cartesian space
+C onto a 2-D plane from the viewpoint of an observer
+C at a specified location. This can be used to view
+C a 3-D object (described by a set of x,y,z points)
+C by projecting the points into a set of x,y points for
+C plotting on a planar 2-D graphics screen.
+C The viewing plane, which has its own x,y coordinate
+C system, always faces the observer but can be turned
+C around the viewing axis, thus simulating the observer
+C tilting his head while looking at the object. This tilt
+C is specified by giving a vector in x,y,z space which
+C "points up" relative to the observer.
+C The distance of the observer from the object is specified
+C explicitly. This does not affect much the size of the viewed
+C object, since the viewing plane contains the 3-D space origin
+C and hence is at or near the object. It does however affect the
+C apparent distortion of the object due to perspective. This
+C is very useful to convey the 3-dimensionality of the object.
+C If the observer is very very far away, there is no distortion
+C (as in a mechanical drawing).
+C
+C X,Y,Z Cartesian point coordinates (input)
+C N number of points (input)
+C XP,YP projected point coordinates on viewing plane (output)
+C XOB,YOB,ZOB Cartesian vector pointing towards observer (input)
+C (magnitude irrelevant)
+C ROBINV 1/(distance to observer) (input)
+C XUP,YUP,ZUP Cartesian vector which points "up" from the
+C observer's viewpoint (magnitude irrelevant) (input)
+C
+C Mark Drela July 1988
+C........................................................................
+C
+C---- unit view vector perpendicular to viewing plane (towards observer)
+ XOBN = XOB/SQRT(XOB**2 + YOB**2 + ZOB**2)
+ YOBN = YOB/SQRT(XOB**2 + YOB**2 + ZOB**2)
+ ZOBN = ZOB/SQRT(XOB**2 + YOB**2 + ZOB**2)
+C
+C---- vector along plane's local x coordinate: (up vector)X(view vector)
+ XIP = YUP*ZOBN - ZUP*YOBN
+ YIP = ZUP*XOBN - XUP*ZOBN
+ ZIP = XUP*YOBN - YUP*XOBN
+C
+C---- normalize plane's x coordinate vector
+ XIHAT = XIP/SQRT(XIP**2 + YIP**2 + ZIP**2)
+ YIHAT = YIP/SQRT(XIP**2 + YIP**2 + ZIP**2)
+ ZIHAT = ZIP/SQRT(XIP**2 + YIP**2 + ZIP**2)
+C
+C---- unit vector along plane's y coordinate: (view vector)X(x unit vector)
+ XJHAT = YOBN*ZIHAT - ZOBN*YIHAT
+ YJHAT = ZOBN*XIHAT - XOBN*ZIHAT
+ ZJHAT = XOBN*YIHAT - YOBN*XIHAT
+C
+C---- go over all points
+ DO 10 I=1, N
+C
+ RDOTR = X(I)*XOBN + Y(I)*YOBN + Z(I)*ZOBN
+C
+C------ viewing-axis component of vector
+ DRX = RDOTR*XOBN
+ DRY = RDOTR*YOBN
+ DRZ = RDOTR*ZOBN
+C
+C------ projected vector scaling factor due to perspective
+ VSCAL = 1.0 / SQRT( (XOBN-ROBINV*DRX)**2
+ & + (YOBN-ROBINV*DRY)**2
+ & + (ZOBN-ROBINV*DRZ)**2 )
+C
+C------ dot vector into plane coordinate system unit vectors, and scale
+ XP(I) = (XIHAT*X(I) + YIHAT*Y(I) + ZIHAT*Z(I))*VSCAL
+ YP(I) = (XJHAT*X(I) + YJHAT*Y(I) + ZJHAT*Z(I))*VSCAL
+C
+ 10 CONTINUE
+C
+ RETURN
+ END
+
+
+ subroutine VIEWR(R,N,RP,ROB,ROBINV,RUP)
+ DIMENSION R(3,N)
+ DIMENSION RP(3,N)
+ DIMENSION ROB(3), RUP(3)
+C........................................................................
+C
+C Same as VIEW, but vectors are passed in R(1:3,...) array form.
+C The out-of-plane RP(3...) value is also returned.
+C
+C R(..) Cartesian point coordinates (input)
+C N number of points (input)
+C RP(..) projected point coordinates on viewing plane (output)
+C ROB(.) Cartesian vector pointing towards observer (input)
+C (magnitude irrelevant)
+C ROBINV 1/(distance to observer) (input)
+C RUP(.) Cartesian vector which points "up" from the
+C observer's viewpoint (magnitude irrelevant) (input)
+C
+C Mark Drela July 2007
+C........................................................................
+ REAL RIP(3), RIN(3), RJN(3), RKN(3), DR(3)
+C
+C---- unit view vector perpendicular to viewing plane (towards observer)
+ SOB = SQRT(ROB(1)**2 + ROB(2)**2 + ROB(3)**2)
+ RKN(1) = ROB(1)/SOB
+ RKN(2) = ROB(2)/SOB
+ RKN(3) = ROB(3)/SOB
+C
+C---- vector along plane's local x coordinate: (up vector)X(view vector)
+ RIP(1) = RUP(2)*RKN(3) - RUP(3)*RKN(2)
+ RIP(2) = RUP(3)*RKN(1) - RUP(1)*RKN(3)
+ RIP(3) = RUP(1)*RKN(2) - RUP(2)*RKN(1)
+C
+C---- normalize plane's x coordinate vector
+ SIP = SQRT(RIP(1)**2 + RIP(2)**2 + RIP(3)**2)
+ RIN(1) = RIP(1)/SIP
+ RIN(2) = RIP(2)/SIP
+ RIN(3) = RIP(3)/SIP
+C
+C---- unit vector along plane's y coordinate: (view vector)X(x unit vector)
+ RJN(1) = RKN(2)*RIN(3) - RKN(3)*RIN(2)
+ RJN(2) = RKN(3)*RIN(1) - RKN(1)*RIN(3)
+ RJN(3) = RKN(1)*RIN(2) - RKN(2)*RIN(1)
+C
+C---- go over all points
+ DO 10 I=1, N
+ RDOTR = R(1,I)*RKN(1) + R(2,I)*RKN(2) + R(3,I)*RKN(3)
+C
+C------ viewing-axis component of vector
+ DR(1) = RDOTR*RKN(1)
+ DR(2) = RDOTR*RKN(2)
+ DR(3) = RDOTR*RKN(3)
+C
+C------ projected vector scaling factor due to perspective
+ VSCAL = 1.0 / SQRT( (RKN(1)-ROBINV*DR(1))**2
+ & + (RKN(2)-ROBINV*DR(2))**2
+ & + (RKN(3)-ROBINV*DR(3))**2 )
+C
+C------ dot vector into plane coordinate system unit vectors, and scale
+ RP(1,I) = (RIN(1)*R(1,I) + RIN(2)*R(2,I) + RIN(3)*R(3,I))*VSCAL
+ RP(2,I) = (RJN(1)*R(1,I) + RJN(2)*R(2,I) + RJN(3)*R(3,I))*VSCAL
+ RP(3,I) = (RKN(1)*R(1,I) + RKN(2)*R(2,I) + RKN(3)*R(3,I))*VSCAL
+ 10 CONTINUE
+C
+ RETURN
+ END ! VIEWR
+
+
+
+ SUBROUTINE PROJMATRIX3 (ROTZ,ROTY,RMAT)
+C...Purpose: To define rotation and transformation matrix. The input
+C pair of angles ROTZ,ROTY specify the viewpoint by
+C an angle about the Z axis (CCW) and an angle about
+C the newly rotated Y axis (CCW). Both angles are
+C right-handed in a conventional sense about each axis.
+C
+C...Input: ROTZ rotation of viewpoint about Z axis (deg)
+C ROTY rotation of viewpoint about new Y axis (deg)
+C
+C...Output: RMAT 3x3 rotation and perspective matrix
+C
+ DIMENSION RMAT(3,3)
+C
+ DTR = 4.0*ATAN(1.0)/180.0
+ COSZ = COS(ROTZ*DTR)
+ SINZ = SIN(ROTZ*DTR)
+ COSY = COS(ROTY*DTR)
+ SINY = SIN(ROTY*DTR)
+C
+C---Rotation matrix (rotation about Z, then rotation about Y)
+c xx = -( SINZ*X + COSZ*Y)
+ RMAT(1,1) = -SINZ
+ RMAT(2,1) = -COSZ
+ RMAT(3,1) = 0.0
+C yy = SINY*COSZ*X - SINY*SINZ*Y + COSY*Z
+ RMAT(1,2) = SINY*COSZ
+ RMAT(2,2) = -SINY*SINZ
+ RMAT(3,2) = COSY
+c zz = -(COSY*COSZ*X - COSY*SINZ*Y - SINY*Z)
+ RMAT(1,3) = -COSY*COSZ
+ RMAT(2,3) = COSY*SINZ
+ RMAT(3,3) = SINY
+C
+c xx = -( SINZ*X + COSZ*Y)
+c yy = SINY*COSZ*X - SINY*SINZ*Y + COSY*Z
+c zz = -(COSY*COSZ*X - COSY*SINZ*Y - SINY*Z)
+C
+c write(*,*) 'Rmatrix row1 ', (RMAT(1,L),L=1,3)
+c write(*,*) 'Rmatrix row2 ', (RMAT(2,L),L=1,3)
+c write(*,*) 'Rmatrix row3 ', (RMAT(3,L),L=1,3)
+c read(*,*)
+C
+ RETURN
+ END
+
+
+ SUBROUTINE PROJMATRIX4 (ROTZ,ROTY,RDIST,RMAT)
+C...Purpose: To define rotation and perspective matrix. The input
+C pair of angles ROTZ,ROTY specify the viewpoint by
+C an angle about the Z axis (CCW) and an angle about
+C the newly rotated Y axis (CCW). Both angles are
+C right-handed in a conventional sense about each axis.
+C The observer distance RDIST specifies the distance from
+C the origin to the observer along the viewpoint direction.
+C
+C...Input: ROTZ rotation of viewpoint about Z axis (deg)
+C ROTY rotation of viewpoint about new Y axis (deg)
+C RDIST distance from origin to observer along viewpoint
+C
+C...Output: RMAT 4x4 rotation and perspective matrix
+C
+ DIMENSION AMAT(4,4),PMAT(4,4), RMAT(4,4)
+C
+ DTR = 4.0*ATAN(1.0)/180.0
+ COSZ = COS(ROTZ*DTR)
+ SINZ = SIN(ROTZ*DTR)
+ COSY = COS(ROTY*DTR)
+ SINY = SIN(ROTY*DTR)
+C
+C---Rotation matrix (rotation about Z, then rotation about Y)
+c xx = -( SINZ*X + COSZ*Y)
+ AMAT(1,1) = -SINZ
+ AMAT(2,1) = -COSZ
+ AMAT(3,1) = 0.0
+ AMAT(4,1) = 0.0
+C yy = SINY*COSZ*X - SINY*SINZ*Y + COSY*Z
+ AMAT(1,2) = SINY*COSZ
+ AMAT(2,2) = -SINY*SINZ
+ AMAT(3,2) = COSY
+ AMAT(4,2) = 0.0
+c zz = -(COSY*COSZ*X - COSY*SINZ*Y - SINY*Z)
+ AMAT(1,3) = -COSY*COSZ
+ AMAT(2,3) = COSY*SINZ
+ AMAT(3,3) = SINY
+ AMAT(4,3) = 0.0
+C
+ AMAT(1,4) = 0.0
+ AMAT(2,4) = 0.0
+ AMAT(3,4) = 0.0
+ AMAT(4,4) = 1.0
+C
+c xx = -( SINZ*X + COSZ*Y)
+c yy = SINY*COSZ*X - SINY*SINZ*Y + COSY*Z
+c zz = -(COSY*COSZ*X - COSY*SINZ*Y - SINY*Z)
+c
+C---Perspective matrix with projection on Z plane
+ PMAT(1,1) = 1.0
+ PMAT(2,1) = 0.0
+ PMAT(3,1) = 0.0
+ PMAT(4,1) = 0.0
+C
+ PMAT(1,2) = 0.0
+ PMAT(2,2) = 1.0
+ PMAT(3,2) = 0.0
+ PMAT(4,2) = 0.0
+C
+ PMAT(1,3) = 0.0
+ PMAT(2,3) = 0.0
+ PMAT(3,3) = 1.0
+ PMAT(4,3) = 0.0
+C
+ PMAT(1,4) = 0.0
+ PMAT(2,4) = 0.0
+ PMAT(3,4) = -RDIST
+ PMAT(4,4) = 1.0
+C
+C---Product of matrices is perspective matrix
+ DO J=1, 4
+ DO K=1, 4
+ TMP = 0.0
+ DO L=1, 4
+ TMP = TMP + AMAT(J,L)*PMAT(L,K)
+ END DO
+ RMAT(J,K) = TMP
+ END DO
+ END DO
+C
+c write(*,*) 'Rmatrix row1 ', (RMAT(1,L),L=1,4)
+c write(*,*) 'Rmatrix row2 ', (RMAT(2,L),L=1,4)
+c write(*,*) 'Rmatrix row3 ', (RMAT(3,L),L=1,4)
+c write(*,*) 'Rmatrix row4 ', (RMAT(4,L),L=1,4)
+c read(*,*)
+C
+ RETURN
+ END
+
+
+ SUBROUTINE ROTPTS3 (RMAT,PTS_IN,NPTS,PTS_OUT)
+C...Purpose: To rotate array of points to a new viewpoint by
+C parallel projection. The input rotation matrix
+C contains the transformation data in a 3x3 matrix.
+C
+C...Input: RMAT 3x3 transformation matrix
+C PTS_IN array (3xNPTS) of input points
+C NPTS number of points in arrays
+C
+C...Output: PTS_OUT array (3xNPTS) of transformed points
+C
+ DIMENSION PTS_IN(3,NPTS), PTS_OUT(3,NPTS)
+ DIMENSION RMAT(4,4)
+C
+ DO I = 1, NPTS
+C
+ DO J=1, 3
+ TMP = 0.0
+ DO K=1, 3
+ TMP = TMP + PTS_IN(K,I)*RMAT(K,J)
+ END DO
+ PTS_OUT(J,I) = TMP
+ END DO
+C
+ END DO
+C
+ RETURN
+ END
+
+
+ SUBROUTINE ROTPTS4 (RMAT,PTS_IN,NPTS,PTS_OUT)
+C...Purpose: To rotate array of points to a new viewpoint by
+C perspective projection. The input rotation matrix
+C contains the transformation and perspective data in
+C a 4x4 matrix in homogeneous coordinates. Note that input
+C coordinates may need to be z-clipped if the user trans.
+C moves the points behind the observer. A check is made
+C for a singular perspective point (at observer).
+C
+C...Input: RMAT 4x4 rotation and perspective matrix
+C PTS_IN array (3xNPTS) of input points
+C NPTS number of points in arrays
+C
+C...Output: PTS_OUT array (3xNPTS) of transformed points
+C
+C...Note: You may need to translate your points to recenter them
+C about the origin to get good perspective views. Points
+C off to the side get pretty distorted...
+C
+ DIMENSION PTS_IN(3,NPTS), PTS_OUT(3,NPTS)
+ DIMENSION RMAT(4,4), PTMP(4), PPTMP(4)
+C
+ DO I = 1, NPTS
+C
+ PTMP(1) = PTS_IN(1,I)
+ PTMP(2) = PTS_IN(2,I)
+ PTMP(3) = PTS_IN(3,I)
+ PTMP(4) = 1.0
+C
+ DO J=1, 4
+ TMP = 0.0
+ DO K=1, 4
+ TMP = TMP + PTMP(K)*RMAT(K,J)
+ END DO
+ PPTMP(J) = TMP
+ END DO
+C
+ IF(PPTMP(4).NE.0.0) THEN
+ PTS_OUT(1,I) = PPTMP(1)/PPTMP(4)
+ PTS_OUT(2,I) = PPTMP(2)/PPTMP(4)
+ PTS_OUT(3,I) = PPTMP(3)/PPTMP(4)
+ ELSE
+ WRITE(*,*) 'Homogeneous coordinate singular for pt ',I
+ ENDIF
+C
+ END DO
+C
+ RETURN
+ END