SUBROUTINE INIT(pn,V)
C
C  ROUTINE FOR THE CALCULATION OF THE PHASE VELOCITIES FOR THE
C  SPECIFIED UNIT NORMAL TO THE PHASE FRONT 
C  SOLUTION OF THE EIGENVALUE PROBLEM
C
      DIMENSION C(3,3),pn(3),V(3)
C
      CALL CHRM1(C,pn,pn)
C
C     COMPUTATION OF THE ROOTS OF THE CHARACTERISTIC EQUATION WITH
C     CARDANS FORMULA
C
C     COMPUTATION OF THE COEFFICIENTS OF THE CHARACTERISTIC EQUATION
C     X**3+R*X**2+S*X+T=0
C
      c11=c(1,1)
      c12=c(1,2)
      c13=c(1,3)
      c22=c(2,2)
      c23=c(2,3)
      c33=c(3,3)
      a1=c22*c33-c23*c23
      a2=c11*c33-c13*c13
      a3=c11*c22-c12*c12
      R=-(C11+C22+C33)
      S=a1+a2+a3
      a2=c12*c33-c13*c23
      a3=c12*c23-c13*c22
      T=-c11*a1+c12*a2-c13*a3
C
C     SOLUTION OF CUBIC EQUATION FOLLOWING THE ALGORITHM
C     FROM NUMERICAL RECIPES     
C
      q=(r*r-3.*s)/9.
      p=(2.*r*r*r-9.*r*s+27.*t)/54.
      a1=q*q*q-p*p
      pi=3.14159265
      if(a1.le.0.)then
        th=pi
      end if
      if(a1.gt.0.)then
        th=q*q*q
        th=p/sqrt(th)
        th=acos(th)
      end if
      a1=-2.*sqrt(q)
      a2=-r/3.
      a3=th/3.
      v(1)=a1*cos(a3)+a2
      a3=(th+2.*pi)/3.
      v(2)=a1*cos(a3)+a2
      a3=(th+4.*pi)/3.
      v(3)=a1*cos(a3)+a2
      RETURN
      END