C This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/perspec.for
C
      SUBROUTINE ROTATN(XDIST,DIST,TH,PH)
C==============F. Y. Borden - PENN STATE FORESTRY SCHOOL - MAY 1966.
C  ROTATN DEFINES PARAMETERS FOR USING THE CENTRAL PROJECTION
C  PERSPECTIVE ROUTINE PERSPC. THE EYE IS AT (XDIST,0,0)
C  AND THE LINE OF SIGHT IS FROM THE EYE TO (0,0,0).
C  THE PROJECTION PLANE IS PARALLEL TO THE Y-Z PLANE,
C  INTERSECTING THE X-AXIS AT DIST FROM THE EYE TOWARDS
C  THE ORIGIN.  BEFORE PROJECTION THE FIGURE IS ROTATED
C  BY TH(RADIANS) ABOUT THE Z-AXIS THEN BY PH(RADIANS)
C  ABOUT THE Y-AXIS.
C
C     PARAMETER DESCRIPTION:
C XDIST  IS GIVEN AS THE LOCATION OF THE EYE ON THE X-AXIS.
C DIST   IS GIVEN AS THE DISTANCE FROM THE DYE TO THE PROJECTION
C        PLANE, MEASURED FROM THE EYE TOWARDS THE ORGIN.
C TH     IS GIVEN AS THE ANGLE, MEASURED COUNTERCLOCKWISE IN RADIANS,
C        THROUGH WHICH THE 3-SPACE FIGURE IS ROTATED ABOUT THE Z-AXIS
C        BEFORE PROJECTION.  COUNTERCLOCKWISE ROTATION IN THIS CASE
C        CARRIES THE X-AXIS INTO THE Y-AXIS.
C PH     IS GIVEN AS THE ANGLE OF COUNTERCLOCKWISE ROTATION ABOUT THE
C        Y-AXIS. COUNTERCLOCKWISE IN THIS CASE CARRIES THE Z-AXIS INTO
C        THE X-AXIS.
C     USES: SIN,COS
      COMMON /HDK/ XX,YY,ZZ,XY,XZ,YX,ZX,ZY,SAVED,SAVEX
      YX=SIN(TH)
      YY=COS(TH)
      XZ=SIN(PH)
      ZZ=COS(PH)
      XX=YY*ZZ
      XY=-1.*YX*ZZ
      ZX=-1.*YY*XZ
      ZY=YX*XZ
      SAVED=DIST
      SAVEX=XDIST
      RETURN
      END
      SUBROUTINE PERSPC(X,Y,Z,XP,YP)
C=========F. YATES BORDEN - PENN STATE COMPUTATION CENTER - MAY 1968.
C  PERSPC TAKES A POINT (X,Y,Z) AND FINDS ITS TWO DIMENSIONAL IMAGE
C  (XP,YP) ON THE PROJECTION PLANE.  THE PARAMETERS RELATED
C  TO THIS TRANSFORMATION MUST HAVE BEEN DEFINED BY A PREVIOUS
C  CALL TO EITHER SETP OR ROTATN.
C
C     PARAMETER DESCRIPTION:
C X,Y,Z   IS GIVEN AS THE ORDER 3-TUPLE REPRESENTING THE 3-SPACE
C         CARTESIAN COORDINATES.
C XP,YP   IS RETURNED AS THE ORDERED 2-TUPLE REPRESENTING THE
C         CARTESIAN COORDINATES OF THE 2-SPACE CENTRAL PROJECTION.
C
      COMMON /HDK/ XX,YY,ZZ,XY,XZ,YX,ZX,ZY,SAVED,SAVEX
      X2=XX*X+XY*Y+XZ*Z
      Y2=YX*X+YY*Y
      Z2=ZX*X+ZY*Y+ZZ*Z
      FAC=SAVED/(SAVEX-X2)
      XP=Y2*FAC
      YP=Z2*FAC
      RETURN
      END
      SUBROUTINE SETP (XE,YE,ZE,DIST)
C=========F. YATES BORDEN - PENN STATE COMPUTATION CENTER - MAY 1968.
C  SETP  DEFINES PARAMETERS FOR USING THE CENTRAL PROJECTION
C  PERSPECTIVE ROUTINE PERSPC.  (XE,YE,ZE) IS THE LOCATION
C  OF THE EYE.  THE LINE OF SIGHT EXTENDS FROM THE EYE TO
C  (0,0,0).  THE PROJECTION PLANE IS PERPENDICULAR TO THE
C  LINE OF SIGHT, INTERSECTING IT AT DIST FROM THE EYE.
C
C     PARAMETER DESCRIPTION:
C XE,YE,ZE   IS GIVEN AS THE CARTESIAN COORDINATES OF THE EYE, THE
C            POINT FROM WHICH A 3-SPACE FIGURE IS TO BE VIEWED.
C DIST       IS GIVEN AS THE DISTANCE FROM THE EYE TO THE PROJECTION
C            PLANE.
C
C     USES: SQRT,ARSIN,ROTATN
      D1=XE*XE+YE*YE
      XDIST=SQRT(D1+ZE*ZE)
      D1=SQRT(D1)
      TH=-1.*ARSIN(YE/D1)
      PH=ARSIN(ZE/XDIST)
      IF(XE.GE.0) GO TO 1
      TH=3.14159-TH
    1 CALL  ROTATN(XDIST,DIST,TH,PH)
      RETURN
      END