PROGRAM DW2DS C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * ****** DOUBLE PRECISION ****** * C * PROGRAM DW2DS: (THIRD VERSION 11-12-89; CONTINUOUS EIGEN-VECTORS) * C * CALCULATES DATA SETS FOR DEFINING THE FULL 3-D FIGURES OF * C * THE PHASE VELOCITY SURFACES BY SWEEPING THRU NPHI PLANES * C * WHICH ARE ORIENTED FROM THE AXIS IXI BY THE ANGLE PHI. * C * WITHIN EACH PLANE THREE CURVES ARE CALCULATED BY INCREMENT- * C * ALLY SOLVING CHRISTOFFEL'S EQUATIONS FOR THE EIGEN-VALUES * C * (WAVE PHASE VELOCITIES) AT NTHETA POINTS STARTING AT THE * C * IORTH AXIS. RESULTS ARE STORED IN THREE FILES CORRESPONDING * C * TO EACH OF THE THREE SURFACES WHERE THE PARTICLE DISPLACE- * C * MENT DIRECTION COSINES ARE CONTINUOUS, EXCEPT ALONG SPECIAL * C * LINES WITHIN THE PRINCIPAL MATERIAL PLANES (SEE SUBROUTINE * C * WENA W-AVE E-XCHANGE N-OT A-LLOWED). THE RCC COORDINATES ARE * C * CALCULATED FOR EACH OF THE THREE WAVES SURFACES AND AT EACH * C * POINT ALONG A SURFACE THE DISPLACEMENT DIRECTION COSINES * C * ARE CALCULATED AS A TILT PLANE ANGLE FROM THE ROTATED PLANE * C * AND A POLARIZATION WITHIN THE ROTATED PLANE. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C IMPLICIT REAL*8 (A-H,O-Z) COMMON /PAR/NPHI,NTHETA,NTOL,IREG REAL*8 C(6,6),N(3),V(3),P(3,3),B(3,3),OP(3,3),DELTA(3),ZETA(3) REAL*8 TILT(3),LT(3),MONE C C OPEN FILES FOR TAPE5=INPUT AND (TAPE6 THRU TAPE10)=OUTPUT C OPEN(5,FILE='dw2ds.data',STATUS='OLD',ERR=8888) OPEN(6,FILE='dw2ds.out',STATUS='UNKNOWN',ERR=8888) OPEN(7,FILE='dsurf1',STATUS='UNKNOWN',ERR=8888) OPEN(8,FILE='dsurf2',STATUS='UNKNOWN',ERR=8888) OPEN(9,FILE='dsurf3',STATUS='UNKNOWN',ERR=8888) OPEN(10,FILE='dw2ds.err',STATUS='UNKNOWN',ERR=8888) OPEN(11,FILE='g_vert1.dat',STATUS='UNKNOWN',ERR=8888) OPEN(12,FILE='g_vert2.dat',STATUS='UNKNOWN',ERR=8888) OPEN(13,FILE='g_vert3.dat',STATUS='UNKNOWN',ERR=8888) C C READ MATERIAL TYPE AND PROPERTIES FROM INPUT FILE C CALL IN(RHO,C) C C WRITE MATERIAL TYPE, DENSITY AND STIFFNESSES ON OUTPUT FILE C CALL OUT(RHO,C) C C DEFINE ROTATION PLANE (POSITIVE PHI FROM IXI AND POSITIVE THETA FROM IORTH) C IXI=1 IXJ=2 C C CALCULATE ORTHOGONAL AXIS OF ROTATION (IORH PERPENDICULAR TO IXI-IXJ PLANE) C C IORTH=DABS(IXI*IXJ-5.0) IORTH=3 C C WITHIN EACH PLANE LOOP THRU NTHETA ANGULAR INCREMENTS STARTING C AT IORTH WITH THETAR=0 RADIANS AND ENDING AT THETAR=PI RADIANS. C EACH PLANE DEFINED ABOVE IS ROTATED AWAY FROM THE IXI AXIS C STARTING WITH PHI=0 RADIANS AND ENDING AT PHI=2*PI RADIANS. C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * POLAR COORDINATES: * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C DEFINE PI RADIANS C C PI=DACOS(-1.) MONE=-1.0 PI=DACOS(MONE) C C INITIALIZE TOLERANCE FOR MAXIMUM ANGLE SEPERATING TWO CONSEQUTIVE C EIGEN-VECTORS: 0.10=84.3 DEG. OR GREATER, 0.20=78.5 DEG. OR GREATER C CTOL=0.0 C C SET TOLERANCE FOR OSCILLATING EIGENVECTORS NEAR PRINCIPAL AXES C NTOL=5 C C WRITE TOLERANCES TO OUTPUT FILE C WRITE(6,676)CTOL,NTOL 676 FORMAT(80(1H#),/,' CALC. TOLERANCES, CTOL=',F5.2,' NTOL=',I2) C C INITIALIZE PARAMETERS THAT DEFINE REGION FOR PLOTTING C (AN OCTANT 90X90 DEG., IREG=1) (FULL SPHERE 360X180 DEG., IREG=4) C IREG=4 NPHI=121 RPHI=360.0 NTHETA=61 RTHETA=180.0 C c IREG=4 c NPHI=361 c RPHI=360.0 c NTHETA=181 c RTHETA=180.0 C C IREG=4 C NPHI=721 C RPHI=360.0 C NTHETA=361 C RTHETA=180.0 C C IREG=1 C NPHI=301 C RPHI=90.0 C NTHETA=301 C RTHETA=90.0 C C IREG=4 C NPHI=901 C RPHI=360.0 C NTHETA=451 C RTHETA=180.0 C C CALCULATE ROTATIONAL INCREMENTS FOR PHI AND THETA C DPHIR=(RPHI/(NPHI-1))*PI/180.0 DTHETAR=(RTHETA/(NTHETA-1))*PI/180.0 C C INITIALIZE WAVE IDENTITY NUMBERS ALONG THE IORTH-AXIS C (THE SLOWEST WAVE IS LABELED IW1=1, ETC.) C IW1=1 IW2=2 IW3=3 C C LOOP THRU NPHI PLANES THAT ARE PERPENDICULAR TO THE IXI-IXJ PLANE. C WITHIN THE IXI-IXJ PLANE, PHI STARTS (PHI=0 DEG.) AND C FINISHES (PHI=360 DEG.) AT THE IXI AXIS. C PHIR=0.0 DO 990 IPHI=1,NPHI C C WITHIN EACH PLANE DEFINED BY PHI, LOOP THRU NTHETA POINTS C FROM THE IORTH AXIS AT THETA=0 DEG. TO THE -IORTH AXIS AT 180 DEGREES C THETAR=0.0 DO 980 ITHETA=1,NTHETA N(IXI)=DSIN(THETAR)*DCOS(PHIR) N(IXJ)=DSIN(THETAR)*DSIN(PHIR) N(IORTH)=DCOS(THETAR) C C CALL ROUTINE TO SOLVE FOR : C PHASE VELOCITIES, V(J) C DISPL. DIRECT. COSINES, P(1-3,J) C ENERGY FLUX DIRECT. COSINES, B(1-3,J) C NUMERICAL INDEX OF CONVERGENCE, CONV C CALL FLUX(RHO,C,N,V,P,CONV,B) C C TEST FOR DISCONTINUITIES IN EIGEN-VECTORS (PARTICLE DISPL.DIRECT. COSINES) C AND REARRANGE VELOCITY(EIGEN-VALUESS) TO FORM VELOCTIY SURFACES WITH C CONTINUOUS PARTICLE DISPLACEMENT DIRECITON COSINES C CALL DISC(CTOL,IW1,IW2,IW3,IPHI,ITHETA,OP,P,V,B) C C CALCULATE FLUX DEVIATIONS AND PARTICLE DEVIATIONS FROM WAVE NORMAL C CALL DEV(IPHI,ITHETA,PI,N,P,B,DELTA,ZETA) C C CALCULATE TILT ANGLE, TILT(J), AWAY FROM ROTATED PLANE AND POLARIZATION C ANGLE, LT(J), OF THE DISPLACEMENT WITHIN THE ROTATED PLANE C CALL TLT(IPHI,ITHETA,PI,N,P,PHIR,TILT,LT) C C CALCULATE (IXI,IXJ,IXK) COORD.S FOR THE THREE PHASE VELOCITY SURFACE C AND WRITE RESULTS TO OUTPUT FILES 7,8,AND 9 C XV1=V(1)*DSIN(THETAR)*DCOS(PHIR) YV1=V(1)*DSIN(THETAR)*DSIN(PHIR) ZV1=V(1)*DCOS(THETAR) C XV2=V(2)*DSIN(THETAR)*DCOS(PHIR) YV2=V(2)*DSIN(THETAR)*DSIN(PHIR) ZV2=V(2)*DCOS(THETAR) C XV3=V(3)*DSIN(THETAR)*DCOS(PHIR) YV3=V(3)*DSIN(THETAR)*DSIN(PHIR) ZV3=V(3)*DCOS(THETAR) C WRITE(11,971)XV1,YV1,ZV1,ZETA(1) WRITE(7,970)IPHI,ITHETA,XV1,YV1,ZV1,LT(1),TILT(1) * ,DELTA(1),ZETA(1) WRITE(12,971)XV2,YV2,ZV2,ZETA(2) WRITE(8,970)IPHI,ITHETA,XV2,YV2,ZV2,LT(2),TILT(2) * ,DELTA(2),ZETA(2) WRITE(13,971)XV3,YV3,ZV3,ZETA(3) WRITE(9,970)IPHI,ITHETA,XV3,YV3,ZV3,LT(3),TILT(3) * ,DELTA(3),ZETA(3) 971 FORMAT(3(1X,E12.5),1x,',',1x,f6.1) 970 FORMAT(2(1X,I3),3(1X,E12.5),2(1X,F6.1),2(1x,f6.1)) C C WRITE RESULTS TO OUTPUT FILE 6 C PHI=PHIR*180.0/PI THETA=THETAR*180.0/PI WRITE(6,31)PHI,THETA,CONV,(N(I),I=1,3) 31 FORMAT(' PHI=',F5.1,' THETA=',F5.1,/,6X,'WAVE NORMAL COSINES', * 11X,'NUMERICAL INDEX CONVERGED TO ',E10.3,/, * 3X,3(1X,F7.4)) WRITE(6,32) 32 FORMAT(' VELOCITY DISPLACEMENT COSINES ENERGY FLUX COSINES ' * ,' DELTA ZETA TILT LT') DO 35 J=1,3 35 WRITE(6,36)V(J),(P(I,J),I=1,3),(B(I,J),I=1,3),DELTA(J),ZETA(J) * ,TILT(J),LT(J) 36 FORMAT(1X,E10.3,' ',3(F6.3,1X),' ',3(F6.3,1X),' ',F5.1,' ', * F5.1,' ',F5.1,' ',F5.1) WRITE(6,37) 37 FORMAT(1X,94(1H*)) C C INCREMENT THETA IN ROTATED PLANE C THETAR=THETAR+DTHETAR C 980 CONTINUE C END LOOP ON THETA IN ROTATED PLANE C C INCREMENT PHI TO A NEW ROTATED PLANE C PHIR=PHIR+DPHIR C 990 CONTINUE C END LOOP ON ROTATED PLANES C C C WRITE OUT FINAL ROTATIONAL ANGLES C WRITE(6,1234)IPHI,ITHETA 1234 FORMAT(80(1H%),/,' FINAL PHI=',I3,' FINAL THETA=',I3) C C 8888 STOP END C C######################################################################## SUBROUTINE FLUX(RHO,C,N,V,P,CONV,B) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * SUBROUTINE FLUX: * C * 1. SETS UP EIGEN VALUE PROBLEM TO SOLVE FOR * C * EIGEN VALUES, RHO*V(J)**2 (SOLVES FOR VELOCITIES V(J))* C * EIGEN VECTORS (DISPL. DIRECT. COSINES), P(1-3,J) * C * 2. SOLVES FOR ENERGY FLUX DIRECT. COSINES, B(1-3,J) * C * 3. CALCULATES NUMERICALLY CONVERGED PERFORMANCE INDEX, CONV* C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * INPUT : RHO, MATERIAL DENSITY * C * C(1-6,1-6), MATERIAL STIFFNESSES * C * N(1-3), WAVE NORMAL DIRECTION COSINES * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * NOTE: MATERIAL SYMMETRY (MONOCLINIC) * C * * C * C(1,1) C(1,2) C(1,3) 0 0 C(1,6) * C * C(2,2) C(2,3) 0 0 C(2,6) * C * C(3,3) 0 0 C(3,6) * C * C(4,4) C(4,5) 0 * C * SYMMETRIC C(5,5) 0 * C * C(6,6) * C * * * * * * * * * * * * * * * * * * * * * * * * * C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 C(6,6),N(3),A(3,3),V(3),P(3,3),B(3,3) REAL*8 N1,N2,N3,L1,L2,L3 C C WAVE NORMAL DIRECT. COSINES SET EXACTLY TO ZERO OR ONE C ALONG VERTICAL OR HORIZONTAL AXIS LOCATIONS C DO 5 I=1,3 IF(N(I).GT.0.0)SIGN=+1.0 IF(N(I).LT.0.0)SIGN=-1.0 ABSN=DABS(N(I)) IF(ABSN.LT.0.00001)N(I)=+0.0 IF(ABSN.GT.0.99999)N(I)=SIGN 5 CONTINUE N1=N(1) N2=N(2) N3=N(3) C C CHRISTOFFEL TERMS ASSIGN TO SYMMETRIC MATRIX A(3,3) C A(1,1)=C(1,1)*N1*N1+C(6,6)*N2*N2+C(5,5)*N3*N3+C(1,6)*2.*N1*N2 A(2,2)=C(6,6)*N1*N1+C(2,2)*N2*N2+C(4,4)*N3*N3+C(2,6)*2.*N1*N2 A(3,3)=C(5,5)*N1*N1+C(4,4)*N2*N2+C(3,3)*N3*N3+C(4,5)*2.*N1*N2 A(2,3)=(C(2,3)+C(4,4))*N2*N3+(C(3,6)+C(4,5))*N3*N1 A(1,3)=(C(3,6)+C(4,5))*N2*N3+(C(1,3)+C(5,5))*N3*N1 A(1,2)=C(1,6)*N1*N1+C(2,6)*N2*N2+C(4,5)*N3*N3 * +(C(1,2)+C(6,6))*N1*N2 A(2,1)=A(1,2) A(3,1)=A(1,3) A(3,2)=A(2,3) C C EIGENVALUE PROBLEM DEFINED C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C ZERO DISPLACEMENT DIRECTION COSINES P(I,J) (EIGEN VECTORS) C AND PHASE VELOCITIES V(J) (EIGEN VALUES) C DO 10 J=1,3 V(J)=0.0 DO 10 I=1,3 P(I,J)=0.0 10 CONTINUE C C IMSL SUBROUTINES EVCSF AND EPIS C CALL DEVCSF(3,A,3,V,P,3) C C NUMERICALLY CONVERGED PERFORMANCE INDEX, CONV C ( CONV.LT.1.0 ; EXCELLENT ) C ( CONV.GT.1.0 BUT LT.100.0 ; FAIR ) C ( CONV.GT.100.0 ; POOR ) C CONV = DEPISF(3,3,A,3,V,P,3) if(conv.gt.100)write(6,100)conv 100 format(80(1h!),/,'poor convergence, conv=',e11.4) c c normalize the eigenvectors p(i,j) c sum1=dsqrt(p(1,1)**2+p(2,1)**2+p(3,1)**2) sum2=dsqrt(p(1,2)**2+p(2,2)**2+p(3,2)**2) sum3=dsqrt(p(1,3)**2+p(2,3)**2+p(3,3)**2) if(sum1.lt.1.0e-15)write(6,110)sum1 if(sum2.lt.1.0e-15)write(6,120)sum2 if(sum3.lt.1.0e-15)write(6,130)sum3 110 format(80(1h@),'sum1=',e11.4) 120 format(80(1h@),'sum2=',e11.4) 130 format(80(1h@),'sum3=',e11.4) do 15 i=1,3 p(i,1)=p(i,1)/sum1 p(i,2)=p(i,2)/sum2 p(i,3)=p(i,3)/sum3 15 continue C C C ENERGY FLUX DIRECT. COSINES B(1-3,J) C CALCULATED FROM: 1. MATERIAL STIFFNESS C(I,J) C 2. DISPL. DIRECT. COSINES P(1-3,J) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * * C * * C * * C * * C * * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C DO 20 J=1,3 C P1=P(1,J) P2=P(2,J) P3=P(3,J) C L1=C(1,1)*P1*P1*N1+C(1,2)*P1*P2*N2+C(1,3)*P1*P3*N3 * +C(1,6)*(P1*P1*N2+2.*P1*P2*N1)+C(2,6)*P2*P2*N2+C(3,6)*P2*P3*N3 * +C(6,6)*(P1*P2*N2+P2*P2*N1)+C(4,5)*(P2*P3*N3+P3*P3*N2) * +C(5,5)*(P1*P3*N3+P3*P3*N1) L2=C(1,6)*P1*P1*N1+C(2,6)*(2.*P1*P2*N2+P2*P2*N1)+C(3,6)*P1*P3*N3 * +C(6,6)*(P1*P1*N2+P1*P2*N1)+C(1,2)*P1*P2*N1+C(2,2)*P2*P2*N2 * +C(2,3)*P2*P3*N3+C(4,4)*(P2*P3*N3+P3*P3*N2) * +C(4,5)*(P1*P3*N3+P3*P3*N1) L3=C(4,5)*(2.*P1*P2*N3+P1*P3*N2+P2*P3*N1) * +C(5,5)*(P1*P1*N3+P1*P3*N1)+C(4,4)*(P2*P2*N3+P2*P3*N2) * +C(1,3)*P1*P3*N1+C(2,3)*P2*P3*N2+C(3,3)*P3*P3*N3 * +C(3,6)*(P1*P3*N2+P2*P3*N1) C SUM2=L1*L1+L2*L2+L3*L3 RMAG=DSQRT(SUM2) B(1,J)=L1/RMAG B(2,J)=L2/RMAG B(3,J)=L3/RMAG C C PHASE VELOCITIES V(J) CALCULATED FROM EIGENVALUES C V(J)=DSQRT(V(J)/RHO) C 20 CONTINUE RETURN END C C####################################################################### SUBROUTINE DISC(CTOL,IW1,IW2,IW3,IPHI,ITHETA,OP,P,V,B) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * SUBROUTINE DISC: * C * TESTS FOR DISCONTINUITIES IN EIGEN-VECTORS AND REARRANGES * C * THE WAVE VELOCITY (EIGEN-VALUES) TO CORRESPOND TO WAVE * C * SURFACES WITH SMOOTH AND CONTINUOUS EIGEN-VECTORS. EXCEPT * C * IN THE PRINCIPAL MATERIAL PLANES, A SMOOTH AND CONTINUOUS * C * EIGEN-VECTOR FIELD EXISTS ON A SINGLE CONNECTED SURFACE * C * BUT WITH A DISCONTINUOUS (SLOPE UNDEFINED) GEOMETRY. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 FP(3,3),OP(3,3),P(3,3),V(3),B(3,3) C C THE DISCONTINUITY IN EIGEN-VECTORS(PARTICLE DISPL.DIRECT.COSINES) ARE C TESTED FOR BY CONSIDERING THE ORDER OF EIGEN-VALUES CALCULATED BY C SUBROUTINE FLUX WHICH NUMERICALLY CONVERGES ON EIGEN-VALUES (WAVE SPEEDS) C IN ASCENDING ORDER (SUBSCRICPT 1 ASSOCIATED WITH SLOWEST VELOCITY). C WHEN PHASE VELOCITY SURFACES INTERSECT WAVE IDENTITIES ARE INTERCHANGED C SUCH THAT THE SLOWER VELOCITY IS ALWAYS ASSIGNED A SUBSCRIPT WHICH IS C ONE LESS THEN THE SUBSCRIPT IDENTIFIED WITH THE FASTER VELOCITY. THIS C EXCHANGE OF WAVE IDENTITIES OCCURS ONLY IN THE PRINCIPAL MATERIAL PLANES C AND IS TESTED FOR BY MONITORING THE CONTINUITY IN THE EIGEN-VECTORS. C EVERYWHERE ELSE THE THE EIGEN-VECTORS ARE SMOOTH AND CONTINUOUS ON C CONTINUOUS WAVE SURFACES. C C TEST IF THIS THE FIRST PLANE C IF(IPHI.EQ.1)GO TO 28 C C TEST IF THIS IS THE FIRST POINT WITHIN THE ROTATED PLANE. C IF TRUE THEN SAVE THE EIGEN-VECTORS FROM THE PREVIOUS PLANE C AS OLD EIGEN-VECTORS FOR THE FOLLOWING TEST OF CONTINUITY C IF(ITHETA.GT.1)GO TO 30 DO 29 I=1,3 OP(I,1)=FP(I,1) OP(I,2)=FP(I,2) OP(I,3)=FP(I,3) 29 CONTINUE GO TO 30 C C TEST IF THIS IS THE FIRST POINT WITH IN THE ROTATED PLANE C IF ITHETA=1 (THETAR=0.0) C THEN BRANCH TO SAVE EIGEN-VECTORS C 28 IF(ITHETA.EQ.1)GO TO 17 C C INITIALIZE VARIABLES INDICATIVE OF A DISCONTINUITY IN EIGEN-VECTORS C ( 0=NO DISCONTINUITY )( 1= DISCONTINUITY EXIST) C 30 IT1=0 IT2=0 IT3=0 C C TEST FOR DISCONTINUITY IN EIGEN-VECTORS BY COMPARING C OLD (OP(I,J)) AND NEW (P(I,J)) VALUES OF EIGEN-VECTORS C C TEST FOR DISCONTINUITY IN FIRST EIGEN-VECTOR C W1=DABS(OP(1,1)*P(1,1)+OP(2,1)*P(2,1)+OP(3,1)*P(3,1)) IF(W1.LT.CTOL)IT1=1 C C TEST FOR DISCONTINUITY IN SECOND EIGEN-VECTOR C W2=DABS(OP(1,2)*P(1,2)+OP(2,2)*P(2,2)+OP(3,2)*P(3,2)) IF(W2.LT.CTOL)IT2=1 C C TEST FOR DISCONTINUITY IN THIRD EIGEN-VECTOR C W3=DABS(OP(1,3)*P(1,3)+OP(2,3)*P(2,3)+OP(3,3)*P(3,3)) IF(W3.LT.CTOL)IT3=1 C C ASSUME ONLY ONE POSSIBLE WAVE PAIR IDENTITY EXCHANGE CAN EXIST C FOR EACH INCREMENT IN THETAR (IANG TO IANG+1 ) AND TEST WHICH C WAVES HAVE EXCHANGED IDENTITIES. C ITT=IT1+IT2+IT3 C C IF NO DISCONTINUITY IN EIGEN-VECTORS EXIST C THEN BRANCH TO SAVE EIGEN-VECTORS C IF(ITT.EQ.0)GO TO 17 C C FOR EACH INCREMENT IN THETA ( IANG TO IANG+1 ) C TEST FOR NUMBER OF EIGEN-VECTOR DISCONTINUITIES C IF(ITT.EQ.1)GO TO 10 IF(ITT.EQ.2)GO TO 14 IF(ITT.EQ.3)GO TO 12 10 WRITE(6,11)IPHI,ITHETA,W1,W2,W3 WRITE(10,11)IPHI,ITHETA,W1,W2,W3 11 FORMAT(9H**ERROR**,' IDENTITY CHANGED FOR ONLY ONE WAVE',/, * ' PHI=',I3,' ITHETA=',I3,' (W1,W2,W3)=(',E12.5,')') STOP C C FOR (ITT=3) TEST FOR IDENTITY EXCHANGE FOR THREE WAVES. THIS CONDITION C OCCURS MOST FREQUENTLY WHEN EIGEN VECTORS ARE ALL AT 90 DEGREES C TO THE PRINCIPAL COORD.S FOR THE PURE MODE CONDITIONS C 12 WRITE(10,123)IPHI,ITHETA,W1,W2,W3 123 FORMAT(80(1H%),/,' THREE WAVE IDENTITY EXCHANGE',/, * ' PHI=',I3,' THETA=',I3,' (W1,W2,W3)=(',3E12.5,')') C C TEST FOR AN EXCHANGE OF SURFACES 1-2, 2-3, AND 3-1 C W1=DABS(OP(1,1)*P(1,2)+OP(2,1)*P(2,2)+OP(3,1)*P(3,2)) W2=DABS(OP(1,2)*P(1,3)+OP(2,2)*P(2,3)+OP(3,2)*P(3,3)) W3=DABS(OP(1,3)*P(1,1)+OP(2,3)*P(2,1)+OP(3,3)*P(3,1)) TW=W1+W2+W3 DW=DABS(3.0-TW) IF(DW.GT.0.01)GO TO 201 C C EXCHANGE 1-2, 2-3, 3-1 C WRITE(6,456) 456 FORMAT(' EXCHANGE WAVES 1-2, 2-3, AND 3-1') IF(IWENA(IPHI,ITHETA).EQ.1)GO TO 19 C I1=IW1 I2=IW2 I3=IW3 IF(I1.EQ.1)IW1=2 IF(I2.EQ.1)IW2=2 IF(I3.EQ.1)IW3=2 IF(I1.EQ.2)IW1=3 IF(I2.EQ.2)IW2=3 IF(I3.EQ.2)IW3=3 IF(I1.EQ.3)IW1=1 IF(I2.EQ.3)IW2=1 IF(I3.EQ.3)IW3=1 GO TO 17 C C TEST FOR AN EXCHANGE OF SURFACES 1-3, 2-1, AND 3-2 C 201 W1=DABS(OP(1,1)*P(1,3)+OP(2,1)*P(2,3)+OP(3,1)*P(3,3)) W2=DABS(OP(1,2)*P(1,1)+OP(2,2)*P(2,1)+OP(3,2)*P(3,1)) W3=DABS(OP(1,3)*P(1,2)+OP(2,3)*P(2,2)+OP(3,3)*P(3,2)) TW=W1+W2+W3 DW=DABS(3.0-TW) IF(DW.GT.0.01)GO TO 202 C C EXCHANGE 1-3, 2-1, 3-2 C WRITE(6,678) 678 FORMAT(' EXCHANGE WAVES 1-3, 2-1, 3-2') IF(IWENA(IPHI,ITHETA).EQ.1)GO TO 19 C I1=IW1 I2=IW2 I3=IW3 IF(I1.EQ.1)IW1=3 IF(I2.EQ.1)IW2=3 IF(I3.EQ.1)IW3=3 IF(I1.EQ.2)IW1=1 IF(I2.EQ.2)IW2=1 IF(I3.EQ.2)IW3=1 IF(I1.EQ.3)IW1=2 IF(I2.EQ.3)IW2=2 IF(I3.EQ.3)IW3=2 GO TO 17 C C ERROR IN LOGIC. PRINT RESULTS AND TERMINATE C 202 WRITE(6,203)IPHI,ITHETA,W1,W2,W3 WRITE(10,203)IPHI,ITHETA,W1,W2,W3 203 FORMAT(' ERROR IN LOGIC FOR DETERMINING IDENTITY EXCHANGE',/, * ' FOR ALL THREE WAVES OVER ONE INCREMENT IN THETA',/, * ' PHI=',I3,' THETA=',I3,' (W1,W2,W3)=(',3E12.5,')') STOP C C FOR (ITT=2) DETERMINE WAVE PAIR IDENTITY EXCHANGE C 14 WRITE(6,890)IPHI,ITHETA,W1,W2,W3 890 FORMAT(80(1H%),/,' WAVE PAIR IDENTITY EXCHANGE',/, *' PHI=',I3,' ITHETA=',I3,' (W1,W2,W3)=(',3E12.5,')') IF((IT1.EQ.1).AND.(IT2.EQ.1))GO TO 15 IF((IT2.EQ.1).AND.(IT3.EQ.1))GO TO 16 IF((IT3.EQ.1).AND.(IT1.EQ.1))GO TO 18 WRITE(6,888)IPHI,ITHETA,W1,W2,W3 WRITE(10,888)IPHI,ITHETA,W1,W2,W3 888 FORMAT(9H**ERROR**,' DATA FOLLOWED INCORRECT LOGICAL PATH',/, * ' IPHI=',I3,' ITHETA=',I3,' (W1,W2,W3)=(',3E12.4,')') STOP C C FIRST AND SECOND WAVES EXCHANGED IDENTITIES C 15 WRITE(6,135) 135 FORMAT(' EXHANGE WAVES 1-2') IF(IWENA(IPHI,ITHETA).EQ.1)GO TO 19 C I1=IW1 I2=IW2 I3=IW3 IF(I1.EQ.1)IW1=2 IF(I2.EQ.1)IW2=2 IF(I3.EQ.1)IW3=2 IF(I1.EQ.2)IW1=1 IF(I2.EQ.2)IW2=1 IF(I3.EQ.2)IW3=1 GO TO 17 C C SECOND AND THIRD WAVES EXCHANGED IDENTITIES C 16 WRITE(6,246) 246 FORMAT(' EXCHANGE WAVES 2-3') IF(IWENA(IPHI,ITHETA).EQ.1)GO TO 19 C I1=IW1 I2=IW2 I3=IW3 IF(I1.EQ.2)IW1=3 IF(I2.EQ.2)IW2=3 IF(I3.EQ.2)IW3=3 IF(I1.EQ.3)IW1=2 IF(I2.EQ.3)IW2=2 IF(I3.EQ.3)IW3=2 GO TO 17 C C FIRST AND THIRD WAVES EXCHANGED IDENTITIES C 18 WRITE(6,357) 357 FORMAT(' EXCHANGE WAVES 1-3') IF(IWENA(IPHI,ITHETA).EQ.1)GO TO 19 C I1=IW1 I2=IW2 I3=IW3 IF(I1.EQ.1)IW1=3 IF(I2.EQ.1)IW2=3 IF(I3.EQ.1)IW3=2 IF(I1.EQ.2)IW1=1 IF(I2.EQ.2)IW2=1 IF(I3.EQ.2)IW3=1 GO TO 17 C C INDICATE THAT WAVE EXCHANGES DO NOT OCCUR WITHIN PRINCIPAL PLANES C 19 WRITE(6,919) 919 FORMAT(' WAVE EXCHANGE DID NOT OCCUR BECAUSE,',/, * ' MODE CHANGE OCCURED WITHIN PRINCIPAL MATERIAL PLANE') C C SAVE EIGEN-VECTORS C 17 DO 24 I=1,3 OP(I,1)=P(I,1) OP(I,2)=P(I,2) OP(I,3)=P(I,3) 24 CONTINUE C C SAVE EIGEN-VECTORS FROM PREVIOUS ROTATED PLANE C IF(ITHETA.GT.1)GO TO 26 DO 27 I=1,3 FP(I,1)=P(I,1) FP(I,2)=P(I,2) FP(I,3)=P(I,3) 27 CONTINUE C C SEQUENCE THE EIGEN-VALUES, EIGEN-VECTORS AND FLUX DEVIATIONS C 26 WV1=V(IW1) WV2=V(IW2) WV3=V(IW3) V(1)=WV1 V(2)=WV2 V(3)=WV3 C DO 25 I=1,3 PI1=P(I,IW1) PI2=P(I,IW2) PI3=P(I,IW3) P(I,1)=PI1 P(I,2)=PI2 P(I,3)=PI3 WB1=B(I,IW1) WB2=B(I,IW2) WB3=B(I,IW3) B(I,1)=WB1 B(I,2)=WB2 B(I,3)=WB3 25 CONTINUE C RETURN END C C####################################################################### SUBROUTINE DEV(IPHI,ITHETA,PI,N,P,B,DELTA,ZETA) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * SUBROUTINE DEV: * C * CALCULATES THE FLUX AND PARTICLE DEVIATIONS FROM THE UNIT * C * VECTOR NORMAL TO THE PLANE WAVE. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 N(3),P(3,3),B(3,3),DELTA(3),ZETA(3) C C CALCULATE FLUX DEVIATION ANGLES, DELTA(J), BETWEEN C ENERGY PROPAGATION DIRECTION AND THE VECTOR NORMAL TO THE PLANE WAVE C DO 26 J=1,3 CDELTA=DABS(N(1)*B(1,J)+N(2)*B(2,J)+N(3)*B(3,J)) DIFF=CDELTA-1.0 IF((CDELTA.GT.1.0).AND.(DIFF.GT.0.00001))WRITE(6,10)IPHI,ITHETA * ,CDELTA,DIFF 10 FORMAT(80(1H#),/,' IPHI=',I3,' ITHETA=',I3,/, *5x,' CDELTA=',E10.3,' DIFF IN DELTA=',E10.3) IF((CDELTA.GT.1.0).AND.(DIFF.GT.0.00001))STOP IF((CDELTA.GT.1.0).AND.(DIFF.LE.0.00001))WRITE(10,10)IPHI,ITHETA * ,CDELTA,DIFF IF((CDELTA.GT.1.0).AND.(DIFF.LE.0.00001))CDELTA=1.0 DELTA(J)=DACOS(CDELTA)*180./PI 26 CONTINUE C C CALCULATE THE DISPLACEMENT DEVIATION ANGLES, ZETA(J), BETWEEN C THE PARTICLE DISPL.S AND THE VECTOR NORMAL TO THE PLANE WAVE C DO 30 J=1,3 CZETA=DABS(N(1)*P(1,J)+N(2)*P(2,J)+N(3)*P(3,J)) DIFF=CZETA-1.0 IF((CZETA.GT.1.0).AND.(DIFF.GT.0.00001))WRITE(6,20)IPHI,ITHETA * ,CZETA,DIFF 20 FORMAT(80(1H#),/,' IPHI=',I3,' ITHETA=',I3,/, *5x,' CZETA=',E10.3,' DIFF IN ZETA=',E10.3) IF((CZETA.GT.1.0).AND.(DIFF.GT.0.00001))STOP IF((CZETA.GT.1.0).AND.(DIFF.LE.0.00001))WRITE(10,20)IPHI,ITHETA * ,CZETA,DIFF IF((CZETA.GT.1.0).AND.(DIFF.LE.0.00001))CZETA=1.0 ZETA(J)=DACOS(CZETA)*180./PI 30 CONTINUE C RETURN END C C####################################################################### SUBROUTINE TLT(IPHI,ITHETA,PI,N,P,PHIR,TILT,LT) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * SUBROUTINE TLT: * C * CALCULATES THE POLARIZATIONS WITHIN THE ROTATED PLANE AND * C * THE TILT OF THE DISPLACEMENT OUT OF THE ROTATED PLANE. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C IMPLICIT REAL*8 (A-H,O-Z) REAL*8 N(3),P(3,3),TILT(3),LT(3) C C CALCULATE SHIFT IN PHIR FOR CALCULATIONS C CPHIR=PI/2.-PHIR C C CALCULATE ANGLE, TILT(J), OF THE DISPLACEMENT VECTOR (EIGEN-VECTOR) C FROM THE ROTATED PLANE THAT CONTAINS THE DISPLACEMENT VECTOR(EIGEN-VECTOR) C DO 100 J=1,3 CTILT=DABS(P(1,J)*DCOS(CPHIR)-P(2,J)*DCOS(PHIR)) DIFF=CTILT-1.0 IF((CTILT.GT.1.0).AND.(DIFF.GT.0.00001))WRITE(6,10)IPHI,ITHETA * ,CTILT,DIFF 10 FORMAT(80(1H#),/,' IPHI=',I3,' ITHETA=',I3,/, *5x,' CTILT=',E10.3,' DIFF IN TILT=',E10.3) IF((CTILT.GT.1.0).AND.(DIFF.GT.0.00001))STOP IF((CTILT.GT.1.0).AND.(DIFF.LE.0.00001))WRITE(10,10)IPHI,ITHETA * ,CTILT,DIFF IF((CTILT.GT.1.0).AND.(DIFF.LE.0.00001))CTILT=1.0 TILT(J)=90.0-180.0*DACOS(CTILT)/PI 100 CONTINUE C C CALCULATE DISPLACEMENT POLARIZATION LT(J) WITHIN THE ROTATED PLANE C DO 200 J=1,3 F1=P(3,J)*DCOS(PHIR) F2=P(3,J)*DCOS(CPHIR) F3=-P(1,J)*DCOS(PHIR)-P(2,J)*DCOS(CPHIR) FMAG=DSQRT(F1**2+F2**2+F3**2) if(fmag.lt.1.0e-6)lt(j)=0.0 if(fmag.lt.1.0e-6)go to 200 c c if(fmag.lt.1.0e-15)write(6,300)iphi,itheta,j,fmag c *,f1,f2,f3,p(3,j),p(2,j),p(1,j) c 300 format(' phi=',i3,' theta=',i3,' j=',i1,' fmag=',e11.4, c *' f1=',e11.4,' f2=',e11.4,' f3=',e11.4, c *' p3=',e11.4,' p2=',e11.4,' p1=',e11.4) c F1=F1/FMAG F2=F2/FMAG F3=F3/FMAG CLT=DABS(F1*N(1)+F2*N(2)+F3*N(3)) DIFF=CLT-1.0 IF((CLT.GT.1.0).AND.(DIFF.GT.0.001))WRITE(6,20)IPHI,ITHETA * ,CLT,DIFF 20 FORMAT(80(1H#),/,' IPHI=',I3,' ITHETA=',I3,/, *5x,' CLT=',E10.3,' DIFF IN LT=',E10.3) IF((CLT.GT.1.0).AND.(DIFF.GT.0.001))STOP IF((CLT.GT.1.0).AND.(DIFF.LE.0.001))WRITE(10,20)IPHI,ITHETA * ,CLT,DIFF IF((CLT.GT.1.0).AND.(DIFF.LE.0.001))CLT=1.0 LT(J)=90.0-180.0*DACOS(CLT)/PI IF(TILT(J).EQ.90.0)LT(J)=90.0 200 CONTINUE C RETURN END C C####################################################################### SUBROUTINE IN(RHO,C) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * SUBROUTINE IN: * C * READS MATERIAL DENSITY AND STIFFNESSES FROM THE INPUT FILE * C * AND STORES MATERIAL INFO. AS ALPHANUMERIC DATA IN LABELED * C * COMMON /K/ TO BE PLOTTED AS A HEADER ABOVE EACH PLOT. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C IMPLICIT REAL*8 (A-H,O-Z) CHARACTER MTYPE*24,ICOM*24,IDEN*20,ISU*10 COMMON /K/MTYPE,ICOM,IDEN,ISU REAL*8 C(6,6) C C ZERO ALL STIFFNESSES C DO 10 I=1,6 DO 10 J=1,6 C(I,J)=0.0 10 CONTINUE C C READ MATERIAL TYPE AND COMMENT FROM INPUT FILE READ(5,1)MTYPE,ICOM 1 FORMAT(/,A24,/,A24) C C READ MATERIAL DENSITY, UNITS OF STIFFNESSES AND C NO. OF INDEPENDENT STIFFNESS TERM FROM INPUT FILE READ(5,4)IDEN,RHO,ISU,NS 4 FORMAT(A20,E10.3,/,16X,A10,/,31X,I2) C C LOOP THRU AND READ NO. OF INDEPENDENT STIFFNESS TERMS C FROM INPUT FILE DO 6 II=1,NS READ(5,5)I,J,C(I,J) 5 FORMAT(1X,I1,I1,1X,E10.3) 6 CONTINUE READ(5,7)MS 7 FORMAT(30X,I2) C C LOOP THRU AND READ NO. OF REMAINING STIFFNESS TERMS C FROM INPUT FILE DO 8 II=1,MS READ(5,5)I,J,C(I,J) 8 CONTINUE C RETURN END C C####################################################################### SUBROUTINE OUT(RHO,C) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * SUBROUTINE OUT * C * WRITES MATERIAL TYPE, DENSITY, AND STIFFNESSES ON OUTPUT FILE * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C IMPLICIT REAL*8 (A-H,O-Z) CHARACTER MTYPE*24,ICOM*24,IDEN*20,ISU*10 COMMON /K/MTYPE,ICOM,IDEN,ISU DIMENSION C(6,6) C C WRITE MATERIAL TYPE, COMMENT, AND DENSITY ON OUTPUT FILE C WRITE(6,10)MTYPE,ICOM,RHO 10 FORMAT(1H1,' MATERIAL: ',a24,/,' COMMENT: ',a24,/, * ' DENSITY:',2X,E10.3,//,18X,'STIFFNESSES',/) c WRITE(6,10)MTYPE,ICOM,IDEN,RHO,ISU c 10 FORMAT(1H1,' MATERIAL: ',A24,/,' COMMENT: ',A24,/, c * 2X,A20,E10.3,//,18X,'STIFFNESSES',A10,/) C C WRITE STIFFNESSES C(I,J) ON OUTPUT FILE C WRITE(6,20) 20 FORMAT(15X,'C11 C12 C13 C14 C15 C16',/, * 15X,' C22 C23 C24 C25 C26',/, * 15X,' C33 C34 C35 C36',/, * 15X,' C44 C45 C46',/, * 15X,' SYMMETRIC C55 C56',/, * 15X,' C66') WRITE(6,21)(C(1,J),J=1,6) 21 FORMAT(1X,6(1X,E11.4)) WRITE(6,22)(C(2,J),J=2,6) 22 FORMAT(13X,5(1X,E11.4)) WRITE(6,23)(C(3,J),J=3,6) 23 FORMAT(25X,4(1X,E11.4)) WRITE(6,24)(C(4,J),J=4,6) 24 FORMAT(37X,3(1X,E11.4)) WRITE(6,25)C(5,5),C(5,6),C(6,6) 25 FORMAT(13X,' SYMMETRIC',24X,2(1X,E11.4),/,62X,E11.4) C RETURN END C C############################################################################### FUNCTION IWENA(IPHI,ITHETA) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * SUBROUTINE WENA:(W-AVE E-XCHANGE N-OT A-LLOWED) * C * TESTS FOR THE LOCATION OF THE DETECTED DISCONTINUITY IN THE * C * DISPLACEMENT DIRECITON COSINES (EIGEN-VECTORS). IF THE * C * DISCONTINUITY FALLS ALONG A FORBIDDEN LINE THEN SET THE * C * VALUE OF WENA TO 1 OTHERWISE SET THE VALUE TO 0. * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C COMMON /PAR/NPHI,NTHETA,NTOL,IREG C C INITIALIZE IWENA TO ZERO C IWENA=0 C C TEST FOR REGION TYPE (IREG=1, 90X90 DEG. GRID; IREG=4, 360X180 DEG. GRID) C IF(IREG.EQ.1)GO TO 5 IF(IREG.EQ.4)GO TO 6 WRITE(6,111)IREG 111 FORMAT(80(1H&),' IREG=',I1,'PROG. STOP, INCORRECT REGION TYPE') STOP C C TEST IF CALCULATION IS ALONG FORBIDDEN LINES OF A 90X90 DEG. GRID C 5 IF(IPHI.EQ.1)GO TO 10 IF(IPHI.EQ.NPHI)GO TO 10 IF(ITHETA.EQ.1)GO TO 20 IF(ITHETA.EQ.NTHETA)GO TO 20 RETURN 10 IF((ITHETA.GT.NTOL).AND.(ITHETA.LT.(NTHETA-NTOL)))IWENA=1 RETURN 20 IF((IPHI.GT.NTOL).AND.(IPHI.LT.(NPHI-NTOL)))IWENA=1 RETURN C C TEST IF CALCULATION IS ALONG FORBIDDEN LINES OF A 360X180 DEG. GRID C 6 IN1=(NPHI-1)/4+1 IN2=(IN1-1)*2+1 IN3=(IN1-1)*3+1 IF(IPHI.EQ.1)GO TO 30 IF(IPHI.EQ.IN1)GO TO 30 IF(IPHI.EQ.IN2)GO TO 30 IF(IPHI.EQ.IN3)GO TO 30 IF(IPHI.EQ.NPHI)GO TO 30 IF(ITHETA.EQ.1)GO TO 40 IF(ITHETA.EQ.IN1)GO TO 40 IF(ITHETA.EQ.NTHETA)GO TO 40 RETURN 30 IF((ITHETA.GT.NTOL).AND.(ITHETA.LT.(IN1-NTOL)))IWENA=1 IF((ITHETA.GT.(IN1+NTOL)).AND.(ITHETA.LT.(NTHETA-NTOL)))IWENA=1 RETURN 40 IF((IPHI.GT.NTOL).AND.(IPHI.LT.(IN1-NTOL)))IWENA=1 IF((IPHI.GT.(IN1+NTOL)).AND.(IPHI.LT.(IN2-NTOL)))IWENA=1 IF((IPHI.GT.(IN2+NTOL)).AND.(IPHI.LT.(IN3-NTOL)))IWENA=1 IF((IPHI.GT.(NTOL+IN3)).AND.(IPHI.LT.(NPHI-NTOL)))IWENA=1 RETURN C END