home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 3 / TheARMClub_PDCD3.iso / hensa / maths / pgplot_1 / Examples / f77 / PGDEMO17 < prev    next >
Text File  |  1997-06-10  |  14KB  |  499 lines
  1. *     This program demonstrates animation and 3D geometry in PGPLOT. It 
  2. *     requires a fast, interactive display, e.g., /XWIN. Do not
  3. *     specify a hardcopy device. The speed of the animation is limited by
  4. *     the cpu speed of the host computer.
  5. *     
  6. *     Thanks to Dr Martin Weisser:
  7. *     Date: Sun, 18 May 1997 16:14:01 CET
  8. *     From: weisser@chclu.chemie.uni-konstanz.de
  9.  
  10.  
  11.       PROGRAM PGDEM17
  12. C-----------------------------------------------------------------------
  13. C     Demonstration program for PGPLOT.
  14. C-----------------------------------------------------------------------
  15.       INTEGER PGOPEN
  16. C
  17.       WRITE (*,*) 'PGPLOT: Demonstration of animation and 3D geometry'
  18.       WRITE (*,*) 'Select a fast, interactive device, e.g., /XWINDOW'
  19.       IF (PGOPEN('?') .LE. 0) STOP
  20.       CALL POLY3D
  21.       CALL PGCLOS
  22. C-----------------------------------------------------------------------
  23.       END
  24.  
  25.       SUBROUTINE POLY3D
  26. C     
  27.       INTEGER NFRAMS
  28. C     
  29.       INTEGER NTOT, NLIN, IPOS, IFIRST          
  30.       REAL T, T1, T2, T3, PI, W, W1, TET, TET1, ROT, ROT1
  31. C     
  32.       PARAMETER (NTOT=34)
  33.       PARAMETER (T=1.618)
  34.       PARAMETER (T1=1.0+T)
  35.       PARAMETER (T2=-T)
  36.       PARAMETER (T3=-T1)
  37.       PARAMETER (W=0.60*T)
  38.       PARAMETER (W1=-W)
  39.       PARAMETER (TET=0.37)
  40.       PARAMETER (TET1=-TET)
  41.       PARAMETER (ROT=0.13)
  42.       PARAMETER (ROT1=-ROT)
  43.       PARAMETER (NLIN=49)
  44. C     
  45.       INTEGER I, J, L, III, ILINE, NTOTM6
  46.       INTEGER ICDFOR, ICCFOR, ICTFOR, ICLFOR
  47.       INTEGER ICDBCK, ICCBCK, ICTBCK, ICLBCK
  48.       INTEGER ITYPE(NTOT), IARRAY(NLIN), JARRAY(NLIN), LITYPE(NLIN)
  49.       REAL RQ, ZZ
  50.       REAL THAXI1, PHAXI1, ALFA1, THAXI2, PHAXI2, ALFA2
  51.       REAL THAXI3, PHAXI3, ALFA3, THAXI4, PHAXI4, ALFA4 
  52.       REAL XOFF, YOFF, ZOFF
  53.       REAL XARRAY(NTOT), YARRAY(NTOT), ZARRAY(NTOT), DISTAN(NLIN)
  54.       REAL POLYS(3,NTOT), X(2), Y(2), C(3), CROT(3), RPOL(3,3)
  55.       PARAMETER (PI=3.14159265359)
  56. C     
  57. C     Cartesian coordinates of the polygons 
  58. C     
  59.       DATA POLYS/ T, T, T,       T, T,T2,
  60.      D     T,T2, T,       T,T2,T2,
  61.      D     T2, T, T,      T2, T,T2,
  62.      D     T2,T2, T,      T2,T2,T2,
  63.      D     T1,1.0,0.0,    T1,-1.0,0.0,
  64.      D     T3,1.0,0.0,    T3,-1.0,0.0,
  65.      D     0.0,T1,1.0,    0.0,T1,-1.0,
  66.      D     0.0,T3,1.0,    0.0,T3,-1.0,
  67.      D     1.0,0.0,T1,    -1.0,0.0,T1,
  68.      D     1.0,0.0,T3,   -1.0,0.0,T3, 
  69.      C     W,    W,    W,     W,    W,   W1,
  70.      C     W,   W1,    W,     W,   W1,   W1,
  71.      C     W1,    W,    W,    W1,   W1,    W,
  72.      C     W1,    W,   W1,    W1,   W1,   W1,                
  73.      T     TET,  TET,  TET, TET1, TET1,  TET,
  74.      T     TET1,  TET, TET1,  TET, TET1, TET1, 
  75.      L     ROT,  0.0,  0.0, ROT1,  0.0,  0.0/
  76. C     
  77.       DATA ITYPE/1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
  78.      C     2,2,2,2,2,2,2,2,
  79.      T     3,3,3,3,
  80.      L     4,4/
  81. C     
  82. C     Initialize the plot (no labels).
  83. C     
  84.       CALL PGENV(-3.2,3.2,-3.2,3.2,1,-2)
  85. C     
  86. C     Switch from page to page without typing return.  
  87. C     
  88.       CALL PGASK(.FALSE.)
  89. C     
  90. C     Rotation axis of the polygons 
  91. C     
  92.       THAXI1 = PI/4.0
  93.       PHAXI1 = PI/4.5
  94.       ALFA1  = 0.0
  95.       THAXI2 = PI/6.0
  96.       PHAXI2 = 0.0
  97.       ALFA2  = 0.02
  98.       THAXI3 = PI/2.0
  99.       PHAXI3 = -PI/3.0
  100.       ALFA3  = 0.0
  101.       THAXI4 = -0.03
  102.       PHAXI4 = PI/7.0
  103.       ALFA4  = 0.9
  104. C     
  105.       XOFF=0.0
  106.       YOFF=0.0
  107.       ZOFF=0.0
  108.       NTOTM6=NTOT-6
  109. C     
  110. C     Colors  
  111. C     
  112.       ICDFOR = 3
  113.       ICDBCK = 10
  114. C     
  115.       ICCFOR = 8
  116.       ICCBCK = 2
  117. C     
  118.       ICTFOR = 5
  119.       ICTBCK = 4
  120. C     
  121.       ICLFOR = 1
  122.       ICLBCK = 7
  123. C     
  124.       IPOS=1
  125. C     
  126.       WRITE(*,*)' Rotation with increasing velocity'
  127. C     
  128.       NFRAMS=3500 
  129. C     
  130.       DO 12 I=1,NTOT
  131.          XARRAY(I) = POLYS(1,I)
  132.          YARRAY(I) = POLYS(2,I)
  133.          ZARRAY(I) = POLYS(3,I)
  134.  12   CONTINUE
  135. C     
  136.       DO 30 L=1,NFRAMS
  137. C     
  138.          CALL PGBBUF
  139.          CALL PGERAS
  140. C     
  141.          CALL SORTPP(NTOT,ITYPE,ZARRAY,YARRAY,XARRAY)
  142. C     
  143.          IFIRST=0
  144.          DO 13 I=1,NTOT
  145.             IF((ZARRAY(I).GE.0.0).AND.(IFIRST.EQ.0)) THEN
  146.                IFIRST = 1
  147.                IPOS = I
  148.             END IF
  149.  13      CONTINUE
  150. C     
  151.          IF(L.EQ.2800) CALL OFFSET (XOFF,YOFF,ZOFF)
  152.          IF (MOD(L,500).EQ.0) THEN 
  153.             CALL CHNAX(THAXI3,PHAXI3,THAXI2,PHAXI2,THAXI1,PHAXI1)
  154.          END IF
  155. C     
  156.          DO 33 I=1,IPOS-1
  157.             IF (ITYPE(I).EQ.1) THEN 
  158.                CALL PGSCI(ICDBCK)
  159.                CALL PGSLW(18)
  160.             ELSE IF (ITYPE(I).EQ.2) THEN 
  161.                CALL PGSCI(ICCBCK)
  162.                CALL PGSLW(17)
  163.             ELSE IF (ITYPE(I).EQ.3) THEN
  164.                CALL PGSCI(ICTBCK)
  165.                CALL PGSLW(15)
  166.             ELSE 
  167.                CALL PGSCI(ICLBCK)
  168.                CALL PGSLW(14)
  169.             END IF
  170.             ZZ = ZARRAY(I)
  171.             CALL PGPT(1,XARRAY(I)+0.2*ZZ,YARRAY(I)+0.3*ZZ,9)
  172.  33      CONTINUE
  173. C     
  174.          DO 44 I=IPOS,NTOT
  175.             IF (ITYPE(I).EQ.1) THEN 
  176.                CALL PGSCI(ICDFOR)
  177.                CALL PGSLW(18)
  178.             ELSE IF (ITYPE(I).EQ.2) THEN 
  179.                CALL PGSCI(ICCFOR)
  180.                CALL PGSLW(17)
  181.             ELSE IF (ITYPE(I).EQ.3) THEN
  182.                CALL PGSCI(ICTFOR)
  183.                CALL PGSLW(15)
  184.             ELSE 
  185.                CALL PGSCI(ICLFOR)
  186.                CALL PGSLW(14)
  187.             END IF
  188.             ZZ = ZARRAY(I)
  189.             CALL PGPT(1,XARRAY(I)+0.2*ZZ,YARRAY(I)+0.3*ZZ,9)
  190.  44      CONTINUE
  191. C     
  192.          ILINE=0
  193. C     
  194.          DO 2000 I=2,NTOT
  195.             DO 1000 J=1,I-1
  196.                IF (ITYPE(I).EQ.ITYPE(J)) THEN
  197.                   RQ = 0.0
  198.                   RQ = RQ + ( (XARRAY(I)-XARRAY(J))**2+
  199.      #                 (YARRAY(I)-YARRAY(J))**2+
  200.      #                 (ZARRAY(I)-ZARRAY(J))**2  )
  201. C     
  202.                   IF ( ((RQ-0.0676)  .LT.0.001).OR.
  203.      #                 ((RQ-1.095199).LT.0.001).OR.          
  204.      #                 ((RQ-3.769809).LT.0.001).OR.
  205.      #                 ((RQ-4.000000).LT.0.001)     ) THEN 
  206.                      ILINE = ILINE + 1
  207.                      DISTAN(ILINE) = ZARRAY(I)+ZARRAY(J)
  208.                      IF(DISTAN(ILINE).LT.0.0) THEN 
  209.                         LITYPE(ILINE) = -ITYPE(I)
  210.                      ELSE  
  211.                         LITYPE(ILINE) =  ITYPE(I)
  212.                      END IF 
  213.                      IARRAY(ILINE) = I
  214.                      JARRAY(ILINE) = J
  215.                   END IF
  216.                END IF
  217.  1000       CONTINUE
  218.  2000    CONTINUE
  219. C     
  220.          CALL SORTLI(ILINE,DISTAN,IARRAY,JARRAY,LITYPE)
  221. C     
  222.          DO 3000 III=1,ILINE
  223.             I=IARRAY(III)
  224.             J=JARRAY(III)               
  225.             ZZ = ZARRAY(I)
  226.             X(1) = XARRAY(I)+0.2*ZZ
  227.             Y(1) = YARRAY(I)+0.3*ZZ
  228.             ZZ = ZARRAY(J)
  229.             X(2) = XARRAY(J)+0.2*ZZ
  230.             Y(2) = YARRAY(J)+0.3*ZZ
  231.             IF (LITYPE(III).GT.0) THEN
  232.                IF(LITYPE(III).EQ.1) THEN  
  233.                   CALL PGSLW(10)
  234.                   CALL PGSCI(ICDFOR)
  235.                ELSE IF (LITYPE(III).EQ.2) THEN          
  236.                   CALL PGSLW(8)
  237.                   CALL PGSCI(ICCFOR)
  238.                ELSE IF (LITYPE(III).EQ.3) THEN
  239.                   CALL PGSLW(6)
  240.                   CALL PGSCI(ICTFOR)
  241.                ELSE
  242.                   CALL PGSLW(4)
  243.                   CALL PGSCI(ICLFOR)
  244.                END IF
  245.             ELSE
  246.                IF(LITYPE(III).EQ.-1) THEN  
  247.                   CALL PGSLW(7)
  248.                   CALL PGSCI(ICDBCK)
  249.                ELSE IF (LITYPE(III).EQ.-2) THEN          
  250.                   CALL PGSLW(4)
  251.                   CALL PGSCI(ICCBCK)
  252.                ELSE IF (LITYPE(III).EQ.-3) THEN
  253.                   CALL PGSLW(3)
  254.                   CALL PGSCI(ICTBCK)
  255.                ELSE
  256.                   CALL PGSLW(2)
  257.                   CALL PGSCI(ICLBCK)
  258.                END IF
  259.             END IF
  260.             CALL PGLINE(2,X,Y)
  261.  3000    CONTINUE
  262. C     
  263.          DO 45 I=NTOTM6,NTOT
  264.             IF (ITYPE(I).EQ.1) THEN 
  265.                CALL PGSCI(ICDFOR)
  266.                CALL PGSLW(19)
  267.                ZZ = ZARRAY(I)
  268.                CALL PGPT(1,XARRAY(I)+0.2*ZZ,YARRAY(I)+0.3*ZZ,9)
  269.             END IF 
  270.  45      CONTINUE          
  271. C     
  272.          DO 4000 III=1,NTOT
  273.             IF (ITYPE(III).EQ.1) THEN 
  274.                CALL POLMAT(RPOL,THAXI1,PHAXI1,ALFA1)
  275.             ELSE IF (ITYPE(III).EQ.2) THEN 
  276.                CALL POLMAT(RPOL,THAXI2,PHAXI2,ALFA2)
  277.             ELSE IF (ITYPE(III).EQ.3) THEN 
  278.                CALL POLMAT(RPOL,THAXI3,PHAXI3,ALFA3)
  279.             ELSE
  280.                CALL POLMAT(RPOL,THAXI4,PHAXI4,ALFA4)
  281.             END IF
  282.             C(1)=XARRAY(III)
  283.             C(2)=YARRAY(III)
  284.             C(3)=ZARRAY(III)
  285.             CALL MATMUL (C,RPOL,CROT)
  286.             XARRAY(III)=CROT(1)+XOFF
  287.             YARRAY(III)=CROT(2)+YOFF
  288.             ZARRAY(III)=CROT(3)+ZOFF
  289.  4000    CONTINUE
  290. C     
  291.          ALFA1 = ALFA1+1.5E-5*(1.0+2.0*L/4000.)
  292.          ALFA2 = ALFA2-2.0E-5*(1.0+4.0*L/4000.)
  293.          ALFA3 = ALFA3-4.0E-5*(1.0+3.0*L/4000.)
  294. C     
  295.          CALL PGEBUF
  296. C     
  297.  30   CONTINUE
  298. C     
  299. C-----------------------------------------------------------------------
  300.       END
  301.  
  302.       SUBROUTINE MATMUL (VECTOR,RMATRX,ROTVEC)
  303. C     
  304. C     Matrix multiplication 
  305. C     
  306.       REAL VECTOR(3)
  307.       REAL ROTVEC(3)
  308.       REAL RMATRX(3,3)
  309. C     
  310.       ROTVEC(1)=RMATRX(1,1)*VECTOR(1)+RMATRX(1,2)*VECTOR(2)+
  311.      #          RMATRX(1,3)*VECTOR(3)
  312.       ROTVEC(2)=RMATRX(2,1)*VECTOR(1)+RMATRX(2,2)*VECTOR(2)+
  313.      #          RMATRX(2,3)*VECTOR(3)
  314.       ROTVEC(3)=RMATRX(3,1)*VECTOR(1)+RMATRX(3,2)*VECTOR(2)+
  315.      #          RMATRX(3,3)*VECTOR(3)
  316. C     
  317.       RETURN
  318.       END        
  319.  
  320.       SUBROUTINE POLMAT(RPOL,THAXI,PHAXI,ALFA)
  321. C     
  322.       REAL THAXI,PHAXI,ALFA
  323.       REAL RPOL(3,3)
  324.       REAL SINT,SINTQ,SINP,SINPQ,SINA
  325.       REAL COST,COSTQ,COSP,COSPQ,COSA,EMCOSA
  326. C     
  327.       SINT = SIN(THAXI)
  328.       COST = COS(THAXI)
  329.       SINP = SIN(PHAXI)
  330.       COSP = COS(PHAXI)
  331.       SINA = SIN(ALFA)
  332.       COSA = COS(ALFA)
  333.       EMCOSA = 1.0-COSA
  334. C     
  335.       SINTQ = SINT*SINT
  336.       COSTQ = COST*COST
  337.       SINPQ = SINP*SINP
  338.       COSPQ = COSP*COSP
  339. C     
  340.       RPOL(1,1) =  COSA+COSPQ*SINTQ*EMCOSA
  341.       RPOL(2,1) =  COST*SINA+SINP*COSP*SINTQ*EMCOSA
  342.       RPOL(3,1) = -SINP*SINT*SINA+SINT*COST*COSP*EMCOSA
  343.       RPOL(1,2) = -COST*SINA+SINP*COSP*SINTQ*EMCOSA
  344.       RPOL(2,2) =  COSA+SINPQ*SINTQ*EMCOSA
  345.       RPOL(2,3) = -COSP*SINT*SINA+SINP*SINT*COST*EMCOSA
  346.       RPOL(1,3) =  SINP*SINT*SINA+SINT*COST*COSP*EMCOSA
  347.       RPOL(3,2) =  COSP*SINT*SINA+COST*SINT*SINP*EMCOSA
  348.       RPOL(3,3) =  COSA+COSTQ*EMCOSA
  349. C     
  350.       RETURN 
  351.       END
  352.  
  353.       SUBROUTINE SORTPP(N,ITYPE,RA1,RA2,RA3)
  354. C     
  355.       REAL RA1, RA2, RA3, RRA1, RRA2, RRA3
  356.       INTEGER ITYPE(*), L, N, IR, I, J, IRRA1
  357.       DIMENSION RA1(*), RA2(*), RA3(*)
  358.       L=N/2+1
  359.       IR=N
  360.  10   CONTINUE
  361.       IF(L.GT.1)THEN
  362.          L=L-1
  363.          RRA1=RA1(L)
  364.          IRRA1=ITYPE(L)
  365.          RRA2=RA2(L)
  366.          RRA3=RA3(L)
  367.       ELSE
  368.          RRA1=RA1(IR)
  369.          IRRA1=ITYPE(IR)
  370.          RRA2=RA2(IR)
  371.          RRA3=RA3(IR)
  372.          RA1(IR)=RA1(1)
  373.          ITYPE(IR)=ITYPE(1)
  374.          RA2(IR)=RA2(1)
  375.          RA3(IR)=RA3(1)
  376.          IR=IR-1
  377.          IF(IR.EQ.1)THEN
  378.             RA1(1)=RRA1
  379.             ITYPE(1)=IRRA1
  380.             RA2(1)=RRA2
  381.             RA3(1)=RRA3
  382.             RETURN
  383.          ENDIF
  384.       ENDIF
  385.       I=L
  386.       J=L+L
  387.  20   IF(J.LE.IR)THEN
  388.          IF(J.LT.IR)THEN
  389.             IF(RA1(J).LT.RA1(J+1))J=J+1
  390.          ENDIF
  391.          IF(RRA1.LT.RA1(J))THEN
  392.             RA1(I)=RA1(J)
  393.             ITYPE(I)=ITYPE(J)
  394.             RA2(I)=RA2(J)
  395.             RA3(I)=RA3(J)
  396.             I=J
  397.             J=J+J
  398.          ELSE
  399.             J=IR+1
  400.          ENDIF
  401.          GO TO 20
  402.       ENDIF
  403.       RA1(I)=RRA1
  404.       ITYPE(I)=IRRA1
  405.       RA2(I)=RRA2
  406.       RA3(I)=RRA3
  407.       GO TO 10
  408.       END
  409. C     
  410.       SUBROUTINE SORTLI(N,RA1,IA1,IA2,IA3)
  411. C     
  412.       REAL RA1, RRA1
  413.       INTEGER L, N, IR, I, J, IRA1, IRA2, IRA3, IA1, IA2, IA3
  414.       DIMENSION RA1(*), IA1(*), IA2(*) , IA3(*)
  415.       L=N/2+1
  416.       IR=N
  417.  10   CONTINUE
  418.       IF(L.GT.1)THEN
  419.          L=L-1
  420.          RRA1=RA1(L)
  421.          IRA1=IA1(L)
  422.          IRA2=IA2(L)
  423.          IRA3=IA3(L)
  424.       ELSE
  425.          RRA1=RA1(IR)
  426.          IRA1=IA1(IR)
  427.          IRA2=IA2(IR)
  428.          IRA3=IA3(IR)
  429.          RA1(IR)=RA1(1)
  430.          IA1(IR)=IA1(1)
  431.          IA2(IR)=IA2(1)
  432.          IA3(IR)=IA3(1)
  433.          IR=IR-1
  434.          IF(IR.EQ.1)THEN
  435.             RA1(1)=RRA1
  436.             IA1(1)=IRA1
  437.             IA2(1)=IRA2
  438.             IA3(1)=IRA3
  439.             RETURN
  440.          ENDIF
  441.       ENDIF
  442.       I=L
  443.       J=L+L
  444.  20   IF(J.LE.IR)THEN
  445.          IF(J.LT.IR)THEN
  446.             IF(RA1(J).LT.RA1(J+1))J=J+1
  447.          ENDIF
  448.          IF(RRA1.LT.RA1(J))THEN
  449.             RA1(I)=RA1(J)
  450.             IA1(I)=IA1(J)
  451.             IA2(I)=IA2(J)
  452.             IA3(I)=IA3(J)
  453.             I=J
  454.             J=J+J
  455.          ELSE
  456.             J=IR+1
  457.          ENDIF
  458.          GO TO 20
  459.       ENDIF
  460.       RA1(I)=RRA1
  461.       IA1(I)=IRA1
  462.       IA2(I)=IRA2
  463.       IA3(I)=IRA3
  464.       GO TO 10
  465.       END
  466.  
  467.       SUBROUTINE OFFSET (XOFF,YOFF,ZOFF)
  468. C
  469.       REAL XOFF,YOFF,ZOFF
  470. C
  471.       WRITE(*,*)' Rotation with shifting'
  472.       XOFF=-0.0002
  473.       YOFF=+0.0004
  474.       ZOFF=-0.0002
  475.       RETURN
  476.       END
  477.  
  478.       SUBROUTINE CHNAX
  479.      #     (THAXI3,PHAXI3,THAXI2,PHAXI2,THAXI1,PHAXI1)
  480. C     
  481.       REAL THAXI1,PHAXI1,PHAXI2,THAXI2,PHAXI3,THAXI3,PI
  482.       PARAMETER (PI=3.14159265359)
  483. C     
  484.       THAXI3 = THAXI3 - PI*0.32
  485.       PHAXI3 = PHAXI3 + PI*0.28
  486.       THAXI2 = THAXI2 + PI*0.18
  487.       PHAXI2 = PHAXI2 - PI*0.14
  488.       THAXI1 = THAXI1 - PI*0.12
  489.       PHAXI1 = PHAXI1 + PI*0.08
  490. C     
  491.       RETURN
  492.       END
  493.  
  494.  
  495.  
  496.  
  497.  
  498.  
  499.