C
      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(4),NINT,
     1   XINTA
C
      READ(LU,101) MPRINT,NINT
      IF(MPRINT.GE.0)WRITE(LOU,101) MPRINT,NINT
      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,101) MX,MY
      IF(MPRINT.GE.0) 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,102)(SX(J),J=MX1,MX2)
      READ(LU,102)(SY(J),J=MY1,MY2)
      IF(MPRINT.GE.0) WRITE(LOU,102)(SX(J),J=MX1,MX2)
      IF(MPRINT.GE.0) WRITE(LOU,102)(SY(J),J=MY1,MY2)
      M1=MXY1
      DO 12 L=1,MX
      M2=M1+MY-1
      READ(LU,102)(Z(J),J=M1,M2)
      IF(MPRINT.GE.0) 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,102) ZMIN,ZMAX
      IF(MPRINT.GE.0) WRITE(LOU,102) ZMIN,ZMAX
      IF(MPRINT.EQ.0)GO TO 50
C
C  NUMERICAL FORM OF 3-D INTERFACES
C
      ZMM=ZMAX-ZMIN
      ZMMM=ZMM/10.
      IF(MPRINT.GE.0) WRITE(LOU,102) ZMIN,ZMAX,ZMMM
      DY=(BRD(4)-BRD(3))/50.
      DX=(BRD(2)-BRD(1))/5.
      Y(2)=BRD(3)-DY
      K1=1
      AUX(1)=BRD(1)
      DO 20 I=2,6
      AUX(I)=AUX(I-1)+DX
   20 CONTINUE
      IF(MPRINT.GE.1) WRITE(LOU,103)(AUX(I),I=1,6)
      DX=(BRD(2)-BRD(1))/100.
      MAUX=0
      NDER=0
      DO 29 K=1,51
      Y(2)=Y(2)+DY
      Y(1)=BRD(1)-DX
      DO 28 L=1,101
      Y(1)=Y(1)+DX
      CALL DISC(Y,DEP)
      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(K1.EQ.K.AND.MPRINT.GE.1) WRITE(LOU,104)Y(2),IPR
      IF(K1.NE.K.AND.MPRINT.GE.1) WRITE(LOU,105)IPR
      IF(K1.EQ.K)K1=K1+10
  29  CONTINUE
C
C  END OF LOOP OVER ALL INTERFACES
C
  50  CONTINUE
C
C  INPUT OF ELASTIC PARAMETERS
C
      READ(LU,101)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
        READ(LU,102)(RHO(L),L=1,NLAY)
        IF(MPRINT.GE.0)WRITE(LOU,102)(RHO(L),L=1,NLAY)
        DO 51 L=1,NLAY
        IF(ABS(RHO(L)).LT..00001)RHO(L)=1.
   51   CONTINUE
      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
      READ(LU,'(I10,6F10.4)')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,111)((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.GT.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.GT.1) WRITE(LOU,111)((A66U(J,K),J=1,6),K=1,6)
        END IF
        READ(LU,111)((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.GT.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.GT.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,111) A66U(1,1),A66U(4,4)
        READ(LU,111) 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.GT.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(5X,5(F7.3,13X),F7.3)
  104 FORMAT(F7.3,1X,101I1)
  105 FORMAT(8X,101I1)
  107 FORMAT(////,3X,A///)
  109 FORMAT(///,' INTERFACE NUMBER ',I5)
  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 ',/,' MATRIX OF
     1 ELASTIC PARAMETERS GIVEN IN (KM/SEC)**2:
     2  IT CONTAINS ELASTIC PARAMETERS/DENSITY',
     3/,' MATRIX IS SPECIFIED IMMEDIATELY BELOW THE UPPER
     4 BOUNDARY OF THE LAYER'/)
  116 FORMAT(/' ROTATION APPLIED AROUND X1 WITH FI1=',F10.5,/,18X,'ARO
     2UND X2 WITH FI2=',F10.5,/,18X,'AROUND X3 WITH FI3=',F10.5,/)
  117 FORMAT(/' MATRIX OF ELASTIC PARAMETERS GIVEN IN (KM/SEC)**2:
     1IT CONTAINS ELASTIC PARAMETERS/DENSITY',
     2/,' MATRIX IS SPECIFIED IMMEDIATELY
     3 ABOVE THE LOWER BOUNDARY OF THE LAYER'/)
  118 FORMAT(/' LAYER',I4,' IS ISOTROPIC'/)
  119 FORMAT(' VP**2=',F10.5,' (KM/SEC)**2      VS=',F10.5,' (KM/SEC)**2
     1 IMMEDIATELY BELOW THE UPPER BOUNDARY OF THE LAYER,',/,
     2' VP**2=',F10.5,' (KM/SEC)**2      VS=',F10.5,' (KM/SEC)*'
     3,'*2 IMMEDIATELY ABOVE THE LOWER BOUNDARY OF THE LAYER')
      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,400),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,400)
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,JJJ
  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
      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
C