      SUBROUTINE POLAR (N1,N2,NN,I)
C
C   ROUTINE COMPUTING POLARIZATION VECTORS
C
C   IF(IPOL.EQ.0) POLARIZATION VECTORS ARE NOT PRINTED
C   IF(IPOL.EQ.1) POLARIZATION VECTORS ARE PRINTED ONLY
C   FOR THE POINTS OF INCIDENCE OF A RAY AT INTERFACES
C   IF(IPOL.EQ.2) POLARIZATION VECTORS ARE PRINTED FOR ALL
C   COMPUTED POINTS ALONG THE RAY
C   N1,N2 - SUCCESSIVE NUMBERS OF THE FIRST AND THE LAST POINT
C   OF AN ELEMENT OF A RAY (IN CASE IPOL.EQ.0, N1=N2 AND THIS
C   NUMBER CORRESPONDS TO THE POINT OF INCIDENCE OF A RAY AT
C   AN INTERFACE; NEGATIVE N2 - INDICATION OF COMPOSED ELEMENT 
C   OF THE RAY
C   NN - TOTAL NUMBER OF POINTS ALONG THE RAY
C   I - NUMBER OF THE ELEMENT OF THE RAY
C
C   CALLED FROM ROUTINES: AMPL
C   ROUTINES CALLED:      CHRM
C
      DIMENSION PXN(3,3),Y(6),B1(3),B2(3),XX(3),E(21),E1(3),E2(3),T(3)
      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 /AUXI/  IANI(20),INTR,INT1,IOUT,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 /DJK/   D11,D12,D13,D22,D23,D33,DTR
      common /gam/   C11,C12,C13,C22,C23,C33
      COMPLEX PS
      COMMON /RAY/   AY(28,400),DS(20,20),KINT(20),HHH(3,3),tmax,
     1               PS(3,7,20),IS(8,20),DINC,DTOLD,N,IREF,IND,IND1
      COMMON /ZERO/ RNULL
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)
C
      IF(N2.LT.0)THEN
        LAY=IS(3,I-1)
        N2=-N2
      ELSE
        LAY=IS(3,I)
      END IF
      IF(ipol.GT.0)
     1WRITE(lou,'(A,6I5)')' LAY,I,N1,N2,ITYPE,IANI',LAY,I,N1,N2,ITYPE,
     2IANI(LAY)
      ns=n2-n1
      if(ns.eq.0)ns=1
      if(ipol.eq.2)ns=1
      IF(IANI(LAY).NE.0) GO TO 2000
C
C  ISOTROPIC CASE
C
      IF(NN.EQ.N2) GO TO 1005
      IF(ITYPE.NE.3) GO TO 1020
C
C  P-WAVE POLARIZATION VECTOR
C
 1005 DO 1010 J=N1,N2,ns
      V=1./SQRT(AY(5,J)*AY(5,J)+AY(6,J)*AY(6,J)+AY(7,J)*AY(7,J))
      t(1)=V*AY(5,J)
      t(2)=V*AY(6,J)
      t(3)=V*AY(7,J)
      if(ipol.ne.0)write(lou,'(a,/,6f10.5)') ' x1,x2,x3,tx,ty,tz',
     1ay(2,j),ay(3,j),ay(4,j),t
 1010 CONTINUE
      DO 1015 K=1,3
      HHH(1,K)=0.
      HHH(2,K)=0.
      HHH(3,K)=t(K)
 1015 CONTINUE
      IF(N2.EQ.NN) GO TO 1020
      RETURN
C
C  S-WAVE POLARIZATION VECTORS
C
 1020 continue
      DO 1100 J=N1,N2
      M1=5
      M2=6
      M3=7
      DVS1=AY(9,J)
      DVS2=AY(10,J)
 1030 P1=AY(M1,J)
      P2=AY(M2,J)
      P3=AY(M3,J)
      XL=SQRT(P1*P1+P2*P2)
      IF(XL.LT.RNULL)THEN
        IF(M1.EQ.4)THEN
          M1=6
          M2=7
          M3=5
          DVS1=AY(10,J)
          DVS2=AY(11,J)
          GO TO 1030
        END IF
        IF(M1.EQ.5)THEN
          M1=7
          M2=5
          M3=6
          DVS1=AY(11,J)
          DVS2=AY(9,J)
          GO TO 1030
        END IF
      END IF
C
C  COMPUTATION OF ANGLE XI SPECIFYING POLARIZATION VECTORS
C
      IF(ipol.Eq.2)WRITE(lou,'(A,/,6F12.6)')
     1' P1,P2,P3,XL,D1,D2',P1,P2,P3,XL,DVS1,DVS2
      IF(J.EQ.N1)THEN
        XO=(P3*(DVS1*P2-DVS2*P1))/(XL*XL)
        XI=0.
        GO TO 1040
      END IF
      Time=AY(1,J)-AY(1,J-1)
      X1=(P3*(DVS1*P2-DVS2*P1))/(XL*XL)
      XI=.5*(X1+XO)*Time+XI
 1040 IF(ipol.eq.2)WRITE(lou,'(A,/,6F12.6)')
     1' XI0,XI1,TSTEP,XISTEP',XO,X1,Time,XI
      XO=X1
      IF(IPOL.EQ.0.AND.J.NE.N2) GO TO 1100
      if(ipol.eq.1.and.(j.ne.n1.and.j.ne.n2))go to 1100
      P1=AY(5,J)
      P2=AY(6,J)
      P3=AY(7,J)
      XL=SQRT(P1*P1+P2*P2)
      if(xl.lt..000001)then
        b1(1)=1.
        b1(2)=0.
        b1(3)=0.
        b2(1)=0.
        b2(2)=1.
        if(p3.lt.0.)b2(2)=-1.
        b2(3)=0.
      else
        U=SQRT(P1*P1+P2*P2+P3*P3)
        P13=P1*P3/xl
        P23=P2*P3/xl
        B1(1)=P13/u
        B1(2)=P23/u
        B1(3)=-xl/u
        B2(1)=-P2/xl
        B2(2)=P1/xl
        B2(3)=0.
      end if
      SN=SIN(XI)
      CS=COS(XI)
      DO 1050 K=1,3
      E1(K)=CS*B1(K)-SN*B2(K)
      E2(K)=SN*B1(K)+CS*B2(K)
 1050 CONTINUE
      if(ipol.ne.0)write(lou,'(a,/,4f10.5/6f10.5)')
     1' TIME,X1,X2,X3/  E1X,E1Y,E1Z,E2X,E2Y,E2Z',(ay(k,j),k=1,4),e1,e2
 1100 CONTINUE
      DO 1110 K=1,3
      HHH(1,K)=E1(K)
      HHH(2,K)=E2(K)
      IF(N2.EQ.NN) GO TO 1110
      HHH(3,K)=0.
 1110 CONTINUE
      RETURN
C
C  ANISOTROPIC CASE
C
C
C  COMPUTATION OF POLARIZATION VECTORS
C
 2000 DO 2300 J=N1,N2,ns
      ISINGU=0
      DO 2110 K=1,21
      E(K)=AY(K+7,J)
 2110 CONTINUE
      DO 2120 K=1,6
      Y(K)=AY(K+1,J)
 2120 CONTINUE
      A2546=A25+A46
      A1266=A12+A66
      A1355=A13+A55
      A1456=A14+A56
      A3645=A36+A45
      A2344=A23+A44
      CALL CHRM(Y)
 2125 if(abs(dtr).lt..000001)go to 2135
      PXN(1,1)=D11
      PXN(1,2)=D12
      PXN(1,3)=D13
      XX(1)=D11*D11+D12*D12+D13*D13
      PXN(2,1)=D12
      PXN(2,2)=D22
      PXN(2,3)=D23
      XX(2)=D12*D12+D22*D22+D23*D23
      PXN(3,1)=D13
      PXN(3,2)=D23
      PXN(3,3)=D33
      XX(3)=D13*D13+D23*D23+D33*D33
      XN=0.
      DO 2130 L=1,3
      if(xn.ge.xx(l))go to 2130
      xn=xx(l)
      nx=l
 2130 CONTINUE
      xn=sqrt(xn)
      if(xn.gt.rnull)go to 2140
 2135 continue
      WRITE(lou,'(A,A,5I5)')
     1' XN.LT.RNULL IN POLAR  -  SHEAR-WAVE SINGULARITY ???',
     2'  LAY,N,N1,N2',LAY,J,J,N1,N2
      IND=10
      RETURN
 2140 P1=PXN(NX,1)/xn
      P2=PXN(NX,2)/xn
      P3=PXN(NX,3)/xn
C
C  CHECK OF PRECISION
C
      Z01=(C11-1.)*P1+C12*P2+C13*P3
      Z02=C12*P1+(C22-1.)*P2+C23*P3
      Z03=C13*P1+C23*P2+(C33-1.)*P3
      Z01=ABS(Z01)
      Z02=ABS(Z02)
      Z03=ABS(Z03)
      IF(Z01.GT.0.01.OR.Z02.GT.0.01.OR.Z03.GT.0.01) THEN
        WRITE(lou,'(A,4I5)') ' CHRISTOFFEL EQUATION IS SATISFIED
     1  WITH PRECISION LESS THAN 0.01: LAY,ITYPE,RAY ELEMENT,NPOINT ',
     2  LAY,ITYPE,I,J
      END IF
      if(ipol.ne.0)write(lou,'(a,/,6f10.5)') ' x1,x2,x3,gx,gy,gz',
     1ay(2,j),ay(3,j),ay(4,j),p1,p2,p3
      if(j.eq.n2)then
        HHH(ITYPE,1)=P1
        HHH(ITYPE,2)=P2
        HHH(ITYPE,3)=P3
      end if
 2300 CONTINUE
      RETURN
      END