C 
C     P R O G R A M   V E L P L
C     *************************
C
C     PROGRAM VELPL IS DESIGNED FOR THE PLOTTING OF PLANE SECTIONS
C     OF SLOWNESS, PHASE VELOCITY AND GROUP VELOCITY  SURFACES
C     FOR GENERAL ANISOTROPIC MEDIA
C
C     ************************************************************
C
      CHARACTER*80 TEXT,FILEIN,FILEOU,FILE1
C
C     **************************************************
C
      LIN=5
      LOU=6
      LU=1
      FILEIN='velpl.dat'
      FILEOU='velpl.out'
      FILE1='lu3.out'
      WRITE(*,'(2A)') ' Specify names of input and output files',
     1' LIN, LOU, LU3: '
      READ(*,*) FILEIN,FILEOU,FILE1
      IF(FILE1.EQ.' ') GO TO 99
      OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD')
      OPEN(LOU,FILE=FILEOU,FORM='FORMATTED')
      OPEN(LU,FILE=FILE1,FORM='FORMATTED',STATUS='OLD')
C
C**************************************************
C
      CALL PLOTS(LDUM1,LDUM2,7)
C
      IRUN=0
      read(lin,100)text
      write(lou,100)text
      READ(lin,107)IPRINT,SC,SHIFT
      WRITE(lou,107)IPRINT,SC,SHIFT
C
      xlen=10.
      ylen=10.
      call plot(xlen,ylen,-3)
c
    1 continue
      do 3 i=1,3
      iplot=0
      read(lu,101)itype,npar,inum
      if(itype.eq.0)stop
      if(npar.eq.1)write(lou,103)
      if(npar.eq.2)write(lou,104)
      if(npar.eq.3)write(lou,105)
      do 2 k=1,inum
      read(lu,102)ddelta,xv,yv,zv
      v=xv*xv+yv*yv+zv*zv
      if(npar.le.2)then
        vv=1./v
        v=sqrt(vv)
      end if
      if(npar.eq.1)then
        if(k.eq.1.and.i.eq.1)vconst=vv
        xv=xv*vconst
        yv=yv*vconst
        zv=zv*vconst
      end if
      if(npar.eq.2)then
        xv=xv*vv
        yv=yv*vv
        zv=zv*vv
      end if
      if(npar.eq.3)then
        v=sqrt(v)
      end if
      if(iprint.ne.0)write(lou,106)ddelta,v,xv,yv,zv
      call symbol(xv,zv,.01,'.',0.,-1)
      iplot=1
      vaux=v
    2 continue
    3 continue
      call plot(shift,-yv-ylen,-3)
      go to 1
c
  100 FORMAT(A)
  101 FORMAT(16I5)
  102 FORMAT(4E15.5)
  103 FORMAT(5x,'ANGLE',5x,'SLOW',8x,'PX',8x,'PY',8x,'PZ')
  104 FORMAT(5x,'ANGLE',5x,'PHVEL',7x,'VPX',7x,'VPY',7x,'VPZ')
  105 FORMAT(5x,'ANGLE',5x,'GRVEL',7x,'VGX',7x,'VGY',7x,'VGZ')
  106 format(5f10.5)
  107 FORMAT(I5,2F10.5)
c
   99 call plot(0.,0.,999)
      stop
      end
C
C=======================================================================
C
      INCLUDE 'calcops.for'
C     calcops.for
C
C=======================================================================
C