* 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 * * Modified by Tsuguhiro TAMARIBUCHI for using the GrWin Library * Date: Sat, 21 Dec 2002 12:44:33 JST 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 GWopenx(IR,0,0,0,19,0,-1,' ') CALL GWsleep2(IR,0) CALL POLY3D * CALL PGCLOS CALL GWctime(lpt) WRITE(*,*) ' Elapsed time =',lpt,' (msec)' CALL GWquit(IR) 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 * INTEGER IR, MS 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) INTEGER MFOR(4), MBCK(4) 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) CALL GWvport(IR,0.0,0.0,1.0,1.0) CALL GWsize(IR, 5, LLW, LLH) IF(IR .NE. 0) THEN IF(LLW .GT. LLH) THEN LLW = LLH ELSE LLH = LLW ENDIF CALL GWsize(IR, -5, LLW, LLH) CALL GWsize(IR, -3, 0, 0) ENDIF CALL GWsize(IR,5,IX,IY) CALL GWindow(IR,-3.2,-3.2,3.2,3.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 ICDFOR = krgb( 0, 255, 0) CALL GWbegincmb(IR, 0, 1.0, 1.0) CALL GWsetpen(IR,-3,-1,-1,-1) CALL GWsetbrs(IR,ICDFOR,1,-1) CALL GWellipse(IR,0.1,0.1,0.9,0.9) CALL GWendcmb(MFOR(1)) * ICDBCK = 10 ICDBCK = krgb( 0, 255, 127) CALL GWbegincmb(IR, 0, 1.0, 1.0) CALL GWsetpen(IR,-3,-1,-1,-1) CALL GWsetbrs(IR,ICDBCK,1,-1) CALL GWellipse(IR,0.1,0.1,0.9,0.9) CALL GWendcmb(MBCK(1)) C * ICCFOR = 8 ICCFOR = krgb(255, 127, 0) CALL GWbegincmb(IR, 0, 1.0, 1.0) CALL GWsetpen(IR,-3,-1,-1,-1) CALL GWsetbrs(IR,ICCFOR,1,-1) CALL GWellipse(IR,0.1,0.1,0.9,0.9) CALL GWendcmb(MFOR(2)) * ICCBCK = 2 ICCBCK = krgb(255, 0, 0) CALL GWbegincmb(IR, 0, 1.0, 1.0) CALL GWsetpen(IR,-3,-1,-1,-1) CALL GWsetbrs(IR,ICCBCK,1,-1) CALL GWellipse(IR,0.1,0.1,0.9,0.9) CALL GWendcmb(MBCK(2)) C * ICTFOR = 5 ICTFOR = krgb( 0, 255, 255) CALL GWbegincmb(IR, 0, 1.0, 1.0) CALL GWsetpen(IR,-3,-1,-1,-1) CALL GWsetbrs(IR,ICTFOR,1,-1) CALL GWellipse(IR,0.1,0.1,0.9,0.9) CALL GWendcmb(MFOR(3)) * ICTBCK = 4 ICTBCK = krgb( 0, 0, 255) CALL GWbegincmb(IR, 0, 1.0, 1.0) CALL GWsetpen(IR,-3,-1,-1,-1) CALL GWsetbrs(IR,ICTBCK,1,-1) CALL GWellipse(IR,0.1,0.1,0.9,0.9) CALL GWendcmb(MBCK(3)) C * ICLFOR = 1 ICLFOR = 1 CALL GWbegincmb(IR, 0, 1.0, 1.0) CALL GWsetpen(IR,-3,-1,-1,-1) CALL GWsetbrs(IR,ICLFOR,1,-1) CALL GWellipse(IR,0.1,0.1,0.9,0.9) CALL GWendcmb(MFOR(4)) * ICLBCK = 7 ICLBCK = krgb(255, 255, 0) CALL GWbegincmb(IR, 0, 1.0, 1.0) CALL GWsetpen(IR,-3,-1,-1,-1) CALL GWsetbrs(IR,ICLBCK,1,-1) CALL GWellipse(IR,0.1,0.1,0.9,0.9) CALL GWendcmb(MBCK(4)) * CALL GWanchor(IR,1) 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 CALL GWsetogn(IR, 0, -1) 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) SZ = 0.36 MS = MBCK(1) ELSE IF (ITYPE(I).EQ.2) THEN * CALL PGSCI(ICCBCK) * CALL PGSLW(17) SZ = 0.34 MS = MBCK(2) ELSE IF (ITYPE(I).EQ.3) THEN * CALL PGSCI(ICTBCK) * CALL PGSLW(15) SZ = 0.30 MS = MBCK(3) ELSE * CALL PGSCI(ICLBCK) * CALL PGSLW(14) SZ = 0.28 MS = MBCK(4) END IF ZZ = ZARRAY(I) * CALL PGPT(1,XARRAY(I)+0.2*ZZ,YARRAY(I)+0.3*ZZ,9) CALL GWputcmb(IR,MS,XARRAY(I)+0.2*ZZ,YARRAY(I)+0.3*ZZ + ,SZ,SZ,0) 33 CONTINUE C DO 44 I=IPOS,NTOT IF (ITYPE(I).EQ.1) THEN * CALL PGSCI(ICDFOR) * CALL PGSLW(18) SZ = 0.36 MS = MFOR(1) ELSE IF (ITYPE(I).EQ.2) THEN * CALL PGSCI(ICCFOR) * CALL PGSLW(17) SZ = 0.34 MS = MFOR(2) ELSE IF (ITYPE(I).EQ.3) THEN * CALL PGSCI(ICTFOR) * CALL PGSLW(15) SZ = 0.30 MS = MFOR(3) ELSE * CALL PGSCI(ICLFOR) * CALL PGSLW(14) SZ = 0.28 MS = MFOR(4) END IF ZZ = ZARRAY(I) * CALL PGPT(1,XARRAY(I)+0.2*ZZ,YARRAY(I)+0.3*ZZ,9) CALL GWputcmb(IR,MS,XARRAY(I)+0.2*ZZ,YARRAY(I)+0.3*ZZ + ,SZ,SZ,0) 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) CALL GWsetpen(IR,ICDFOR,-1,10*2,-1) ELSE IF (LITYPE(III).EQ.2) THEN * CALL PGSLW(8) * CALL PGSCI(ICCFOR) CALL GWsetpen(IR,ICCFOR,-1,8*2,-1) ELSE IF (LITYPE(III).EQ.3) THEN * CALL PGSLW(6) * CALL PGSCI(ICTFOR) CALL GWsetpen(IR,ICTFOR,-1,6*2,-1) ELSE * CALL PGSLW(4) * CALL PGSCI(ICLFOR) CALL GWsetpen(IR,ICLFOR,-1,4*2,-1) END IF ELSE IF(LITYPE(III).EQ.-1) THEN * CALL PGSLW(7) * CALL PGSCI(ICDBCK) CALL GWsetpen(IR,ICDBCK,-1,7*2,-1) ELSE IF (LITYPE(III).EQ.-2) THEN * CALL PGSLW(4) * CALL PGSCI(ICCBCK) CALL GWsetpen(IR,ICCBCK,-1,4*2,-1) ELSE IF (LITYPE(III).EQ.-3) THEN * CALL PGSLW(3) * CALL PGSCI(ICTBCK) CALL GWsetpen(IR,ICTBCK,-1,3*2,-1) ELSE * CALL PGSLW(2) * CALL PGSCI(ICLBCK) CALL GWsetpen(IR,ICLBCK,-1,2*2,-1) END IF END IF * CALL PGLINE(2,X,Y) CALL GWline(IR,X(1),Y(1),X(2),Y(2)) 3000 CONTINUE C DO 45 I=NTOTM6,NTOT IF (ITYPE(I).EQ.1) THEN * CALL PGSCI(ICDFOR) * CALL PGSLW(19) SZ = 0.38 MS = MFOR(1) ZZ = ZARRAY(I) * CALL PGPT(1,XARRAY(I)+0.2*ZZ,YARRAY(I)+0.3*ZZ,9) CALL GWputcmb(IR,MS,XARRAY(I)+0.2*ZZ,YARRAY(I)+0.3*ZZ + ,SZ,SZ,0) 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 MMULT (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 CALL GWsetogn(IR, -1, 1) C 30 CONTINUE C C----------------------------------------------------------------------- END SUBROUTINE MMULT (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