SUBROUTINE MODEL (MTEXT,LU) C C APPROXIMATION OF INTERFACES AND VELOCITY DISTRIBUTION C IN INDIVIDUAL LAYERS (ISOVELOCITY DISCONTINUITIES) C CHARACTER*80 MTEXT DIMENSION A66U(6,6),A66L(6,6),ANGU(3),ANGL(3) 1 ,AUX(12),DEP(6),Y(2),IPR(101) COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori COMMON /AUXX/ MMX(20),MMY(20),MMXY(20) COMMON /EPAR/ E66U(6,6,20),E66L(6,6,20) COMMON /DENS/ RHO(20) COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT, 1 XINTA COMMON/VRML/LUBRD,LUGRD,LUIND,LURAY C MPRINT=0 NINT=2 N1=10 N2=10 READ(LU,*) MPRINT,NINT,N1,N2 WRITE(LOU,101) MPRINT,NINT,N1,N2 IF(LUGRD.NE.0)THEN CALL FACETS(N1,N2,NINT) WRITE(LUGRD,108) WRITE(LUGRD,110) END IF NLAY=NINT-1 C C INPUT FOR 3D INTERFACES C MX2=0 MY2=0 MXY2=0 DO 13 I=1,NINT MX1=MX2+1 MY1=MY2+1 MXY1=MXY2+1 READ(LU,*) MX,MY WRITE(LOU,101) MX,MY NX(I)=MX NY(I)=MY MX2=MX1+MX-1 MY2=MY1+MY-1 MXY2=MXY1+MX*MY-1 READ(LU,*)(SX(J),J=MX1,MX2) READ(LU,*)(SY(J),J=MY1,MY2) WRITE(LOU,102)(SX(J),J=MX1,MX2) WRITE(LOU,102)(SY(J),J=MY1,MY2) M1=MXY1 DO 12 L=1,MX M2=M1+MY-1 READ(LU,*)(Z(J),J=M1,M2) WRITE(LOU,102)(Z(J),J=M1,M2) 12 M1=M2+1 C C DETERMINATION OF COEFFICIENTS OF BICUBIC SPLINES C APPROXIMATING INTERFACES C CALL BIAP(MX1,MX,MY1,MY,MXY1) MMX(I)=MX1 MMY(I)=MY1 MMXY(I)=MXY1 13 CONTINUE NB1=MMX(1) NB2=MMX(2)-1 BRD(1)=SX(NB1) BRD(2)=SX(NB2) NB1=MMY(1) NB2=MMY(2)-1 BRD(3)=SY(NB1) BRD(4)=SY(NB2) C C INPUT DATA FOR PRINTER PLOT OF 3-D INTERFACES C DO 50 INTR=1,NINT IF(MPRINT.GE.1) WRITE(LOU,109) INTR READ(LU,*) ZMIN,ZMAX IF(INTR.EQ.1)BRD(5)=ZMIN IF(INTR.EQ.NINT)BRD(6)=ZMAX IF(MPRINT.EQ.0)WRITE(LOU,102) ZMIN,ZMAX C C NUMERICAL FORM OF 3-D INTERFACES C ZMM=ZMAX-ZMIN ZMMM=ZMM/10. IF(MPRINT.GE.1) WRITE(LOU,103) ZMIN,ZMAX,ZMMM IF(MPRINT.GE.1) WRITE(LOU,104) BRD(1),BRD(2),BRD(3),BRD(4) DY=(BRD(4)-BRD(3))/FLOAT(N2-1) DX=(BRD(2)-BRD(1))/FLOAT(N1-1) MAUX=0 NDER=1 Y(2)=BRD(3)-DY DO 29 K=1,N2 Y(2)=Y(2)+DY Y(1)=BRD(1)-DX DO 28 L=1,N1 Y(1)=Y(1)+DX CALL DISC(Y,DEP) IF(LUGRD.NE.0)THEN DD=SQRT(1.+DEP(2)*DEP(2)+DEP(3)*DEP(3)) UN1=DEP(2)/DD UN2=DEP(3)/DD UN3=-1./DD KL=L+(K-1)*N1 WRITE(LUGRD,106)KL,Y(1),Y(2),DEP(1),UN1,UN2,UN3,INTR END IF AUX1=10.*(DEP(1)-ZMIN)/ZMM IPR(L)=IFIX(AUX1) IF(AUX1.LT.0.0.OR.AUX1.GT.10) IPR(L)=11 28 CONTINUE C C PRINTER PLOT OF 3-D INTERFACES C IF(MPRINT.GE.1) WRITE(LOU,105)(IPR(L),L=1,N1) 29 CONTINUE C C END OF LOOP OVER ALL INTERFACES C 50 CONTINUE IF(LUGRD.NE.0)WRITE(LUGRD,110) IF(LUGRD.NE.0)WRITE(LUGRD,110) C C INPUT OF ELASTIC PARAMETERS C ISQRT=0 IRHO=0 READ(LU,*)ISQRT,IRHO IF(MPRINT.EQ.0)WRITE(LOU,101)ISQRT,IRHO IF(ISQRT.EQ.0.AND.MPRINT.GT.0) WRITE(LOU,114) IF(ISQRT.NE.0.AND.MPRINT.GT.0) WRITE(LOU,113) IF(IRHO.NE.0)THEN DO 51 L=1,NLAY RHO(L)=1. 51 CONTINUE READ(LU,*)(RHO(L),L=1,NLAY) WRITE(LOU,102)(RHO(L),L=1,NLAY) END IF C C INPUT OF MATRIX OF ELASTIC PARAMETERS IN COMPRESSED FORM C ELEMENTS OF THE MATRIX HAVE TO BE GIVEN IN (KM**2/SEC**2) C WHICH CORRESPONDS TO ELASTIC PARAMETERS DIVIDED BY DENSITY C C IF THE CRYSTAL IS GIVEN IN OTHER COORDINATE SYSTEM THAN C THE COORDINATE SYSTEM IN WHICH RAY TRACING IS PERFORMED, C ROUTINE TRFMAT MAKES THE CORRESPONDING TRANSFORMATION C DO 30 L=1,NLAY IANI(L)=1 ANGU(1)=0. ANGU(2)=0. ANGU(3)=0. ANGL(1)=0. ANGL(2)=0. ANGL(3)=0. READ(LU,*)IANI(L),ANGU(1),ANGU(2),ANGU(3), 1ANGL(1),ANGL(2),ANGL(3) IF(MPRINT.EQ.0) 1WRITE(LOU,'(I10,6F10.4)')IANI(L),ANGU(1),ANGU(2),ANGU(3), 2ANGL(1),ANGL(2),ANGL(3) IROT1=1 IROT2=1 IF(ABS(ANGU(1)+ANGU(2)+ANGU(3)).LT.0.001) IROT1=0 IF(ABS(ANGL(1)+ANGL(2)+ANGL(3)).LT.0.001) IROT2=0 IF(IANI(L).NE.0) THEN READ(LU,*)((A66U(J,K),J=1,6),K=1,6) IF(MPRINT.EQ.0) 1 WRITE(LOU,111)((A66U(J,K),J=1,6),K=1,6) DO 55 K=1,6 DO 55 J=1,6 A66U(K,J)=A66U(J,K) 55 CONTINUE IF(MPRINT.GE.1) THEN WRITE(LOU,115) L WRITE(LOU,111)((A66U(J,K),J=1,6),K=1,6) END IF IF(IROT1.NE.0) THEN CALL TRFMAT(A66U,ANGU) IF(MPRINT.GE.1)WRITE(LOU,116)(ANGU(K),K=1,3) IF(MPRINT.GE.1) WRITE(LOU,111)((A66U(J,K),J=1,6),K=1,6) END IF READ(LU,*)((A66L(J,K),J=1,6),K=1,6) IF(MPRINT.EQ.0) 1 WRITE(LOU,111)((A66L(J,K),J=1,6),K=1,6) DO 60 K=1,6 DO 60 J=1,6 A66L(K,J)=A66L(J,K) 60 CONTINUE IF(MPRINT.GE.1) THEN WRITE(LOU,117) WRITE(LOU,111)((A66L(J,K),J=1,6),K=1,6) END IF IF(IROT2.NE.0) THEN CALL TRFMAT(A66L,ANGL) IF(MPRINT.GE.1)WRITE(LOU,116)(ANGL(K),K=1,3) IF(MPRINT.GE.1) WRITE(LOU,111)((A66L(J,K),J=1,6),K=1,6) END IF ELSE IF(MPRINT.GE.1)WRITE(LOU,118) L READ(LU,*) A66U(1,1),A66U(4,4) READ(LU,*) A66L(1,1),A66L(4,4) IF(MPRINT.EQ.0)WRITE(LOU,102) A66U(1,1),A66U(4,4) IF(MPRINT.EQ.0)WRITE(LOU,102) A66L(1,1),A66L(4,4) IF(MPRINT.GE.1) 1 WRITE(LOU,119)A66U(1,1),A66U(4,4),A66L(1,1),A66L(4,4) END IF DO 40 J=1,6 DO 40 K=1,6 E66U(K,J,L)=A66U(K,J) E66L(K,J,L)=A66L(K,J) 40 CONTINUE 30 CONTINUE WRITE(LOU,107)MTEXT C C FORMATS C 101 FORMAT(14I5) 102 FORMAT(8F10.5) 103 FORMAT(1X,'ZMIN,ZMAX,ZDIF',3F10.5) 104 FORMAT(1X,'(XMIN, XMAX)',2X,2F10.5,5X,'(YMIN, YMAX)',2F10.5) 105 FORMAT(8X,101I1) 106 FORMAT(5H'VRTX,I3,1H',1X,6F10.5,1X,I2,1X,'/') 107 FORMAT(////,3X,A///) 108 FORMAT(10H'VERTICES') 109 FORMAT(///,' INTERFACE NUMBER ',I5) 110 FORMAT('/') 111 FORMAT(6F10.5) 113 FORMAT(//' INTERPOLATION IN SQUARE ROOTS OF ELASTIC PARAMETERES'/) 114 FORMAT(//' INTERPOLATION IN VALUES OF ELASTIC PARAMETERES'/) 115 FORMAT(//' LAYER',I4,' IS ANISOTROPIC ',//,' DENSITY NORMALIZED MAT 1RIX OF ELASTIC PARAMETERS IN (KM/SEC)**2', 2/,' MATRIX IS SPECIFIED IMMEDIATELY BELOW THE UPPER BOUNDARY OF TH 3E LAYER'/) 116 FORMAT(/' ROTATION APPLIED AROUND X1 WITH FI1=',F10.5,/,18X,'AROUN 2D X2 WITH FI2=',F10.5,/,18X,'AROUND X3 WITH FI3=',F10.5,/) 117 FORMAT(/' DENSITY NORMALIZED MATRIX OF ELASTIC PARAMETERS IN (KM/S 1EC)**2',/,' MATRIX IS SPECIFIED IMMEDIATELY ABOVE THE LOWER BOUNDA 2RY OF THE LAYER'/) 118 FORMAT(//' LAYER',I4,' IS ISOTROPIC'/) 119 FORMAT(' IMMEDIATELY BELOW THE UPPER BOUNDARY OF THE LAYER' 1/' VP**2=',F7.5,' (KM/SEC)**2 VS=',F7.5,' (KM/SEC)**2'// 2' IMMEDIATELY ABOVE THE LOWER BOUNDARY OF THE LAYER' 2/' VP**2=',F7.5,' (KM/SEC)**2 VS=',F7.5,' (KM/SEC)**2') RETURN END C C ********************************************************* C SUBROUTINE PARDIS(Y,IAY) C SAVE Z1,DZX1,DZY1,DZXX1,DZXY1,DZYY1,INTR1 DIMENSION Y(18),DEP(6),B(21),E(21),EX(21),EY(21),EZ(21),EXX(21), 1 EYY(21),EZZ(21),EXY(21),EXZ(21),EYZ(21) COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 COMMON /APROX1/ D(21),DX(21),DY(21),DZ(21),DXX(21), 1 DXY(21),DXZ(21),DYY(21),DYZ(21),DZZ(21) COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /EPAR/ E66U(6,6,20),E66L(6,6,20) COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,20),KINT(20),HHH(3,3),TMAX, 1 PS(3,7,20),IS(8,20),N,IREF,IND,IND1 COMMON /RAY2/ DRY(3,2000) COMMON/DENS/RHO(20) C EQUIVALENCE(E(1),A11),(E(2),A12),(E(3),A13),(E(4),A14),(E(5),A15) 1 ,(E(6),A16),(E(7),A22),(E(8),A23),(E(9),A24) 2 ,(E(10),A25),(E(11),A26),(E(12),A33),(E(13),A34),(E(14),A35) 2 ,(E(15),A36),(E(16),A44),(E(17),A45),(E(18),A46),(E(19),A55) 2 ,(E(20),A56),(E(21),A66) 1 ,(EX(1),DXA11),(EX(2),DXA12),(EX(3),DXA13),(EX(4),DXA14) 1 ,(EX(5),DXA15),(EX(6),DXA16),(EX(7),DXA22),(EX(8),DXA23) 1 ,(EX(9),DXA24),(EX(10),DXA25),(EX(11),DXA26),(EX(12),DXA33) 1 ,(EX(13),DXA34),(EX(14),DXA35),(EX(15),DXA36),(EX(16),DXA44) 1 ,(EX(17),DXA45),(EX(18),DXA46),(EX(19),DXA55),(EX(20),DXA56) 1 ,(EX(21),DXA66) EQUIVALENCE 1 (EY(1),DYA11),(EY(2),DYA12),(EY(3),DYA13),(EY(4),DYA14) 1 ,(EY(5),DYA15),(EY(6),DYA16),(EY(7),DYA22),(EY(8),DYA23) 1 ,(EY(9),DYA24),(EY(10),DYA25),(EY(11),DYA26),(EY(12),DYA33) 1 ,(EY(13),DYA34),(EY(14),DYA35),(EY(15),DYA36),(EY(16),DYA44) 1 ,(EY(17),DYA45),(EY(18),DYA46),(EY(19),DYA55),(EY(20),DYA56) 1 ,(EY(21),DYA66) 1 ,(EZ(1),DZA11),(EZ(2),DZA12),(EZ(3),DZA13),(EZ(4),DZA14) 1 ,(EZ(5),DZA15),(EZ(6),DZA16),(EZ(7),DZA22),(EZ(8),DZA23) 1 ,(EZ(9),DZA24),(EZ(10),DZA25),(EZ(11),DZA26),(EZ(12),DZA33) 1 ,(EZ(13),DZA34),(EZ(14),DZA35),(EZ(15),DZA36),(EZ(16),DZA44) 1 ,(EZ(17),DZA45),(EZ(18),DZA46),(EZ(19),DZA55),(EZ(20),DZA56) 1 ,(EZ(21),DZA66) C INTR1=INTR INTR=LAY CALL DISC(Y,DEP) Z1=DEP(1) DZX1=DEP(2) DZY1=DEP(3) IF(NDER.LE.1) GOTO 10 DZXX1=DEP(4) DZXY1=DEP(5) DZYY1=DEP(6) 10 INTR=LAY+1 CALL DISC(Y,DEP) INTR=INTR1 Z2=DEP(1) DZX2=DEP(2) DZY2=DEP(3) IF(NDER.LE.1) GOTO 20 DZXX2=DEP(4) DZXY2=DEP(5) DZYY2=DEP(6) 20 AUX1=Z2-Z1 IF(AUX1.LE.0.) THEN IND=20 WRITE(LOU,'(A,5F12.6,2I5)')' Z1,Z2,X,Y,Z,LAY,IND',Z1,Z2,Y(1), 1Y(2),Y(3),LAY,IND RETURN END IF C C THE NEIGHBOURHOOD INTERFACES INTERSECT EACH OTHER C AUX2=Y(3)-Z1 AUX3=Y(3)-Z2 AUX4X=DZX2-DZX1 AUX4Y=DZY2-DZY1 AU2=1./AUX1/AUX1 AU3=AU2/AUX1 AU3X=DZX1*AUX3 AU3Y=DZY1*AUX3 AU4X=DZX2*AUX2 AU4Y=DZY2*AUX2 AU5X=AU4X-AU3X AU5Y=AU4Y-AU3Y A1=AUX2/AUX1 A2=-AU2*AU5X A3=-AU2*AU5Y A4=1./AUX1 IF(NDER.EQ.1) GOTO 30 A5=AU3*(2.*AUX4X*AU5X+AUX1*(DZXX1*AUX3-DZXX2*AUX2)) A6=AU3*(2.*AUX4Y*AU5X+AUX1*(DZXY1*AUX3-DZXY2*AUX2+DZX2*DZY1- 1DZX1*DZY2)) A7=-AU2*AUX4X A8=AU3*(2.*AUX4Y*AU5Y+AUX1*(DZYY1*AUX3-DZYY2*AUX2)) A9=-AU2*AUX4Y 30 JJ=21 JJJ=1 IF(ISQRT.NE.0.AND.IANI(LAY).EQ.0) GOTO 37 C C INTERPOLATION OF ELASTIC PARAMETERS DIVIDED BY DENSITY C (CORRESPONDS TO THE INTERPOLATION IN SQUARES OF VELOCITY) C IF(IANI(LAY).EQ.0) GOTO 33 J1=0 DO 31 L=1,6 DO 32 J=L,6 K1=J-L+1+J1 E(K1)=E66U(J,L,LAY) B(K1)=E66L(J,L,LAY)-E(K1) 32 CONTINUE J1=K1 31 CONTINUE GOTO 52 33 E(1)=E66U(1,1,LAY) B(1)=E66L(1,1,LAY)-E(1) E(16)=E66U(4,4,LAY) B(16)=E66L(4,4,LAY)-E(16) JJ=16 JJJ=15 GOTO 52 C C INTERPOLATION OF SQUARE ROOTS OF ELASTIC PARAMETERS C (CORRESPONDS TO THE INTERPOLATION OF VELOCITIES) C IT WORKS ONLY IN ISOTROPIC LAYERS C 37 E(1)=SQRT(E66U(1,1,LAY)) B(1)=SQRT(E66L(1,1,LAY))-E(1) E(16)=SQRT(E66U(4,4,LAY)) B(16)=SQRT(E66L(4,4,LAY))-E(16) JJ=16 JJJ=15 C C ELASTIC PARAMETERS AND THEIR DERIVATIVES OBTAINED BY C INTERPOLATION IN VELOCITIES. ELASTIC PARAMETERS C ARE OBTAINED AS SQUARES OF INTERPOLATED QUANTITIES C 40 DO 50 J=1,JJ,JJJ BB=B(J) C C ELASTIC PARAMETERS C E(J)=E(J)+A1*BB EE=2.*E(J) E(J)=E(J)*E(J) C C FIRST DERIVATIVES OF ELASTIC PARAMETERS C EX(J)=A2*BB EEX=EX(J) EX(J)=EX(J)*EE EY(J)=A3*BB EEY=EY(J) EY(J)=EY(J)*EE EZ(J)=A4*BB EEZ=EZ(J) EZ(J)=EZ(J)*EE D(J)=E(J) DX(J)=EX(J) DY(J)=EY(J) DZ(J)=EZ(J) IF(NDER.LE.1) GOTO 50 C C SECOND DERIVATIVES OF ELASTIC PARAMETERS C EXX(J)=A5*BB*EE+2.*EEX*EEX EXY(J)=A6*BB*EE+2.*EEX*EEY EXZ(J)=A7*BB*EE+2.*EEX*EEZ EYY(J)=A8*BB*EE+2.*EEY*EEY EYZ(J)=A9*BB*EE+2.*EEY*EEZ EZZ(J)=2.*EEZ*EEZ DXX(J)=EXX(J) DXY(J)=EXY(J) DXZ(J)=EXZ(J) DYY(J)=EYY(J) DYZ(J)=EYZ(J) DZZ(J)=EZZ(J) 50 CONTINUE GOTO 59 C C ELASTIC PARAMETERS AND THEIR DERIVATIVES OBTAINED BY C INTERPOLATION IN VALUES OF ELASTIC PARAMETERS C 52 DO 55 J=1,JJ,JJJ BB=B(J) C C ELASTIC PARAMETERS C E(J)=E(J)+A1*BB C C FIRST DERIVATIVES OF ELASTIC PARAMETERS C EX(J)=A2*BB EY(J)=A3*BB EZ(J)=A4*BB D(J)=E(J) DX(J)=EX(J) DY(J)=EY(J) DZ(J)=EZ(J) IF(NDER.LE.1) GOTO 55 C C SECOND DERIVATIVES OF ELASTIC PARAMETERS C EXX(J)=A5*BB EXY(J)=A6*BB EXZ(J)=A7*BB EYY(J)=A8*BB EYZ(J)=A9*BB EZZ(J)=0. DXX(J)=EXX(J) DXY(J)=EXY(J) DXZ(J)=EXZ(J) DYY(J)=EYY(J) DYZ(J)=EYZ(J) DZZ(J)=EZZ(J) 55 CONTINUE C 59 IF(IANI(LAY).EQ.0) GOTO 90 A2546=A25+A46 A1266=A12+A66 A1355=A13+A55 A1456=A14+A56 A3645=A36+A45 A2344=A23+A44 IF(IAY.EQ.0)RETURN DO 60 I=1,21 60 AY(I+7,N)=E(I) RETURN 90 IF(IAY.EQ.0)RETURN AY(8,N)=A11 AY(9,N)=DXA11 AY(10,N)=DYA11 AY(11,N)=DZA11 AY(12,N)=A44 AY(13,N)=DXA44 AY(14,N)=DYA44 AY(15,N)=DZA44 RO=1.7+0.2*SQRT(A11) IF(IRHO.NE.0)RO=RHO(LAY) AY(16,N)=RO RETURN END C C ********************************************************* C SUBROUTINE TRFMAT(A66,ANG) C C SUBROUTINE FOR THE TRANSFORMATION OF THE CRYSTAL FROM ITS OWN C COORDINATE SYSTEM TO THE COORDINATE SYSTEM USED BY THE RAY C TRACING PROGRAM C C ANG(1-3) ROTATION ANGLES C COMMON/ZERRO/RNULL DIMENSION A66(6,6),ANG(3),A1N(3,3,3,3),D(3,3), 1DA1(3,3,3,3) C C COMPUTATION OF THE MATRIX A1N (MATRIX OF ELASTIC COEFICIENTS C IN ITS OWN COORDINATE SYSTEM) C DO 30 I=1,3 DO 30 J=1,3 DO 30 K=1,3 DO 30 L=1,3 CALL INDEX1(I,J,I1) CALL INDEX1(K,L,K1) A1N(I,J,K,L)=A66(I1,K1) 30 CONTINUE C C COMPUTATION OF MATRIX D (THE MATRIX FOR TRANSFORMATION) C CALL TRANSF(ANG,D) C C COMPUTATION OF THE MATRIX OF ELASTIC COEFFICIENTS C IN THE COORDINATE SYSTEM FOR RAY COMPUTATION C DO 40 I=1,3 DO 40 N=1,3 DO 40 IR=1,3 DO 40 IS=1,3 DA1(I,N,IR,IS)=0.0 DO 40 M=1,3 DA1(I,N,IR,IS)=DA1(I,N,IR,IS)+D(I,M)*A1N(M,N,IR,IS) 40 CONTINUE C DO 41 I=1,3 DO 41 J=1,3 DO 41 IR=1,3 DO 41 IS=1,3 A1N(I,J,IR,IS)=0.0 DO 41 N=1,3 A1N(I,J,IR,IS)=A1N(I,J,IR,IS)+D(J,N)*DA1(I,N,IR,IS) 41 CONTINUE C DO 43 I=1,3 DO 43 J=1,3 DO 43 K=1,3 DO 43 IS=1,3 DA1(I,J,K,IS)=0.0 DO 43 IR=1,3 DA1(I,J,K,IS)=DA1(I,J,K,IS)+D(K,IR)*A1N(I,J,IR,IS) 43 CONTINUE C DO 44 I=1,3 DO 44 J=1,3 DO 44 K=1,3 DO 44 L=1,3 A1N(I,J,K,L)=0.0 DO 44 IS=1,3 A1N(I,J,K,L)=A1N(I,J,K,L)+D(L,IS)*DA1(I,J,K,IS) 44 CONTINUE C C C COMPUTATION OF THE MATRIX A66 (THE MATRIX A1N IN THE C COMPRESSED FORM) C DO 50 I=1,6 DO 50 J=1,6 CALL INDEX2(I,I1,J1) CALL INDEX2(J,K1,L1) A66(I,J)=A1N(I1,J1,K1,L1) IF(ABS(A66(I,J)).LT.RNULL) A66(I,J)=0. 50 CONTINUE RETURN END C C ************************************************************ C SUBROUTINE INDEX1(I,J,I1) C C SUBROUTINE FOR DETERMINING THE INDEX I1 FOR THE SYMETRIC C TENSOR OF SECOND RANK IN COMPRESSED FORM FROM TWO C INDICES I,J OF THE SAME TENSOR ELEMENT IN NONCOMPRESSED C FORM C C IF(I.NE.J) GO TO 10 I1=I RETURN C 10 CONTINUE IJ=I+J-2 GO TO (20,30,40),IJ C 20 CONTINUE I1=6 RETURN C 30 CONTINUE I1=5 RETURN C 40 CONTINUE I1=4 RETURN END C C *********************************************************** C SUBROUTINE INDEX2(I,I1,J1) C C SUBROUTINE FOR DETERMINING THE TWO INDEXES I1,J1 WHICH C CORRESPOND TO INDEX I OF THE ELEMENT OF THE SYMETRIC C TENSOR OF SECOND RANK IN THE COMPRESSED FORM C GO TO (10,10,10,20,30,40),I C 10 CONTINUE I1=I J1=I RETURN C 20 CONTINUE I1=2 J1=3 RETURN C 30 CONTINUE I1=1 J1=3 RETURN C 40 CONTINUE I1=1 J1=2 RETURN END C C *********************************************************** C SUBROUTINE TRANSF(AN,D) C C SUBROUTINE FOR COMPUTING THE MATRIX OF ROTATION FOR THE C CRYSTAL AXES, WHICH INITIALLY COINCIDE WITH THE MODEL AXES C (X1, X2 SITUATED IN THE HORIZONTAL PLANE, X3 VERTICAL), C INTO THEIR PROPER POSITION IN THE MODEL. C C THE MATRIX IS SPECIFIED BY THREE ANGLES, PHI=AN(1), DE=AN(2), C GA=AN(3). THE ANGLES PHI AND DE SPECIFY THE ORIENTATION OF C THE CRYSTAL AXIS X3 IN THE MODEL, THE ANGLE GA SPECIFIES C THE ORIENATION OF CRYSTAL AXES X1, X2 IN THE PLANE PERPEN- C DICULAR TO THE CRYSTAL AXIS X3. C PHI... AZIMUTH (IN DEGREES). ANGLE BETWEEN POSITIVE DIREC- C TIONS OF THE X1 MODEL AXIS AND HORIZONTAL PROJECTION C OF X3 CRYSTAL AXIS. IT IS POSITIVE IF MEASURED FROM C POSITIVE DIRECTION OF X1 AXIS TOWARDS POSITIVE DIREC- C TION OF X2 AXIS OF MODEL COORDINATES. C DE ... INCLINATION (0-90 DEGREES). ANGLE BETWEEN THE POSITI- C VE DIRECTION OF THE MODEL X3 AXIS AND X3 AXIS OF C CRYSTAL COORDINATE SYSTEM. C GA ... ROTATION ANGLE (IN DEGREES) IN THE PLANE PERPEN- C DICULAR TO THE CRYSTAL X3 AXIS. FOR GA=0, CRYSTAL C X2 AXIS IS HORIZONTAL, X1 AXIS IS PERPENDICULAR C TO IT AND POINTS DOWN. FOR GA NONZERO, THE CRYSTAL C X1, X2 AXES ARE ROTATED POSITIVELY FROM THEIR POSI- C TION FOR GA=0. C NOTE1: FOR PSI=0, DE=0 AND GA=0 CRYSTAL AXES COINCIDE WITH C AXES OF MODEL COORDINATES, I.E. POSITIVE CRYSTAL AXIS C POINTS DOWN. C NOTE2: BOTH MODEL AND CRYSTAL COORDINATE SYSTEMS ARE RIGHT- C HANDED. C NOTE3: COMPONENTS OF THE UNIT BASE VECTOR OF THE CRYSTAL C AXIS X3 IN MODEL COORDINATES ARE: C (COS(PHI)SIN(DE), SIN(PHI)SIN(DE), COS(DE)). C C DECLARATIONS C DIMENSION AN(3),D(3,3) C DATA PI180/0.0174532925/ C C C CONVERSION INTO RAD MODE C PHI=AN(1)*PI180 DE=AN(2)*PI180 GA=AN(3)*PI180 C C SF=SIN(PHI) CF=COS(PHI) SD=SIN(DE) CD=COS(DE) SG=SIN(GA) CG=COS(GA) C IF(ABS(GA).GE..000001)GO TO 10 C C C ONLY X3 AXIS OF CRYSTAL COORDINATE SYSTEM IS C SPECIFIED. ORIENTATION OF X1 AND X2 CRYSTAL AXES C CAN BE ARBITRARY - CASE OF MEDIUM WITH HEXAGONAL C SYMMETRY C D(1,1)=CF*CD D(1,2)=-SF D(1,3)=CF*SD D(2,1)=SF*CD D(2,2)=CF D(2,3)=SF*SD D(3,1)=-SD D(3,2)=0. D(3,3)=CD RETURN C C C ALL THREE AXES OF CRYSTAL COORDINATE SYSTEM ARE C SPECIFIED C 10 D(1,1)=-SF*SG+CF*CD*CG D(1,2)=-SF*CG-CF*CD*SG D(1,3)=CF*SD D(2,1)=CF*SG+SF*CD*CG D(2,2)=CF*CG-SF*CD*SG D(2,3)=SF*SD D(3,1)=-SD*CG D(3,2)=SD*SG D(3,3)=CD C RETURN END