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,PSTEXT,FILEIN,FILEOU,FILE1
      DATA RTOD/57.29578122/
C
C     **************************************************
C
      LIN=5
      LOU=6
      LU3=1
      FILEIN='velpl.dat'
      FILEOU='velpl.out'
      FILE1='lu3.dat'
      WRITE(*,'(2A)') ' (VELPL) SPECIFY NAMES OF INPUT AND OUTPUT',
     1' FILES 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(LU3,FILE=FILE1,FORM='FORMATTED',STATUS='OLD')
C
C**************************************************
C
      IRUN=0
      TEXT='VELPL'
      PSTEXT=' '
      IPLOT=0
      IPRINT=0
      SC=1.
      READ(LIN,*)TEXT
      WRITE(LOU,100)TEXT
      READ(LIN,*)IPLOT,IPRINT,SC,PSTEXT
      WRITE(LOU,107)IPLOT,IPRINT,SC,PSTEXT
      IF(IPRINT.LT.0)THEN
        CALL PLOTN(PSTEXT,0)
        IPRINT=-IPRINT
      END IF
      CALL PLOTS(LDUM1,LDUM2,7)
      IF(IPLOT.GE.0)THEN
        READ(LIN,*)XM,XLEN,DTICX,YM,YLEN,DTICY
        WRITE(LOU,106)XM,XLEN,DTICX,YM,YLEN,DTICY
        XMER=.5*XLEN/XM
        YMER=.5*YLEN/YM
      ELSE
        READ(LIN,*)XLEN,DTICX,YMIN,YMAX,YLEN,DTICY
        WRITE(LOU,106)XLEN,DTICX,YMIN,YMAX,YLEN,DTICY
        XMER=XLEN/90.
        YMER=YLEN/(YMAX-YMIN)
      END IF
      NTICX=1
      NTICY=1
      NDX=0
      NDY=0
      READ(LIN,*)NTICX,NTICY,NDX,NDY
      WRITE(LOU,101)NTICX,NTICY,NDX,NDY
C
      CALL PLOT(XLEN,YLEN,-3)
      IF(IPLOT.GE.0)THEN
        CALL BORDER(XLEN,DTICX,YLEN,DTICY,SC,TEXT,-1,-XM,XM,-YM,YM,
     1  NTICX,NTICY,NDX,NDY)
      ELSE
        CALL BORDER(XLEN,DTICX,YLEN,DTICY,SC,TEXT,-1,0.,90.,YMIN,YMAX,
     1  NTICX,NTICY,NDX,NDY)
       END IF
C
    1 CONTINUE
      DO 3 I=1,3
      READ(LU3,101)ITYPE,NPAR,INUM
      IF(ITYPE.EQ.0)GO TO 99
      IF(IPRINT.NE.0)THEN
        IF(NPAR.EQ.1)WRITE(LOU,103)
        IF(NPAR.EQ.2)WRITE(LOU,104)
        IF(NPAR.EQ.3)WRITE(LOU,105)
      END IF
      DO 2 K=1,INUM
      READ(LU3,102)DDELTA,XV,YV,ZV
      V=XV*XV+YV*YV+ZV*ZV
      HV=SQRT(XV*XV+YV*YV)
      IF(XV.LT.0.)HV=-HV
      V=SQRT(V)
      IF(NPAR.EQ.2)THEN
        IF(V.NE.0.)V=1./V
        HV=HV*V*V
        ZV=ZV*V*V
      END IF
      IF(IPRINT.NE.0)WRITE(LOU,106)DDELTA,V,XV,YV,ZV
      IF(IPLOT.GE.0)THEN
        XX=XLEN/2.+HV*XMER
        ZZ=YLEN/2.+ZV*YMER
        IF(ABS(HV*XMER).LE..5*XLEN.AND.ABS(ZV*YMER).LE..5*YLEN)
     *  CALL SYMBOL(XX,ZZ,.01,'.',0.,-1)
      ELSE
        DEG=DDELTA*RTOD
        XX=DEG*XMER
        ZZ=(V-YMIN)*YMER
        IF(XX.LE.XLEN.AND.ZZ.GE.0..AND.ZZ.LE.YLEN)
     *  CALL SYMBOL(XX,ZZ,.01,'.',0.,-1)
      END IF
    2 CONTINUE
    3 CONTINUE
      GO TO 1
C
  100 FORMAT(A)
  101 FORMAT(16I5)
  102 FORMAT(6E15.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(8F10.5)
  107 FORMAT(2I5,F10.5,1X,A)
c
   99 CALL PLOT(0.,0.,999)
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'border.for'
C     border.for
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'calcops.for'
C     calcops.for
C
C=======================================================================
C