*     This program demonstrates animation and 3D geometry in PGPLOT. It 
*     requires a fast, interactive display, e.g., /XWIN. Do not
*     specify a hardcopy device. The speed of the animation is limited by
*     the cpu speed of the host computer.
*     
*     Thanks to Dr Martin Weisser:
*     Date: Sun, 18 May 1997 16:14:01 CET
*     From: weisser@chclu.chemie.uni-konstanz.de


      PROGRAM PGDEM17
C-----------------------------------------------------------------------
C     Demonstration program for PGPLOT.
C-----------------------------------------------------------------------
      INTEGER PGOPEN
C
      WRITE (*,*) 'PGPLOT: Demonstration of animation and 3D geometry'
      WRITE (*,*) 'Select a fast, interactive device, e.g., /XWINDOW'
      IF (PGOPEN('?') .LE. 0) STOP
      CALL POLY3D
      CALL PGCLOS
C-----------------------------------------------------------------------
      END

      SUBROUTINE POLY3D
C     
      INTEGER NFRAMS
C     
      INTEGER NTOT, NLIN, IPOS, IFIRST          
      REAL T, T1, T2, T3, PI, W, W1, TET, TET1, ROT, ROT1
C     
      PARAMETER (NTOT=34)
      PARAMETER (T=1.618)
      PARAMETER (T1=1.0+T)
      PARAMETER (T2=-T)
      PARAMETER (T3=-T1)
      PARAMETER (W=0.60*T)
      PARAMETER (W1=-W)
      PARAMETER (TET=0.37)
      PARAMETER (TET1=-TET)
      PARAMETER (ROT=0.13)
      PARAMETER (ROT1=-ROT)
      PARAMETER (NLIN=49)
C     
      INTEGER I, J, L, III, ILINE, NTOTM6
      INTEGER ICDFOR, ICCFOR, ICTFOR, ICLFOR
      INTEGER ICDBCK, ICCBCK, ICTBCK, ICLBCK
      INTEGER ITYPE(NTOT), IARRAY(NLIN), JARRAY(NLIN), LITYPE(NLIN)
      REAL RQ, ZZ
      REAL THAXI1, PHAXI1, ALFA1, THAXI2, PHAXI2, ALFA2
      REAL THAXI3, PHAXI3, ALFA3, THAXI4, PHAXI4, ALFA4 
      REAL XOFF, YOFF, ZOFF
      REAL XARRAY(NTOT), YARRAY(NTOT), ZARRAY(NTOT), DISTAN(NLIN)
      REAL POLYS(3,NTOT), X(2), Y(2), C(3), CROT(3), RPOL(3,3)
      PARAMETER (PI=3.14159265359)
C     
C     Cartesian coordinates of the polygons 
C     
      DATA POLYS/ T, T, T,       T, T,T2,
     D     T,T2, T,       T,T2,T2,
     D     T2, T, T,      T2, T,T2,
     D     T2,T2, T,      T2,T2,T2,
     D     T1,1.0,0.0,    T1,-1.0,0.0,
     D     T3,1.0,0.0,    T3,-1.0,0.0,
     D     0.0,T1,1.0,    0.0,T1,-1.0,
     D     0.0,T3,1.0,    0.0,T3,-1.0,
     D     1.0,0.0,T1,    -1.0,0.0,T1,
     D     1.0,0.0,T3,   -1.0,0.0,T3, 
     C     W,    W,    W,     W,    W,   W1,
     C     W,   W1,    W,     W,   W1,   W1,
     C     W1,    W,    W,    W1,   W1,    W,
     C     W1,    W,   W1,    W1,   W1,   W1,                
     T     TET,  TET,  TET, TET1, TET1,  TET,
     T     TET1,  TET, TET1,  TET, TET1, TET1, 
     L     ROT,  0.0,  0.0, ROT1,  0.0,  0.0/
C     
      DATA ITYPE/1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
     C     2,2,2,2,2,2,2,2,
     T     3,3,3,3,
     L     4,4/
C     
C     Initialize the plot (no labels).
C     
      CALL PGENV(-3.2,3.2,-3.2,3.2,1,-2)
C     
C     Switch from page to page without typing return.  
C     
      CALL PGASK(.FALSE.)
C     
C     Rotation axis of the polygons 
C     
      THAXI1 = PI/4.0
      PHAXI1 = PI/4.5
      ALFA1  = 0.0
      THAXI2 = PI/6.0
      PHAXI2 = 0.0
      ALFA2  = 0.02
      THAXI3 = PI/2.0
      PHAXI3 = -PI/3.0
      ALFA3  = 0.0
      THAXI4 = -0.03
      PHAXI4 = PI/7.0
      ALFA4  = 0.9
C     
      XOFF=0.0
      YOFF=0.0
      ZOFF=0.0
      NTOTM6=NTOT-6
C     
C     Colors  
C     
      ICDFOR = 3
      ICDBCK = 10
C     
      ICCFOR = 8
      ICCBCK = 2
C     
      ICTFOR = 5
      ICTBCK = 4
C     
      ICLFOR = 1
      ICLBCK = 7
C     
      IPOS=1
C     
      WRITE(*,*)' Rotation with increasing velocity'
C     
      NFRAMS=3500 
C     
      DO 12 I=1,NTOT
         XARRAY(I) = POLYS(1,I)
         YARRAY(I) = POLYS(2,I)
         ZARRAY(I) = POLYS(3,I)
 12   CONTINUE
C     
      DO 30 L=1,NFRAMS
C     
         CALL PGBBUF
         CALL PGERAS
C     
         CALL SORTPP(NTOT,ITYPE,ZARRAY,YARRAY,XARRAY)
C     
         IFIRST=0
         DO 13 I=1,NTOT
            IF((ZARRAY(I).GE.0.0).AND.(IFIRST.EQ.0)) THEN
               IFIRST = 1
               IPOS = I
            END IF
 13      CONTINUE
C     
         IF(L.EQ.2800) CALL OFFSET (XOFF,YOFF,ZOFF)
         IF (MOD(L,500).EQ.0) THEN 
            CALL CHNAX(THAXI3,PHAXI3,THAXI2,PHAXI2,THAXI1,PHAXI1)
         END IF
C     
         DO 33 I=1,IPOS-1
            IF (ITYPE(I).EQ.1) THEN 
               CALL PGSCI(ICDBCK)
               CALL PGSLW(18)
            ELSE IF (ITYPE(I).EQ.2) THEN 
               CALL PGSCI(ICCBCK)
               CALL PGSLW(17)
            ELSE IF (ITYPE(I).EQ.3) THEN
               CALL PGSCI(ICTBCK)
               CALL PGSLW(15)
            ELSE 
               CALL PGSCI(ICLBCK)
               CALL PGSLW(14)
            END IF
            ZZ = ZARRAY(I)
            CALL PGPT(1,XARRAY(I)+0.2*ZZ,YARRAY(I)+0.3*ZZ,9)
 33      CONTINUE
C     
         DO 44 I=IPOS,NTOT
            IF (ITYPE(I).EQ.1) THEN 
               CALL PGSCI(ICDFOR)
               CALL PGSLW(18)
            ELSE IF (ITYPE(I).EQ.2) THEN 
               CALL PGSCI(ICCFOR)
               CALL PGSLW(17)
            ELSE IF (ITYPE(I).EQ.3) THEN
               CALL PGSCI(ICTFOR)
               CALL PGSLW(15)
            ELSE 
               CALL PGSCI(ICLFOR)
               CALL PGSLW(14)
            END IF
            ZZ = ZARRAY(I)
            CALL PGPT(1,XARRAY(I)+0.2*ZZ,YARRAY(I)+0.3*ZZ,9)
 44      CONTINUE
C     
         ILINE=0
C     
         DO 2000 I=2,NTOT
            DO 1000 J=1,I-1
               IF (ITYPE(I).EQ.ITYPE(J)) THEN
                  RQ = 0.0
                  RQ = RQ + ( (XARRAY(I)-XARRAY(J))**2+
     #                 (YARRAY(I)-YARRAY(J))**2+
     #                 (ZARRAY(I)-ZARRAY(J))**2  )
C     
                  IF ( ((RQ-0.0676)  .LT.0.001).OR.
     #                 ((RQ-1.095199).LT.0.001).OR.          
     #                 ((RQ-3.769809).LT.0.001).OR.
     #                 ((RQ-4.000000).LT.0.001)     ) THEN 
                     ILINE = ILINE + 1
                     DISTAN(ILINE) = ZARRAY(I)+ZARRAY(J)
                     IF(DISTAN(ILINE).LT.0.0) THEN 
                        LITYPE(ILINE) = -ITYPE(I)
                     ELSE  
                        LITYPE(ILINE) =  ITYPE(I)
                     END IF 
                     IARRAY(ILINE) = I
                     JARRAY(ILINE) = J
                  END IF
               END IF
 1000       CONTINUE
 2000    CONTINUE
C     
         CALL SORTLI(ILINE,DISTAN,IARRAY,JARRAY,LITYPE)
C     
         DO 3000 III=1,ILINE
            I=IARRAY(III)
            J=JARRAY(III)               
            ZZ = ZARRAY(I)
            X(1) = XARRAY(I)+0.2*ZZ
            Y(1) = YARRAY(I)+0.3*ZZ
            ZZ = ZARRAY(J)
            X(2) = XARRAY(J)+0.2*ZZ
            Y(2) = YARRAY(J)+0.3*ZZ
            IF (LITYPE(III).GT.0) THEN
               IF(LITYPE(III).EQ.1) THEN  
                  CALL PGSLW(10)
                  CALL PGSCI(ICDFOR)
               ELSE IF (LITYPE(III).EQ.2) THEN          
                  CALL PGSLW(8)
                  CALL PGSCI(ICCFOR)
               ELSE IF (LITYPE(III).EQ.3) THEN
                  CALL PGSLW(6)
                  CALL PGSCI(ICTFOR)
               ELSE
                  CALL PGSLW(4)
                  CALL PGSCI(ICLFOR)
               END IF
            ELSE
               IF(LITYPE(III).EQ.-1) THEN  
                  CALL PGSLW(7)
                  CALL PGSCI(ICDBCK)
               ELSE IF (LITYPE(III).EQ.-2) THEN          
                  CALL PGSLW(4)
                  CALL PGSCI(ICCBCK)
               ELSE IF (LITYPE(III).EQ.-3) THEN
                  CALL PGSLW(3)
                  CALL PGSCI(ICTBCK)
               ELSE
                  CALL PGSLW(2)
                  CALL PGSCI(ICLBCK)
               END IF
            END IF
            CALL PGLINE(2,X,Y)
 3000    CONTINUE
C     
         DO 45 I=NTOTM6,NTOT
            IF (ITYPE(I).EQ.1) THEN 
               CALL PGSCI(ICDFOR)
               CALL PGSLW(19)
               ZZ = ZARRAY(I)
               CALL PGPT(1,XARRAY(I)+0.2*ZZ,YARRAY(I)+0.3*ZZ,9)
            END IF 
 45      CONTINUE          
C     
         DO 4000 III=1,NTOT
            IF (ITYPE(III).EQ.1) THEN 
               CALL POLMAT(RPOL,THAXI1,PHAXI1,ALFA1)
            ELSE IF (ITYPE(III).EQ.2) THEN 
               CALL POLMAT(RPOL,THAXI2,PHAXI2,ALFA2)
            ELSE IF (ITYPE(III).EQ.3) THEN 
               CALL POLMAT(RPOL,THAXI3,PHAXI3,ALFA3)
            ELSE
               CALL POLMAT(RPOL,THAXI4,PHAXI4,ALFA4)
            END IF
            C(1)=XARRAY(III)
            C(2)=YARRAY(III)
            C(3)=ZARRAY(III)
            CALL MATMUL (C,RPOL,CROT)
            XARRAY(III)=CROT(1)+XOFF
            YARRAY(III)=CROT(2)+YOFF
            ZARRAY(III)=CROT(3)+ZOFF
 4000    CONTINUE
C     
         ALFA1 = ALFA1+1.5E-5*(1.0+2.0*L/4000.)
         ALFA2 = ALFA2-2.0E-5*(1.0+4.0*L/4000.)
         ALFA3 = ALFA3-4.0E-5*(1.0+3.0*L/4000.)
C     
         CALL PGEBUF
C     
 30   CONTINUE
C     
C-----------------------------------------------------------------------
      END

      SUBROUTINE MATMUL (VECTOR,RMATRX,ROTVEC)
C     
C     Matrix multiplication 
C     
      REAL VECTOR(3)
      REAL ROTVEC(3)
      REAL RMATRX(3,3)
C     
      ROTVEC(1)=RMATRX(1,1)*VECTOR(1)+RMATRX(1,2)*VECTOR(2)+
     #          RMATRX(1,3)*VECTOR(3)
      ROTVEC(2)=RMATRX(2,1)*VECTOR(1)+RMATRX(2,2)*VECTOR(2)+
     #          RMATRX(2,3)*VECTOR(3)
      ROTVEC(3)=RMATRX(3,1)*VECTOR(1)+RMATRX(3,2)*VECTOR(2)+
     #          RMATRX(3,3)*VECTOR(3)
C     
      RETURN
      END        

      SUBROUTINE POLMAT(RPOL,THAXI,PHAXI,ALFA)
C     
      REAL THAXI,PHAXI,ALFA
      REAL RPOL(3,3)
      REAL SINT,SINTQ,SINP,SINPQ,SINA
      REAL COST,COSTQ,COSP,COSPQ,COSA,EMCOSA
C     
      SINT = SIN(THAXI)
      COST = COS(THAXI)
      SINP = SIN(PHAXI)
      COSP = COS(PHAXI)
      SINA = SIN(ALFA)
      COSA = COS(ALFA)
      EMCOSA = 1.0-COSA
C     
      SINTQ = SINT*SINT
      COSTQ = COST*COST
      SINPQ = SINP*SINP
      COSPQ = COSP*COSP
C     
      RPOL(1,1) =  COSA+COSPQ*SINTQ*EMCOSA
      RPOL(2,1) =  COST*SINA+SINP*COSP*SINTQ*EMCOSA
      RPOL(3,1) = -SINP*SINT*SINA+SINT*COST*COSP*EMCOSA
      RPOL(1,2) = -COST*SINA+SINP*COSP*SINTQ*EMCOSA
      RPOL(2,2) =  COSA+SINPQ*SINTQ*EMCOSA
      RPOL(2,3) = -COSP*SINT*SINA+SINP*SINT*COST*EMCOSA
      RPOL(1,3) =  SINP*SINT*SINA+SINT*COST*COSP*EMCOSA
      RPOL(3,2) =  COSP*SINT*SINA+COST*SINT*SINP*EMCOSA
      RPOL(3,3) =  COSA+COSTQ*EMCOSA
C     
      RETURN 
      END

      SUBROUTINE SORTPP(N,ITYPE,RA1,RA2,RA3)
C     
      REAL RA1, RA2, RA3, RRA1, RRA2, RRA3
      INTEGER ITYPE(*), L, N, IR, I, J, IRRA1
      DIMENSION RA1(*), RA2(*), RA3(*)
      L=N/2+1
      IR=N
 10   CONTINUE
      IF(L.GT.1)THEN
         L=L-1
         RRA1=RA1(L)
         IRRA1=ITYPE(L)
         RRA2=RA2(L)
         RRA3=RA3(L)
      ELSE
         RRA1=RA1(IR)
         IRRA1=ITYPE(IR)
         RRA2=RA2(IR)
         RRA3=RA3(IR)
         RA1(IR)=RA1(1)
         ITYPE(IR)=ITYPE(1)
         RA2(IR)=RA2(1)
         RA3(IR)=RA3(1)
         IR=IR-1
         IF(IR.EQ.1)THEN
            RA1(1)=RRA1
            ITYPE(1)=IRRA1
            RA2(1)=RRA2
            RA3(1)=RRA3
            RETURN
         ENDIF
      ENDIF
      I=L
      J=L+L
 20   IF(J.LE.IR)THEN
         IF(J.LT.IR)THEN
            IF(RA1(J).LT.RA1(J+1))J=J+1
         ENDIF
         IF(RRA1.LT.RA1(J))THEN
            RA1(I)=RA1(J)
            ITYPE(I)=ITYPE(J)
            RA2(I)=RA2(J)
            RA3(I)=RA3(J)
            I=J
            J=J+J
         ELSE
            J=IR+1
         ENDIF
         GO TO 20
      ENDIF
      RA1(I)=RRA1
      ITYPE(I)=IRRA1
      RA2(I)=RRA2
      RA3(I)=RRA3
      GO TO 10
      END
C     
      SUBROUTINE SORTLI(N,RA1,IA1,IA2,IA3)
C     
      REAL RA1, RRA1
      INTEGER L, N, IR, I, J, IRA1, IRA2, IRA3, IA1, IA2, IA3
      DIMENSION RA1(*), IA1(*), IA2(*) , IA3(*)
      L=N/2+1
      IR=N
 10   CONTINUE
      IF(L.GT.1)THEN
         L=L-1
         RRA1=RA1(L)
         IRA1=IA1(L)
         IRA2=IA2(L)
         IRA3=IA3(L)
      ELSE
         RRA1=RA1(IR)
         IRA1=IA1(IR)
         IRA2=IA2(IR)
         IRA3=IA3(IR)
         RA1(IR)=RA1(1)
         IA1(IR)=IA1(1)
         IA2(IR)=IA2(1)
         IA3(IR)=IA3(1)
         IR=IR-1
         IF(IR.EQ.1)THEN
            RA1(1)=RRA1
            IA1(1)=IRA1
            IA2(1)=IRA2
            IA3(1)=IRA3
            RETURN
         ENDIF
      ENDIF
      I=L
      J=L+L
 20   IF(J.LE.IR)THEN
         IF(J.LT.IR)THEN
            IF(RA1(J).LT.RA1(J+1))J=J+1
         ENDIF
         IF(RRA1.LT.RA1(J))THEN
            RA1(I)=RA1(J)
            IA1(I)=IA1(J)
            IA2(I)=IA2(J)
            IA3(I)=IA3(J)
            I=J
            J=J+J
         ELSE
            J=IR+1
         ENDIF
         GO TO 20
      ENDIF
      RA1(I)=RRA1
      IA1(I)=IRA1
      IA2(I)=IRA2
      IA3(I)=IRA3
      GO TO 10
      END

      SUBROUTINE OFFSET (XOFF,YOFF,ZOFF)
C
      REAL XOFF,YOFF,ZOFF
C
      WRITE(*,*)' Rotation with shifting'
      XOFF=-0.0002
      YOFF=+0.0004
      ZOFF=-0.0002
      RETURN
      END

      SUBROUTINE CHNAX
     #     (THAXI3,PHAXI3,THAXI2,PHAXI2,THAXI1,PHAXI1)
C     
      REAL THAXI1,PHAXI1,PHAXI2,THAXI2,PHAXI3,THAXI3,PI
      PARAMETER (PI=3.14159265359)
C     
      THAXI3 = THAXI3 - PI*0.32
      PHAXI3 = PHAXI3 + PI*0.28
      THAXI2 = THAXI2 + PI*0.18
      PHAXI2 = PHAXI2 - PI*0.14
      THAXI1 = THAXI1 - PI*0.12
      PHAXI1 = PHAXI1 + PI*0.08
C     
      RETURN
      END






