C
C *** CALCULATION OF Q, V, W AND ,X ,X ,Z USING 
C     THE NEWTON'S (TRAPEZOIDAL)INTEGRAL FORMULA ***
C
      SUBROUTINE QVWXYZ(A,B,N,DNS,VP,VS,THK,ML,IER,L1,DEP1,L2,DEP2
     *                 ,ISO,CR,SR,CRX,SRX,CRY,SRY
     *                 ,QX,QY,QZ,VX,VY,WX,WY,WZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 VP(ML),VS(ML)
     *          ,QX(0:3),QY(0:3),QZ(0:3),VX(0:3),VY(0:3)
     *          ,WX(0:3),WY(0:3),WZ(0:3)
     *          ,FQX(0:3),FVX(0:3),FWX(0:3),FQZ(0:3),FWZ(0:3)
     *          ,YQX(0:3),YVX(0:3),YWX(0:3),YQZ(0:3),YWZ(0:3)
     *          ,QQX(0:3),QVX(0:3),QWX(0:3),QQZ(0:3),QWZ(0:3)
     *          ,QRX,VRX,WRX,QRZ,WRZ,QLX,VLX
     *          ,QRXZ,VRXZ,WRXZ,QRZZ,WRZZ,QLXZ,VLXZ
      DIMENSION DNS(ML),THK(ML)
      COMMON /FDAT/ OM,NL,RR
C
      IF(N.LT.1)THEN
        IER=1000
        DO 100 I=0,3
          QX(I)=(0.D0,0.D0)
          QY(I)=(0.D0,0.D0)
          QZ(I)=(0.D0,0.D0)
          VX(I)=(0.D0,0.D0)
          VY(I)=(0.D0,0.D0)
          WX(I)=(0.D0,0.D0)
          WY(I)=(0.D0,0.D0)
          WZ(I)=(0.D0,0.D0)
  100   CONTINUE
        GO TO 9999
      END IF
C
      IF(A.EQ.B)THEN
        IER=2000
        DO 150 I=0,3
          QX(I)=(0.D0,0.D0)
          QY(I)=(0.D0,0.D0)
          QZ(I)=(0.D0,0.D0)
          VX(I)=(0.D0,0.D0)
          VY(I)=(0.D0,0.D0)
          WX(I)=(0.D0,0.D0)
          WY(I)=(0.D0,0.D0)
          WZ(I)=(0.D0,0.D0)
  150   CONTINUE
        GO TO 9999
      END IF
C
      IF(ISO.EQ.1)THEN
C
        H=(B-A)/N
        CALL UVWRL1(A,DNS,VP,VS,THK,ML,L1,DEP1,L2,DEP2
     *             ,QRX,VRX,WRX,QRZ,WRZ,QLX,VLX
     *             ,QRXZ,VRXZ,WRXZ,QRZZ,WRZZ,QLXZ,VLXZ)
C
        CALL QVWXZ(A,CR,SR,YQX,YVX,YWX,YQZ,YWZ
     *             ,QRX,VRX,WRX,QRZ,WRZ,QLX,VLX
     *             ,QRXZ,VRXZ,WRXZ,QRZZ,WRZZ,QLXZ,VLXZ)
C
        DO 200 I=0,3
          QQX(I)=(0.D0,0.D0)
          QVX(I)=(0.D0,0.D0)
          QWX(I)=(0.D0,0.D0)
          QQZ(I)=(0.D0,0.D0)
          QWZ(I)=(0.D0,0.D0)
  200   CONTINUE
C
        DO 300 I=1,N
C
          X=A+H*I
          CALL UVWRL1(X,DNS,VP,VS,THK,ML,L1,DEP1,L2,DEP2
     *               ,QRX,VRX,WRX,QRZ,WRZ,QLX,VLX
     *               ,QRXZ,VRXZ,WRXZ,QRZZ,WRZZ,QLXZ,VLXZ)
C
          CALL QVWXZ(X,CR,SR,FQX,FVX,FWX,FQZ,FWZ
     *               ,QRX,VRX,WRX,QRZ,WRZ,QLX,VLX
     *               ,QRXZ,VRXZ,WRXZ,QRZZ,WRZZ,QLXZ,VLXZ)
C
          DO 310 J=0,3
            QQX(J)= QQX(J)+H*(YQX(J)+FQX(J))/2.D0
            QVX(J)= QVX(J)+H*(YVX(J)+FVX(J))/2.D0
            QWX(J)= QWX(J)+H*(YWX(J)+FWX(J))/2.D0
            QQZ(J)= QQZ(J)+H*(YQZ(J)+FQZ(J))/2.D0
            QWZ(J)= QWZ(J)+H*(YWZ(J)+FWZ(J))/2.D0
C
            YQX(J)= FQX(J)
            YVX(J)= FVX(J)
            YWX(J)= FWX(J)
            YQZ(J)= FQZ(J)
            YWZ(J)= FWZ(J)
  310    CONTINUE
C
  300   CONTINUE
C
      ELSE
C
        H=(B-A)/N
        CALL UVWRL2(A,DNS,VP,VS,THK,ML,L1,DEP1,L2,DEP2
     *             ,QRX,VRX,WRX,QRZ,WRZ,QLX,VLX
     *             ,QRXZ,VRXZ,WRXZ,QRZZ,WRZZ,QLXZ,VLXZ)
C
        CALL QVWXZ(A,CR,SR,YQX,YVX,YWX,YQZ,YWZ
     *             ,QRX,VRX,WRX,QRZ,WRZ,QLX,VLX
     *             ,QRXZ,VRXZ,WRXZ,QRZZ,WRZZ,QLXZ,VLXZ)
C
        DO 400 I=0,3
          QQX(I)=(0.D0,0.D0)
          QVX(I)=(0.D0,0.D0)
          QWX(I)=(0.D0,0.D0)
          QQZ(I)=(0.D0,0.D0)
          QWZ(I)=(0.D0,0.D0)
  400   CONTINUE
C
        DO 500 I=1,N
C
          X=A+H*I
          CALL UVWRL2(X,DNS,VP,VS,THK,ML,L1,DEP1,L2,DEP2
     *               ,QRX,VRX,WRX,QRZ,WRZ,QLX,VLX
     *               ,QRXZ,VRXZ,WRXZ,QRZZ,WRZZ,QLXZ,VLXZ)
C
          CALL QVWXZ(X,CR,SR,YQX,YVX,YWX,YQZ,YWZ
     *               ,QRX,VRX,WRX,QRZ,WRZ,QLX,VLX
     *               ,QRXZ,VRXZ,WRXZ,QRZZ,WRZZ,QLXZ,VLXZ)
C
          DO 510 J=0,3
            QQX(J)= QQX(J)+H*(YQX(J)+FQX(J))/2.D0
            QVX(J)= QVX(J)+H*(YVX(J)+FVX(J))/2.D0
            QWX(J)= QWX(J)+H*(YWX(J)+FWX(J))/2.D0
            QQZ(J)= QQZ(J)+H*(YQZ(J)+FQZ(J))/2.D0
            QWZ(J)= QWZ(J)+H*(YWZ(J)+FWZ(J))/2.D0
C
            YQX(J)= FQX(J)
            YVX(J)= FVX(J)
            YWX(J)= FWX(J)
            YQZ(J)= FQZ(J)
            YWZ(J)= FWZ(J)
  510    CONTINUE
C
  500   CONTINUE
C
      END IF
C
      DO 600 I=0,3,3
        QX(I)= QQX(I)*CR
        QY(I)= QQX(I)*SR
        QZ(I)= QQZ(I)
        VX(I)= QVX(I)*SR
        VY(I)=-QVX(I)*CR
        WX(I)= QWX(I)*CR
        WY(I)= QWX(I)*SR
        WZ(I)= QWZ(I)
  600 CONTINUE
C
      QX(1)= QQX(1)*CR +QQX(0)*CRX
      QY(1)= QQX(1)*SR +QQX(0)*SRX
      QZ(1)= QQZ(1)
      VX(1)= QVX(1)*SR +QVX(0)*SRX
      VY(1)=-QVX(1)*CR -QVX(0)*CRX
      WX(1)= QWX(1)*CR +QWX(0)*CRX
      WY(1)= QWX(1)*SR +QWX(0)*SRX
      WZ(1)= QWZ(1)
C
      QX(2)= QQX(2)*CR +QQX(0)*CRY
      QY(2)= QQX(2)*SR +QQX(0)*SRY
      QZ(2)= QQZ(2)
      VX(2)= QVX(2)*SR +QVX(0)*SRY
      VY(2)=-QVX(2)*CR -QVX(0)*CRY
      WX(2)= QWX(2)*CR +QWX(0)*CRY
      WY(2)= QWX(2)*SR +QWX(0)*SRY
      WZ(2)= QWZ(2)
C
      IER=0
C
 9999 RETURN
C
      END
C
C
C
C
      SUBROUTINE QVWXZ(WN,CR,SR,QX,VX,WX,QZ,WZ
     *                ,QRX,VRX,WRX,QRZ,WRZ,QLX,VLX
     *                ,QRXZ,VRXZ,WRXZ,QRZZ,WRZZ,QLXZ,VLXZ)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 QX(0:3),VX(0:3),WX(0:3),QZ(0:3),WZ(0:3)
     *          ,QRX,VRX,WRX,QRZ,WRZ,QLX,VLX
     *          ,QRXZ,VRXZ,WRXZ,QRZZ,WRZZ,QLXZ,VLXZ
      DIMENSION BJ0(0:2),BJ1(0:2),BJD(0:2),BJR(0:2)
C
      CALL BES01D(WN,CR,SR,BJ0,BJ1,BJD,BJR)
C
      DO 100 I=0,2
        QX(I)= QRX*BJD(I)+QLX*BJR(I)
        VX(I)= VRX*BJR(I)+VLX*BJD(I)    
        WX(I)= WRX*BJ1(I)    
        QZ(I)= QRZ*BJ1(I)
        WZ(I)= WRZ*BJ0(I)
  100 CONTINUE
C
      QX(3)= QRXZ*BJD(0)+QLXZ*BJR(0)    
      VX(3)= VRXZ*BJR(0)+VLXZ*BJD(0)    
      WX(3)= WRXZ*BJ1(0)
      QZ(3)= QRZZ*BJ1(0)
      WZ(3)= WRZZ*BJ0(0)
C
      RETURN
C
      END
