C     *******************
C     *  PROGRAM AMAIN  *
C     *******************
C
      PROGRAM AMAIN
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (ND1=202,ND3= 80,ND2=6)
      COMPLEX*16 PM1(ND1,ND1),PM2(ND1,ND1),VM1(ND1,ND1),VM2(ND1,ND1)
     *          ,VPM(ND1,ND1),V2(ND1,ND3,ND2),P2(ND1,ND3,ND2)
     *          ,VF1(ND1),VF(ND1,ND2)
      DIMENSION SX2(ND1),SZ2(ND1),WK(ND1),IWK(ND1),NCF(10)
     *         ,FX0(ND2),FZ0(ND2),DA(ND2)
      COMMON /FREQ/ FT,DT /TIME/NO,NT
     *       / XZ / FW
     *       / NS / NS2
C
      WRITE(6,6000)ND1,ND3,ND2
C
C      CALL CLOCK
C
      CALL INDATA(SX2,SZ2,FX0,FZ0,DA,NCF,ND1,ND3,ND2)
C
      DO 100 KO=1,NO
C
        OM=FT*DFLOAT(KO)
C
        CALL CALMT (SX2,SZ2,ND1,ND2,VM1,VM2,PM1,PM2,VPM,OM,VF,WK,IWK
     *             ,FX0,FZ0,DA,NCF)
        CALL CALVP2(ND1,VM2,PM1,PM2,VPM,VF,VF1,V2,P2,WK,IWK,NCF,KO,ND3)
C
  100 CONTINUE
C
      CALL OUTDAT(V2,P2,ND1,ND3)
C
C      CALL CLOCK(TIME)
C      WRITE(6,6010)TIME
C
      STOP
C
 6000 FORMAT(///10X,'DIMENSION DATA'//
     *    10X,'ND1( 2 * Elem.Num. )= ',I5/
     *   ,10X,'ND3( Freq.Num. )= ',I5,'  ND2( Fault Num. )= ',I5/)
 6010 FORMAT(///10X,'CPU TIME =',F15.7/)
      END
C     ***********************
C     *  SUBROUTINE INDATA  *
C     ***********************
C
      SUBROUTINE INDATA(SX2,SZ2,FX0,FZ0,DA,NCF,ND1,ND3,ND2)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION SX2(*),SZ2(*),FX0(*),FZ0(*),DA(*),NCF(10)
      COMMON /FREQ/ FT,DT /TIME/ NO,NT
     *       /CONS/ RT,VR,DV,VS1,VS2,DN1,DN2
     *       / XZ / FW  / NFA / NFAI
     *       / NS / NS2  /QVAL/Q1,Q2 /NQOUT/NQO
     *       / GAUSS / NG / RADPATN / NRAD
      CHARACTER FILE1*10
C
C     *****  READ DATA  *****
C
      WRITE(5,*)' ENTER INPUT FILE ;'
      READ(5,4999)FILE1
 4999 FORMAT(A10)
      OPEN(4,FILE=FILE1,STATUS='OLD')
      READ(4,5000) (NCF(I),I=1,9)
      READ(4,5010) FT,NO,DT,NT,NG,NQO,NRAD
      READ(4,5020) VS1,DN1,Q1
      READ(4,5020) VS2,DN2,Q2
      READ(4,5030) RT,VR,DV,FW
      READ(4,5040) NFAI
      DO 100 N=1,NFAI
      READ(4,5050) FX0(N),FZ0(N),DA(N)
      DA(N)=DA(N)*DATAN(1.D0)/45.D0
  100 CONTINUE
      READ(4,5060) NS2
      READ(4,5070) (SX2(I),SZ2(I),SX2(I+1),SZ2(I+1)
     *    ,SX2(I+2),SZ2(I+2),I=1,NS2+1,3)
      CLOSE(4)
C
      IF(ND1.LE.NS2)WRITE(6,*)' ND1 < NS2 !!!!!!!!!!!!!!!!!!!!!!!!!!'
      IF(ND3.LE.NO )WRITE(6,*)' ND3 < NO  !!!!!!!!!!!!!!!!!!!!!!!!!!'
      IF(ND2.LE.NFAI)WRITE(6,*)' ND2 < NFAI !!!!!!!!!!!!!!!!!!!!!!!!!!'
C
C     *****  WRITE DATA  *****
C
      WRITE(6,6000) FT,NO,DT,NT,NG
      WRITE(6,6010) RT,VR,DV,VS1,VS2,DN1,DN2,FW
      WRITE(6,6020) NFAI
      DO 200 N=1,NFAI
      WRITE(6,6030) FX0(N),FZ0(N),DA(N)
  200 CONTINUE
      WRITE(6,6040) NQO,Q1,Q2,NRAD,NS2
      WRITE(6,6050) (SX2(I),SZ2(I),SX2(I+1),SZ2(I+1)
     *   ,SX2(I+2),SZ2(I+2),I=1,NS2+1,3)
C
      RETURN
C
 5000 FORMAT(9I5)
 5010 FORMAT(F10.0,I5,F10.0,4I5)
 5020 FORMAT(3F10.0)
 5030 FORMAT(4F10.0)
 5040 FORMAT(I10)
 5050 FORMAT(3F10.0)
 5060 FORMAT(////I15)
 5070 FORMAT(6F10.0)
 6000 FORMAT(///2X,'FT                 =',F15.7
     *         ,2X,'NUMBER OF FREQUENCY=',I8/
     *         ,2X,'DT                 =',F15.7
     *         ,2X,'NUMBER OF TIME     =',I8/
     *         ,2X,'NUMBER OF GAUSS IN.=',I8)
 6010 FORMAT(  /2X,'RISE TIME          =',F15.7
     *         ,2X,'RUPTURE VELOCITY   =',F15.7/
     *         ,2X,'DISPLACEMENT OF F. =',F15.7
     *         ,2X,'VS      (1)        =',F15.7/
     *         ,2X,'VS      (2)        =',F15.7
     *         ,2X,'DENSITY (1)        =',F15.7/
     *         ,2X,'DENSITY (2)        =',F15.7
     *         ,2X,'FAULT WIDTH        =',F15.7)
 6020 FORMAT(  /2X,'NUMBER OF FAULT    =',I8/)
 6030 FORMAT(   2X,'FXO,FZO,DA         =',3F15.7)
 6040 FORMAT(  /2X,' INPUT OF Q VALURE ( 1: IN , 0: OUT) =',I5,
     *         /2X,' Q VALUE ; Q1, Q2 =',2F15.7
     *         /2X,' RADIATION PATERN  ( 1: OUT, 0: IN ) =',I5,
     *         /2X,' NUMBER OF BOUNDARY2 ELEM.',I4
     *       ///2X,' (SX2,SZ2)*3 -------------'/)
 6050 FORMAT(2X,2F10.2,2X,2F10.2,2X,2F10.2)
C
      END
C     ***********************
C     *  SUBROUTINE CALMT   *
C     ***********************
C
      SUBROUTINE CALMT (SX2,SZ2,ND1,ND2,VM1,VM2,PM1,PM2,VPM,OM,
     *                 VF,WK,IWK,FX0,FZ0,DA,NCF)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 PM1(ND1,ND1),PM2(ND1,ND1),VM1(ND1,ND1),VM2(ND1,ND1)
     *          ,VPM(ND1,ND1),VF(ND1,*),RDET
      DIMENSION SX2(ND1),SZ2(ND1),NCF(10),WK(ND1),IWK(ND1)
     *         ,FX0(*),FZ0(*),DA(*)
C
      COMMON / NS / NS2 / NFA / NFAI
C
C     *****  INITIALIZATION  *****
C
      DO 150 I=1,ND1
       DO 110 N=1,ND2
       VF(I,N)=(0.D0,0.D0)
  110 CONTINUE
        DO 100 J=1,ND1
         VM1(I,J)=(0.D0,0.D0)
         VM2(I,J)=(0.D0,0.D0)
         PM1(I,J)=(0.D0,0.D0)
         PM2(I,J)=(0.D0,0.D0)
         VPM(I,J)=(0.D0,0.D0)
  100   CONTINUE
  150 CONTINUE
C
C     *****  MAKE MATRIX TM  *****
C
      CALL PV22 (SX2,SZ2,ND1,PM1,VM1,OM,NS2,1,NCF)
      CALL PV22 (SX2,SZ2,ND1,PM2,VM2,OM,NS2,2,NCF)
C
      CALL ZIF1M(VM1,NS2,ND1,1.D-10,RDET,IDET,WK,IWK,IER)
      WRITE(6,1000)OM,RDET,IDET
 1000 FORMAT(' OM,RDET(RE,IM),IDET)=',3F8.4,I5)
C
      DO 290 I=1,NS2
        DO 280 J=1,NS2
        RDET=(0.D0,0.D0)
          DO 285 K=1,NS2
          RDET=RDET+VM1(I,K)*PM1(K,J)
  285     CONTINUE
        VPM(I,J)=RDET
  280   CONTINUE
  290 CONTINUE
C
C
C     *****  MAKE VECTOR VF *****
C
      CALL MKVCF (SX2,SZ2,OM,VF,NS2,FX0,FZ0,DA,ND1,NCF)
C
C     *****  CHECKER'S FLAG NO.3  *****
C
      IF(NCF(3).EQ.1) THEN
      WRITE(6,*) '*****  FROM CALMAT (NCF=3) MATRIX VM1 *****'
      WRITE(6,*) ((VM1(I,J),J=1,NS2),I=1,NS2)
      WRITE(6,*) '*****  MATRIX VM2 *****'
      WRITE(6,*) ((VM2(I,J),J=1,NS2),I=1,NS2)
      WRITE(6,*) '*****  MATRIX PM1 *****'
      WRITE(6,*) ((PM1(I,J),J=1,NS2),I=1,NS2)
      WRITE(6,*) '*****  MATRIX PM2 *****'
      WRITE(6,*) ((PM2(I,J),J=1,NS2),I=1,NS2)
      WRITE(6,*) '*****  MATRIX VPM *****'
      WRITE(6,*) ((VPM(I,J),J=1,NS2),I=1,NS2)
      WRITE(6,*) '*****  VECTOR VF  *****'
      WRITE(6,*) ((VF(I,N),I=1,NS2),N=1,NFAI)
      END IF
C
      RETURN
      END
C     ***********************
C     *  SUBROUTINE CALVP2  *
C     ***********************
C
      SUBROUTINE CALVP2(ND1,VM2,PM1,PM2,VPM,VF,VF1,V2,P2,
     *                  WK,IWK,NCF,KO,ND3)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 PM1(ND1,ND1),PM2(ND1,ND1),VM2(ND1,ND1),DD
     *    ,VPM(ND1,ND1),VF(ND1,*),VF1(*),V2(ND1,ND3,*),P2(ND1,ND3,*)
      DIMENSION NCF(10),WK(ND1),IWK(ND1)
C
      COMMON / NS / NS2 / NFA / NFAI
C
      ND3=ND3
      DO 110 I=1,NS2
      DO 100 J=1,NS2
      PM1(I,J)=(0.D0,0.D0)
  100 CONTINUE
  110 CONTINUE
C
      DO 220 I=1,NS2
        DO 210 J=1,NS2
        DD=(0.D0,0.D0)
          DO 200 K=1,NS2
          DD=DD+VM2(I,K)*VPM(K,J)
  200     CONTINUE
        PM1(I,J)=PM2(I,J)+DD
  210   CONTINUE
  220 CONTINUE
C
      DO 330 N=1,NFAI
C
       DO 300 K=1,NS2
       VF1(K)=VF(K,N)
        DO 300 L=1,NS2
        PM2(K,L)=PM1(K,L)
  300  CONTINUE
C
      CALL ZLF1M(PM2,NS2,ND1,VF1,1.D-10,1,IWK,WK,IER)
C
       DO 320 I=1,NS2
        P2(I,KO,N)=(0.D0,0.D0)
        V2(I,KO,N)=VF1(I)
        DO 310 J=1,NS2
          P2(I,KO,N)=P2(I,KO,N)+VPM(I,J)*VF1(J)
  310   CONTINUE
  320  CONTINUE
C
  330 CONTINUE
C
C     *****  CHECKER'S FLAG NO.4  *****
C
      IF(NCF(4).EQ.1) THEN
      WRITE(6,*) '*****  FROM CALVP2 (NCF=4) MATRIX V2  *****'
      WRITE(6,*)((V2(I,KO,N),I=1,NS2),N=1,NFAI)
      WRITE(6,*) '*****  MATRIX P2  *****'
      WRITE(6,*) ((P2(I,KO,N),I=1,NS2),N=1,NFAI)
      END IF
C
      RETURN
      END
      SUBROUTINE GAUSPV(A,B,N,PSPR,PSPI,VSPR,VSPI,IER)
C
C   S-1511-1 COPYRIGHT HITACHI,LTD. 1980
C
C   STATUS              - S-1511-1 01-00
C
C   HISTORY             - DATE.    1980.3
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /EXTE/ SS,XN,ZN,AA1,AA2,AA3,AA4
      DIMENSION Y(6,12),X(6,12),AA(24),U(24)
     1,A1(19),A2(19),A3(19),A4(15)
     2,X1(19),X2(19),X3(19),X4(15)
C      DIMENSION NAME$(2)
      EQUIVALENCE (Y(1),A1(1)),(Y(20),A2(1)),(Y(39),A3(1)),(Y(58),A4(1))
     1,           (X(1),X1(1)),(X(20),X2(1)),(X(39),X3(1)),(X(58),X4(1))
      DATA  X1    /               0.86113631159405260D+00
     1,  0.96028985649753620D+00, 0.98156063424671920D+00
     2,  0.98940093499164990D+00, 0.99312859918509490D+00
     3,  0.99518721999702140D+00, 0.33998104358485630D+00
     4,  0.79666647741362670D+00, 0.90411725637047490D+00
     5,  0.94457502307323260D+00, 0.96397192727791380D+00
     6,  0.97472855597130950D+00, 0.00000000000000000D-40
     7,  0.52553240991632900D+00, 0.76990267419430470D+00
     8,  0.86563120238783170D+00, 0.91223442825132590D+00
     9,  0.93827455200273280D+00, 0.00000000000000000D-40 /
      DATA  X2     /              0.18343464249564980D+00
     1,  0.58731795428661740D+00, 0.75540440835500300D+00
     2,  0.83911697182221880D+00, 0.88641552700440100D+00
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.36783149899818020D+00, 0.61787624440264370D+00
     5,  0.74633190646015080D+00, 0.82000198597390290D+00
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.12523340851146890D+00, 0.45801677765722740D+00
     8,  0.63605368072651500D+00, 0.74012419157855440D+00
     9,  0.00000000000000000D-40, 0.00000000000000000D-40 /
      DATA  X3     /              0.00000000000000000D-40
     1,  0.28160355077925890D+00, 0.51086700195082710D+00
     2,  0.64809365193697560D+00, 0.00000000000000000D-40
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.95012509837637440D-01, 0.37370608871541960D+00
     5,  0.54542147138883950D+00, 0.00000000000000000D-40
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.00000000000000000D-40, 0.22778585114164510D+00
     8,  0.43379350762604510D+00, 0.00000000000000000D-40
     9,  0.00000000000000000D-40, 0.00000000000000000D-40 /
      DATA  X4     /              0.51086700195082710D+00
     1,  0.76526521133497330D-01, 0.31504267969616340D+00
     2,  0.00000000000000000D-40, 0.00000000000000000D-40
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.00000000000000000D-40, 0.19111886747361630D+00
     5,  0.00000000000000000D-40, 0.00000000000000000D-40
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.00000000000000000D-40, 0.64056892862605620D-01 /
      DATA  A1    /               0.34785484513745390D+00
     1,  0.10122853629037630D+00, 0.47175336386511830D-01
     2,  0.27152459411754100D-01, 0.17614007139152120D-01
     3,  0.12341229799987200D-01, 0.65214515486254620D+00
     4,  0.22238103445337450D+00, 0.10693932599531840D+00
     5,  0.62253523938647890D-01, 0.40601429800386940D-01
     6,  0.28531388628933660D-01, 0.00000000000000000D-40
     7,  0.31370664587788730D+00, 0.16007832854334620D+00
     8,  0.95158511682492780D-01, 0.62672048334109070D-01
     9,  0.44277438817419810D-01, 0.00000000000000000D-40 /
      DATA  A2     /              0.36268378337836200D+00
     1,  0.20316742672306590D+00, 0.12462897125553390D+00
     2,  0.83276741576704750D-01, 0.59298584915436780D-01
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.23349253653835480D+00, 0.14959598881657670D+00
     5,  0.10193011981724040D+00, 0.73346481411080310D-01
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.24914704581340280D+00, 0.16915651939500250D+00
     8,  0.11819453196151840D+00, 0.86190161531953270D-01
     9,  0.00000000000000000D-40, 0.00000000000000000D-40 /
      DATA  A3     /              0.00000000000000000D-40
     1,  0.18260341504492360D+00, 0.13168863844917660D+00
     2,  0.97618652104113890D-01, 0.00000000000000000D-40
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.18945061045506850D+00, 0.14209610931838210D+00
     5,  0.10744427011596560D+00, 0.00000000000000000D-40
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.00000000000000000D-40, 0.14917298647260380D+00
     8,  0.11550566805372560D+00, 0.00000000000000000D-40
     9,  0.00000000000000000D-40, 0.00000000000000000D-40 /
      DATA  A4     /              0.13168863844917660D+00
     1,  0.15275338713072580D+00, 0.12167047292780340D+00
     2,  0.00000000000000000D-40, 0.00000000000000000D-40
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.00000000000000000D-40, 0.12583745634682830D+00
     5,  0.00000000000000000D-40, 0.00000000000000000D-40
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.00000000000000000D-40, 0.12793819534675220D+00 /
C      DATA NAME$ /4H$DQB,4HGM  /
      IER=0
      K=N
      IF(N-1) 600,100,100
  100 IF(N-6) 200,200,600
  200 NN=4*N
      C1=(B+A)/2.D0
      C2=(B-A)/2.D0
      N2=N+N
      DO 1 J=1,N2
      U(J)=C1-C2*X(N,J)
      U(NN-J+1)  = C1+C2*X(N,J)
      AA(J)=Y(N,J)
    1 AA(NN-J+1) = Y(N,J)
      PSPR=0.D0
      PSPI=0.D0
      VSPR=0.D0
      VSPI=0.D0
      DO 2 J=1,NN
      SA=U(J)
      RXRS=AA1*SA+AA2
      RZRS=AA3*SA+AA4
      RS=DSQRT(RXRS*RXRS+RZRS*RZRS)
      PXZN=RXRS/RS*XN+RZRS/RS*ZN
      SR=SS*RS
      CALL BESJN0(SR,BJ0,BY0)
      CALL BESJN1(SR,BJ1,BY1)
      PSPR=PSPR+AA(J)*BY1*PXZN
      PSPI=PSPI+AA(J)*BJ1*PXZN
      VSPR=VSPR+AA(J)*BY0
      VSPI=VSPI+AA(J)*BJ0
    2 CONTINUE
      PSPR=C2*PSPR
      PSPI=C2*PSPI
      VSPR=C2*VSPR
      VSPI=C2*VSPI
      N=K
C      CALL $SURTM(NAME$,IER)
      RETURN
  600 CONTINUE
C     CALL $SUERM(NAME$,1200,IER)
      N=6
      GO TO 200
      END
      SUBROUTINE GASPVM(A,B,N,PSPR,PSPI,VSPR,VSPI,IER)
C
C   S-1511-1 COPYRIGHT HITACHI,LTD. 1980
C
C   STATUS              - S-1511-1 01-00
C
C   HISTORY             - DATE.    1980.3
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /EXTE/ SS,XN,ZN,AA1,AA2,AA3,AA4
     *       /EXTEM/ AA1M,AA2M,AA3M,AA4M
      DIMENSION Y(6,12),X(6,12),AA(24),U(24)
     1,A1(19),A2(19),A3(19),A4(15)
     2,X1(19),X2(19),X3(19),X4(15)
C      DIMENSION NAME$(2)
      EQUIVALENCE (Y(1),A1(1)),(Y(20),A2(1)),(Y(39),A3(1)),(Y(58),A4(1))
     1,           (X(1),X1(1)),(X(20),X2(1)),(X(39),X3(1)),(X(58),X4(1))
      DATA  X1    /               0.86113631159405260D+00
     1,  0.96028985649753620D+00, 0.98156063424671920D+00
     2,  0.98940093499164990D+00, 0.99312859918509490D+00
     3,  0.99518721999702140D+00, 0.33998104358485630D+00
     4,  0.79666647741362670D+00, 0.90411725637047490D+00
     5,  0.94457502307323260D+00, 0.96397192727791380D+00
     6,  0.97472855597130950D+00, 0.00000000000000000D-40
     7,  0.52553240991632900D+00, 0.76990267419430470D+00
     8,  0.86563120238783170D+00, 0.91223442825132590D+00
     9,  0.93827455200273280D+00, 0.00000000000000000D-40 /
      DATA  X2     /              0.18343464249564980D+00
     1,  0.58731795428661740D+00, 0.75540440835500300D+00
     2,  0.83911697182221880D+00, 0.88641552700440100D+00
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.36783149899818020D+00, 0.61787624440264370D+00
     5,  0.74633190646015080D+00, 0.82000198597390290D+00
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.12523340851146890D+00, 0.45801677765722740D+00
     8,  0.63605368072651500D+00, 0.74012419157855440D+00
     9,  0.00000000000000000D-40, 0.00000000000000000D-40 /
      DATA  X3     /              0.00000000000000000D-40
     1,  0.28160355077925890D+00, 0.51086700195082710D+00
     2,  0.64809365193697560D+00, 0.00000000000000000D-40
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.95012509837637440D-01, 0.37370608871541960D+00
     5,  0.54542147138883950D+00, 0.00000000000000000D-40
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.00000000000000000D-40, 0.22778585114164510D+00
     8,  0.43379350762604510D+00, 0.00000000000000000D-40
     9,  0.00000000000000000D-40, 0.00000000000000000D-40 /
      DATA  X4     /              0.51086700195082710D+00
     1,  0.76526521133497330D-01, 0.31504267969616340D+00
     2,  0.00000000000000000D-40, 0.00000000000000000D-40
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.00000000000000000D-40, 0.19111886747361630D+00
     5,  0.00000000000000000D-40, 0.00000000000000000D-40
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.00000000000000000D-40, 0.64056892862605620D-01 /
      DATA  A1    /               0.34785484513745390D+00
     1,  0.10122853629037630D+00, 0.47175336386511830D-01
     2,  0.27152459411754100D-01, 0.17614007139152120D-01
     3,  0.12341229799987200D-01, 0.65214515486254620D+00
     4,  0.22238103445337450D+00, 0.10693932599531840D+00
     5,  0.62253523938647890D-01, 0.40601429800386940D-01
     6,  0.28531388628933660D-01, 0.00000000000000000D-40
     7,  0.31370664587788730D+00, 0.16007832854334620D+00
     8,  0.95158511682492780D-01, 0.62672048334109070D-01
     9,  0.44277438817419810D-01, 0.00000000000000000D-40 /
      DATA  A2     /              0.36268378337836200D+00
     1,  0.20316742672306590D+00, 0.12462897125553390D+00
     2,  0.83276741576704750D-01, 0.59298584915436780D-01
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.23349253653835480D+00, 0.14959598881657670D+00
     5,  0.10193011981724040D+00, 0.73346481411080310D-01
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.24914704581340280D+00, 0.16915651939500250D+00
     8,  0.11819453196151840D+00, 0.86190161531953270D-01
     9,  0.00000000000000000D-40, 0.00000000000000000D-40 /
      DATA  A3     /              0.00000000000000000D-40
     1,  0.18260341504492360D+00, 0.13168863844917660D+00
     2,  0.97618652104113890D-01, 0.00000000000000000D-40
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.18945061045506850D+00, 0.14209610931838210D+00
     5,  0.10744427011596560D+00, 0.00000000000000000D-40
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.00000000000000000D-40, 0.14917298647260380D+00
     8,  0.11550566805372560D+00, 0.00000000000000000D-40
     9,  0.00000000000000000D-40, 0.00000000000000000D-40 /
      DATA  A4     /              0.13168863844917660D+00
     1,  0.15275338713072580D+00, 0.12167047292780340D+00
     2,  0.00000000000000000D-40, 0.00000000000000000D-40
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.00000000000000000D-40, 0.12583745634682830D+00
     5,  0.00000000000000000D-40, 0.00000000000000000D-40
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.00000000000000000D-40, 0.12793819534675220D+00 /
C      DATA NAME$ /4H$DQB,4HGM  /
      IER=0
      K=N
      IF(N-1) 600,100,100
  100 IF(N-6) 200,200,600
  200 NN=4*N
      C1=(B+A)/2.D0
      C2=(B-A)/2.D0
      N2=N+N
      DO 1 J=1,N2
      U(J)=C1-C2*X(N,J)
      U(NN-J+1)  = C1+C2*X(N,J)
      AA(J)=Y(N,J)
    1 AA(NN-J+1) = Y(N,J)
      PSPR=0.D0
      PSPI=0.D0
      VSPR=0.D0
      VSPI=0.D0
      DO 2 J=1,NN
      SA=U(J)
      RXRS=AA1*SA+AA2
      RZRS=AA3*SA+AA4
      RS=DSQRT(RXRS*RXRS+RZRS*RZRS)
      PXZN=RXRS/RS*XN+RZRS/RS*ZN
      RXRSM=AA1M*SA+AA2M
      RZRSM=AA3M*SA+AA4M
      RSM=DSQRT(RXRSM*RXRSM+RZRSM*RZRSM)
      PXZNM=RXRSM/RSM*XN+RZRSM/RSM*ZN
      SR =SS*RS
      SRM=SS*RSM
      CALL BESJN0(SR,BJ0,BY0)
      CALL BESJN1(SR,BJ1,BY1)
      CALL BESJN0(SRM,BJ0M,BY0M)
      CALL BESJN1(SRM,BJ1M,BY1M)
      PSPR=PSPR+AA(J)*(BY1*PXZN+BY1M*PXZNM)
      PSPI=PSPI+AA(J)*(BJ1*PXZN+BJ1M*PXZNM)
      VSPR=VSPR+AA(J)*(BY0+BY0M)
      VSPI=VSPI+AA(J)*(BJ0+BJ0M)
    2 CONTINUE
      PSPR=C2*PSPR
      PSPI=C2*PSPI
      VSPR=C2*VSPR
      VSPI=C2*VSPI
      N=K
C      CALL $SURTM(NAME$,IER)
      RETURN
  600 CONTINUE
C     CALL $SUERM(NAME$,1200,IER)
      N=6
      GO TO 200
      END
      SUBROUTINE GASPVQ(A,B,N,PR,PI,VR,VI,PRM,PIM,VRM,VIM,IER)
C
C   S-1511-1 COPYRIGHT HITACHI,LTD. 1980
C
C   STATUS              - S-1511-1 01-00
C
C   HISTORY             - DATE.    1980.3
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /EXTE/ SS,XN,ZN,AA1,AA2,AA3,AA4
     *       /EXTEM/ AA1M,AA2M,AA3M,AA4M
      DIMENSION Y(6,12),X(6,12),AA(24),U(24)
     1,A1(19),A2(19),A3(19),A4(15)
     2,X1(19),X2(19),X3(19),X4(15)
C      DIMENSION NAME$(2)
      EQUIVALENCE (Y(1),A1(1)),(Y(20),A2(1)),(Y(39),A3(1)),(Y(58),A4(1))
     1,           (X(1),X1(1)),(X(20),X2(1)),(X(39),X3(1)),(X(58),X4(1))
      DATA  X1    /               0.86113631159405260D+00
     1,  0.96028985649753620D+00, 0.98156063424671920D+00
     2,  0.98940093499164990D+00, 0.99312859918509490D+00
     3,  0.99518721999702140D+00, 0.33998104358485630D+00
     4,  0.79666647741362670D+00, 0.90411725637047490D+00
     5,  0.94457502307323260D+00, 0.96397192727791380D+00
     6,  0.97472855597130950D+00, 0.00000000000000000D-40
     7,  0.52553240991632900D+00, 0.76990267419430470D+00
     8,  0.86563120238783170D+00, 0.91223442825132590D+00
     9,  0.93827455200273280D+00, 0.00000000000000000D-40 /
      DATA  X2     /              0.18343464249564980D+00
     1,  0.58731795428661740D+00, 0.75540440835500300D+00
     2,  0.83911697182221880D+00, 0.88641552700440100D+00
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.36783149899818020D+00, 0.61787624440264370D+00
     5,  0.74633190646015080D+00, 0.82000198597390290D+00
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.12523340851146890D+00, 0.45801677765722740D+00
     8,  0.63605368072651500D+00, 0.74012419157855440D+00
     9,  0.00000000000000000D-40, 0.00000000000000000D-40 /
      DATA  X3     /              0.00000000000000000D-40
     1,  0.28160355077925890D+00, 0.51086700195082710D+00
     2,  0.64809365193697560D+00, 0.00000000000000000D-40
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.95012509837637440D-01, 0.37370608871541960D+00
     5,  0.54542147138883950D+00, 0.00000000000000000D-40
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.00000000000000000D-40, 0.22778585114164510D+00
     8,  0.43379350762604510D+00, 0.00000000000000000D-40
     9,  0.00000000000000000D-40, 0.00000000000000000D-40 /
      DATA  X4     /              0.51086700195082710D+00
     1,  0.76526521133497330D-01, 0.31504267969616340D+00
     2,  0.00000000000000000D-40, 0.00000000000000000D-40
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.00000000000000000D-40, 0.19111886747361630D+00
     5,  0.00000000000000000D-40, 0.00000000000000000D-40
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.00000000000000000D-40, 0.64056892862605620D-01 /
      DATA  A1    /               0.34785484513745390D+00
     1,  0.10122853629037630D+00, 0.47175336386511830D-01
     2,  0.27152459411754100D-01, 0.17614007139152120D-01
     3,  0.12341229799987200D-01, 0.65214515486254620D+00
     4,  0.22238103445337450D+00, 0.10693932599531840D+00
     5,  0.62253523938647890D-01, 0.40601429800386940D-01
     6,  0.28531388628933660D-01, 0.00000000000000000D-40
     7,  0.31370664587788730D+00, 0.16007832854334620D+00
     8,  0.95158511682492780D-01, 0.62672048334109070D-01
     9,  0.44277438817419810D-01, 0.00000000000000000D-40 /
      DATA  A2     /              0.36268378337836200D+00
     1,  0.20316742672306590D+00, 0.12462897125553390D+00
     2,  0.83276741576704750D-01, 0.59298584915436780D-01
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.23349253653835480D+00, 0.14959598881657670D+00
     5,  0.10193011981724040D+00, 0.73346481411080310D-01
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.24914704581340280D+00, 0.16915651939500250D+00
     8,  0.11819453196151840D+00, 0.86190161531953270D-01
     9,  0.00000000000000000D-40, 0.00000000000000000D-40 /
      DATA  A3     /              0.00000000000000000D-40
     1,  0.18260341504492360D+00, 0.13168863844917660D+00
     2,  0.97618652104113890D-01, 0.00000000000000000D-40
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.18945061045506850D+00, 0.14209610931838210D+00
     5,  0.10744427011596560D+00, 0.00000000000000000D-40
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.00000000000000000D-40, 0.14917298647260380D+00
     8,  0.11550566805372560D+00, 0.00000000000000000D-40
     9,  0.00000000000000000D-40, 0.00000000000000000D-40 /
      DATA  A4     /              0.13168863844917660D+00
     1,  0.15275338713072580D+00, 0.12167047292780340D+00
     2,  0.00000000000000000D-40, 0.00000000000000000D-40
     3,  0.00000000000000000D-40, 0.00000000000000000D-40
     4,  0.00000000000000000D-40, 0.12583745634682830D+00
     5,  0.00000000000000000D-40, 0.00000000000000000D-40
     6,  0.00000000000000000D-40, 0.00000000000000000D-40
     7,  0.00000000000000000D-40, 0.12793819534675220D+00 /
C      DATA NAME$ /4H$DQB,4HGM  /
      IER=0
      K=N
      IF(N-1) 600,100,100
  100 IF(N-6) 200,200,600
  200 NN=4*N
      C1=(B+A)/2.D0
      C2=(B-A)/2.D0
      N2=N+N
      DO 1 J=1,N2
      U(J)=C1-C2*X(N,J)
      U(NN-J+1)  = C1+C2*X(N,J)
      AA(J)=Y(N,J)
    1 AA(NN-J+1) = Y(N,J)
      PR=0.D0
      PRM=0.D0
      PI=0.D0
      PIM=0.D0
      VR=0.D0
      VRM=0.D0
      VI=0.D0
      VIM=0.D0
      DO 2 J=1,NN
      SA=U(J)
      RXRS=AA1*SA+AA2
      RZRS=AA3*SA+AA4
      RS=DSQRT(RXRS*RXRS+RZRS*RZRS)
      PXZN=RXRS/RS*XN+RZRS/RS*ZN
      RXRSM=AA1M*SA+AA2M
      RZRSM=AA3M*SA+AA4M
      RSM=DSQRT(RXRSM*RXRSM+RZRSM*RZRSM)
      PXZNM=RXRSM/RSM*XN+RZRSM/RSM*ZN
      SR =SS*RS
      SRM=SS*RSM
      CALL BESJN0(SR,BJ0,BY0)
      CALL BESJN1(SR,BJ1,BY1)
      CALL BESJN0(SRM,BJ0M,BY0M)
      CALL BESJN1(SRM,BJ1M,BY1M)
      PR=PR+AA(J)*BY1*PXZN
      PRM=PRM+AA(J)*BY1M*PXZNM
      PI=PI+AA(J)*BJ1*PXZN
      PIM=PIM+AA(J)*BJ1M*PXZNM
      VR=VR+AA(J)*BY0
      VRM=VRM+AA(J)*BY0M
      VI=VI+AA(J)*BJ0
      VIM=VIM+AA(J)*BJ0M
    2 CONTINUE
      PR=C2*PR
      PRM=C2*PRM
      PI=C2*PI
      PIM=C2*PIM
      VR=C2*VR
      VRM=C2*VRM
      VI=C2*VI
      VIM=C2*VIM
      N=K
C      CALL $SURTM(NAME$,IER)
      RETURN
  600 CONTINUE
C      CALL $SUERM(NAME$,1200,IER)
      N=6
      GO TO 200
      END
C     ***********************
C     *  SUBROUTINE MKVCF   *
C     ***********************
C
      SUBROUTINE MKVCF (SX,SZ,OM,F,NI,FX0,FZ0,DA,ND1,NCF)
C
      IMPLICIT REAL*8(A-H,M,O-Z)
      COMPLEX*16 F(ND1,*),M0W,H1,H2
      DIMENSION SX(*),SZ(*),NCF(10),FX0(*),FZ0(*),DA(*)
      COMMON /CONS/ RT,VR,DV,VS1,VS2,DN1,DN2 / RADPATN / NRAD
     *       / XZ / FW / QVAL / Q1,Q2 / NQOUT / NQO / NFA/ NFAI
C
      ND1=ND1
      VS=VS2
      DN=DN2
      QV=Q2
C
      SS=OM/VS
C
      TW=FW/VR
      XW=OM*TW/2.D0
      XT=OM*RT/2.D0
      M0=DN*VS*VS*DV*FW
      M0W=M0*DSIN(XT)/XT*DSIN(XW)/XW
     *    *CDEXP(-DCMPLX(0.D0,1.D0)*(XT+XW))
C
      DO 110 N=1,NFAI
C
      XF=FX0(N)
      ZF=FZ0(N)
      XN=+DSIN(DA(N))
      ZN=-DCOS(DA(N))
C
      DO 100 I=1,NI
C
        XI=(SX(I)+SX(I+1))/2.D0
        ZI=(SZ(I)+SZ(I+1))/2.D0
C
        RX=XF-XI
        RZ=ZF-ZI
        RF1=DSQRT(RX*RX+RZ*RZ)
C
        CALL BESJN1(SS*RF1,Y2,Y1)
        R1=(RX*XN+RZ*ZN)/RF1
        IF(NRAD.NE.0)R1=-1.D0
        H1=DCMPLX(Y2,-Y1)*R1
        IF(NQO.NE.0)THEN
        H1=H1*DEXP(-OM*RF1/(2.D0*QV*VS))
        END IF
C
        RZ=ZF+ZI
        RF2=DSQRT(RX*RX+RZ*RZ)
C
        CALL BESJN1(SS*RF2,Y2,Y1)
        R2=(RX*XN+RZ*ZN)/RF2
        IF(NRAD.NE.0)R2=-1.D0
        H2=DCMPLX(Y2,-Y1)*R2
        IF(NQO.NE.0)THEN
        H2=H2*DEXP(-OM*RF2/(2.D0*QV*VS))
        END IF
C
        F(I,N)=(H2+H1)*.25D0/DN/VS**3*M0W
C
  100 CONTINUE
  110 CONTINUE
C
C     *****  CHECKER'S FLAG NO.2  *****
C
      IF(NCF(2).EQ.1) THEN
      WRITE(6,*) ' *****  FROM MKVCF  (NCF=2) VECTOR VF  *****'
      WRITE(6,*) ((F(I,N),I=1,NI),N=1,NFAI)
      END IF
C
      RETURN
      END
C     ***********************
C     *  SUBROUTINE OUTDAT  *
C     ***********************
C
      SUBROUTINE OUTDAT(V2,P2,ND1,ND3)
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 V2(ND1,ND3,*),P2(ND1,ND3,*)
      COMMON /TIME/NO,NT
     *       / NS / NS2 / NFA/ NFAI
      CHARACTER FILE1*10
C
      WRITE(5,*)' ENTER OUTPUT FILE;'
      READ(5,999)FILE1
  999 FORMAT(A10)
      OPEN(7,FILE=FILE1,STATUS='NEW')
      DO 120 N=1,NFAI
      DO 110 KO=1,NO
      DO 100 I=1,NS2,3
        WRITE(7,1000)V2(I,KO,N),V2(I+1,KO,N),V2(I+2,KO,N)
        WRITE(7,1000)P2(I,KO,N),P2(I+1,KO,N),P2(I+2,KO,N)
  100  CONTINUE
  110 CONTINUE
  120 CONTINUE
      CLOSE(7)
C
      RETURN
C
 1000 FORMAT(6D13.5)
C
      END
C     ***********************
C     *  SUBROUTINE PV22    *
C     ***********************
C
      SUBROUTINE PV22 (SX2,SZ2,ND1,P,V,OM,NI,IND,NCF)
C
      IMPLICIT REAL*8(A-H,L,O-Z)
      COMPLEX*16 P(ND1,ND1),V(ND1,ND1)
      DIMENSION SX2(ND1),SZ2(ND1),NCF(10)
      COMMON /CONS/ RT,VRR,DV,VS1,VS2,DN1,DN2
     *       /EXTE/ SS,XN,ZN,A1,A2,A3,A4
     *       /EXTEM/ A1M,A2M,A3M,A4M
     *       /GAUSS/NG /QVAL/Q1,Q2 /NQOUT/NQO
C
      ND1=ND1
      IF(IND.EQ.1) THEN
        VS=VS1
        DN=DN1
        QV=Q1
      ELSE
        VS=VS2
        DN=DN2
        QV=Q2
      END IF
C
      SS=OM/VS
C
      PAI=3.1415926535D0
      Q  =0.5772156649D0
C
      DO 200 I=1,NI
C
        XI=(SX2(I)+SX2(I+1))/2.D0
        ZI=(SZ2(I)+SZ2(I+1))/2.D0
C
        DO 100 J=1,NI
C
          XJ1=SX2(J)
          ZJ1=SZ2(J)
          XJ2=SX2(J+1)
          ZJ2=SZ2(J+1)
C
          LXJ=XJ2-XJ1
          LZJ=ZJ2-ZJ1
          LJ=DSQRT(LXJ**2+LZJ**2)
C
          RX1=XJ1-XI
          RX2=XJ2-XI
C
          IF(IND.EQ.1) THEN
            XN=-LZJ/LJ
            ZN=+LXJ/LJ
          ELSE
            XN=+LZJ/LJ
            ZN=-LXJ/LJ
          END IF
C
          IF(I.EQ.J) THEN
            B=LJ*SS/2.D0
            IF(B.GT.2.D0) WRITE(6,*) B,' B.GT.2'
            B3 =B**3
            B5 =B3*B*B
            B7 =B5*B*B
            YR=2.D0/PAI*((Q+DLOG(B/2.D0))*(B-B3/12.D0+B5/320.D0
     *         -B7/16128.D0)-B+B3/9.D0-17.D0/3200.D0*B5
     *         +83.D0/677376.D0*B7)
            YI=B-B3/12.D0+B5/320.D0-B7/16128.D0
            V(I,J)=-.5D0/SS/VS/VS/DN*DCMPLX(YR,YI)
            P(I,J)=(0.5D0,0.D0)
C
            RZ1=ZJ1+ZI
            R1 =DSQRT(RX1**2+RZ1**2)
            RZ2=ZJ2+ZI
            R2 =DSQRT(RX2**2+RZ2**2)
C
            A1=(XJ2-XJ1)/LJ
            A2=RX1
            A3=(ZJ2-ZJ1)/LJ
            A4=RZ1
C
            CALL GAUSPV(0.D0,LJ,NG,PSPR,PSPI,VSPR,VSPI,IER)
            PY1=SS/4.D0*PSPR
            PY2=SS/4.D0*PSPI
            VY1=-.25D0/VS/VS/DN*VSPR
            VY2=-.25D0/VS/VS/DN*VSPI
C
            AA=1.D0
            IF(NQO.NE.0)THEN
              RR=(R1+R2)/2.D0
              AA=DEXP(-OM*RR/(2.D0*QV*VS))
            END IF
            P(I,J)=P(I,J)+DCMPLX(PY1,PY2)*AA
            V(I,J)=V(I,J)+DCMPLX(VY1,VY2)*AA
C
          ELSE
C
            RZ1=ZJ1-ZI
            R1 =DSQRT(RX1**2+RZ1**2)
            RZ2=ZJ2-ZI
            R2 =DSQRT(RX2**2+RZ2**2)
C
            A1=(XJ2-XJ1)/LJ
            A2=RX1
            A3=(ZJ2-ZJ1)/LJ
            A4=RZ1
C
            RZ1M=ZJ1+ZI
            R1M=DSQRT(RX1**2+RZ1M**2)
            RZ2M=ZJ2+ZI
            R2M=DSQRT(RX2**2+RZ2M**2)
C
            A1M=(XJ2-XJ1)/LJ
            A2M=RX1
            A3M=(ZJ2-ZJ1)/LJ
            A4M=RZ1M
C
            IF(NQO.NE.0)THEN
              CALL GASPVQ(0.D0,LJ,NG,PR,PI,VR,VI,PRM,PIM,VRM,VIM,IER)
              PY1=SS/4.D0*PR
              PY1M=SS/4.D0*PRM
              PY2=SS/4.D0*PI
              PY2M=SS/4.D0*PIM
              VY1=-.25D0/VS/VS/DN*VR
              VY1M=-.25D0/VS/VS/DN*VRM
              VY2=-.25D0/VS/VS/DN*VI
              VY2M=-.25D0/VS/VS/DN*VIM
              RR=(R1+R2)/2.D0
              RRM=(R1M+R2M)/2.D0
              AA =DEXP(-OM*RR/(2.D0*QV*VS))
              AAM=DEXP(-OM*RRM/(2.D0*QV*VS))
              P(I,J)=DCMPLX(PY1,PY2)*AA+DCMPLX(PY1M,PY2M)*AAM
              V(I,J)=DCMPLX(VY1,VY2)*AA+DCMPLX(VY1M,VY2M)*AAM
C
            ELSE
C
              CALL GASPVM(0.D0,LJ,NG,PSPR,PSPI,VSPR,VSPI,IER)
              PY1=SS/4.D0*PSPR
              PY2=SS/4.D0*PSPI
              VY1=-.25D0/VS/VS/DN*VSPR
              VY2=-.25D0/VS/VS/DN*VSPI
C
              P(I,J)=DCMPLX(PY1,PY2)
              V(I,J)=DCMPLX(VY1,VY2)
C
            ENDIF
C
          ENDIF
C
  100   CONTINUE
  200 CONTINUE
C
C     *****  CHECKER'S FLAG NO.1  *****
C
      IF(NCF(1).EQ.1) THEN
      WRITE(6,*) ' *****  FROM PV22  (NCF=1) MATRIX P22    *****'
      WRITE(6,*) ((P(I,J),J=1,NI),I=1,NI)
      WRITE(6,*) ' *****  MATRIX V22    *****'
      WRITE(6,*) ((V(I,J),J=1,NI),I=1,NI)
      END IF
C
      RETURN
      END
C
C *** BESSEL FUNCTION J0(X) AND N0(X) ***
C
      SUBROUTINE BESJN0(X,BJ0,BN0)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DATA C0,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10
     * /  -2.2270614         D-9,  7.85940046        D-8,
     *    -1.9891965786      D-6,  4.03109603859     D-5,
     *    -6.449978225950    D-4,  7.9012341828039   D-3,
     *    -7.11111110046304  D-2,  4.444444444257695 D-1,
     *    -1.7777777777758446D 0,  3.9999999999998966D 0,
     *    -3.9999999999999978D 0                         /
      DATA D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13
     * /   3.315018          D-11, 4.697063          D-11,
     *    -5.66336210        D-9, -4.88429867        D-9,
     *     6.7810966225      D-7,  2.1332909858      D-7,
     *    -5.379500918201    D-5,  1.120221502384    D-5,
     *     2.52021197018629  D-3, -1.70738759690737  D-3,
     *    -5.614869347449754 D-2,  5.604047189802616 D-2,
     *     3.2757913759146509D-1, -1.7759677131433834D-1 /
      DATA E0,E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,E13,E14
     * /  -1.98650           D-12, -1.517802          D-11,
     *     3.9649790         D-10,  2.40512241        D-9,
     *    -5.707143476       D-8,  -2.5797893105      D-7,
     *     5.66074614693     D-6,   1.722612595802    D-5,
     *    -3.5284910023231   D-4,  -5.9314897388202   D-4,
     *     1.179012935710246 D-2,   6.39612989784324  D-3,
     *    -1.5037412265137390D-1,   4.68282348234590  D-3,
     *     3.0007927051955560D-1                          /
      DATA A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10
     * /   2.8228605          D-10, -1.94155134      D-9,
     *     6.51602613         D-9,  -1.533982923     D-8,
     *     3.267679187        D-8,  -8.099284883     D-8,
     *     2.8872150072       D-7,  -1.74249896802   D-6,
     *     2.184678334973     D-5,  -8.7658215906110 D-4,
     *     0.79788456080286526D 0                         /
      DATA B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12
     * /  -8.9331052         D-10,  6.65327442      D-9,
     *    -2.322479945       D- 8,  5.187839807     D-8,
     *    -8.765163373       D- 8,  1.3089752200    D-7,
     *    -2.0632292049      D- 7,  4.0856346318    D-7,
     *    -1.15933062842     D- 6,  5.25865147508   D-6,
     *    -4.423973781374    D- 5,  9.1310641570497 D-4,
     *    -9.973557010035815 D- 2                         /
      DATA H0,H1,H2,H3,H4,H5,H6,H7,H8,H9,H10,H11
     * /   3.4432403         D-9,  -1.175592825       D-7,
     *     2.8512737759      D-6,  -5.49344073691     D-5,
     *     8.276591200190    D-4,  -9.4202584419331   D-3,
     *     7.72372081288392  D-2,  -4.261441276304766 D-1,
     *     1.4216343895861494D 0,  -2.3498510134188587D 0,
     *     1.0766114686838581D 0,   3.6746690519661597D-1  /
      DATA G0,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15
     * /   1.50410           D-12, -5.72117           D-12,
     *     4.050146          D-11, -6.4960948         D-10,
     *     5.3845950         D-10,  5.922131699       D- 8,
     *     6.874770130       D- 8, -6.56197858169     D- 6,
     *     5.9100539363      D- 7,  3.9226240684696   D- 4,
     *    -1.5291708227790   D- 4, -1.348097236217831 D- 2,
     *     1.238842781168767 D- 2,  1.6904512696363842D- 1,
     *    -1.4786314339122686D- 1, -3.0851762524903378D- 1  /
      DATA F0,F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,F13,F14
     * /  -1.14460           D-12,  2.930515          D-11,
     *     1.9780007         D-10, -4.93772861        D- 9,
     *    -2.641221656       D- 8,  5.9964479378      D- 7,
     *     2.24627737107     D- 6, -4.765550355443    D- 5,
     *    -1.1239446152872   D- 4,  2.24134018235420  D- 3,
     *     2.36753744518725  D- 3, -4.900342975478215 D- 2,
     *    -8.64421637526570  D- 3,  3.0266723702418476D- 1,
     *    -2.5949743967209265D- 2                           /
      DATA PI2,PI4
     *   /  0.63661 97723 67581 34D0, 0.78539 81633 97448 31D0 /
C
      AX=DABS(X)
      IF (AX .EQ. 0.0D0 ) GO TO 50
      IF (AX .GE. 1.0D12) GO TO 40
      IF (AX .GE. 8.0D0 ) GO TO 30
      IF (AX .GE. 6.0D0 ) GO TO 20
      IF (AX .GE. 4.0D0 ) GO TO 10
C
C     !X! < 4
C
      Z=X*X*0.0625D0
      BJ0=((((((((((C0*Z +C1)*Z +C2)*Z +C3)*Z +C4)*Z +C5)*Z +C6)*Z
     *    +C7)*Z +C8)*Z +C9)*Z +C10)*Z + 1.0D0
      BN0=((((((((((H0*Z +H1)*Z +H2)*Z +H3)*Z +H4)*Z +H5)*Z +H6)*Z
     *    +H7)*Z +H8)*Z +H9)*Z +H10)*Z +H11
      BN0=PI2*BJ0*DLOG(AX*0.5D0) +BN0
      GO TO 9999
C
C     4 <= !X! < 6
   10 Z=AX-5.0D0
      BJ0=((((((((((((D0*Z+ D1)*Z+ D2)*Z+ D3)*Z+ D4)*Z+ D5)*Z
     *       +D6)*Z+ D7)*Z+ D8)*Z+ D9)*Z+ D10)*Z+ D11)*Z+ D12)*Z+ D13
      BN0=((((((((((((((G0*Z +G1)*Z +G2)*Z +G3)*Z +G4)*Z +G5)*Z
     *       +G6)*Z +G7)*Z +G8)*Z +G9)*Z +G10)*Z +G11)*Z +G12)*Z
     *       +G13)*Z +G14)*Z +G15
      GO TO 9999
C
C     6 <= !X! < 8
C
   20 Z=AX-7.0D0
      BJ0=(((((((((((((E0*Z+ E1)*Z+ E2)*Z+ E3)*Z+ E4)*Z+ E5)*Z
     *       + E6)*Z+ E7)*Z+ E8)*Z+ E9)*Z+ E10)*Z+ E11)*Z+ E12)*Z
     *       + E13)*Z+ E14
      BN0=(((((((((((((F0*Z +F1)*Z +F2)*Z +F3)*Z +F4)*Z +F5)*Z
     *       +F6)*Z +F7)*Z +F8)*Z +F9)*Z +F10)*Z +F11)*Z +F12)*Z
     *       +F13)*Z +F14
      GO TO 9999
C
C     8 <= !X! < 1.0D12
C
   30 Z=64.0D0/(AX*AX)
      A=(((((((((A0*Z +A1)*Z +A2)*Z +A3)*Z +A4)*Z +A5)*Z +A6)*Z
     *  +A7)*Z +A8)*Z +A9)*Z +A10
      B=((((((((((((B0*Z +B1)*Z +B2)*Z +B3)*Z +B4)*Z +B5)*Z +B6)*Z
     *  +B7)*Z +B8)*Z +B9)*Z +B10)*Z +B11)*Z +B12) /AX
      Q=AX-PI4
      CC= DCOS(Q)
      CS= DSIN(Q)
      CQ= DSQRT(AX)
      BJ0=(A*CC - B*CS) /CQ
      BN0=(A*CS + B*CC) /CQ
      GO TO 9999
C
C     PARAMETER ERROR
C
   40 BJ0= 0.0D0
      BN0= 0.0D0
      GO TO 9999
C
   50 BJ0=  1.0D 0
C      BN0= -1.0D75
      BN0= -1.0D35
C
 9999 RETURN
      END
C
C *** BESSEL FUNCTION J1(X), N1(X)
C
      SUBROUTINE BESJN1(X,BJ1,BN1)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DATA C0,C1,C2,C3,C4,C5,C6,C7,C8,C9
     * /  1.22494472        D-8, -3.929965711       D-7,
     *    8.9514968720      D-6, -1.612440490740    D-4,
     *    2.2574925903842   D-3, -2.37037026763209  D-2,
     *    1.777777775582676 D-1, -8.888888888615124 D-1,
     *    2.6666666666649150D 0, -3.9999999999999563D 0  /
      DATA D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11
     * /  2.3353113         D-9, -8.71498135        D-8,
     *    2.3316556575      D-6, -5.00978325295     D-5,
     *    8.533226991090    D-4, -1.11766285197065  D-2,
     *    1.080130262837708 D-1, -7.271013435638566 D-1,
     *    3.1262109002217823D 0, -7.3980378533256488D 0,
     *    6.8529249642057360D 0,  3.9325615210260295D-1  /
      DATA E0,E1,E2,E3,E4,E5,E6,E7,E8,E9,E10,E11,E12,E13,E14,E15
     * /   4.84682           D-12, -2.280090          D-11,
     *     6.313256          D-11, -5.2574110         D-10,
     *     7.81894200        D- 9, -5.92404185        D- 9,
     *    -5.9222983079      D- 7, -6.1872869462      D- 7,
     *     5.249583490136    D- 5, -4.13703795176     D- 6,
     *    -2.35357444227453  D- 3,  7.6458541141893   D- 4,
     *     5.392388944881264 D- 2, -3.716528343506466 D- 2,
     *    -3.3809025392727920D- 1,  1.4786314339122685D- 1  /
      DATA F0,F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12,F13,F14
     * /   1.78737           D-12,  1.601135          D-11,
     *    -3.8677585         D-10, -2.37356163        D- 9,
     *     5.432238759       D- 8,  2.6412212075      D- 7,
     *    -5.39680775206     D- 6, -1.797021894410    D- 5,
     *     3.3358852634720   D- 4,  6.7436676916589   D- 4,
     *    -1.120670091199092 D- 2, -9.47014978074827  D- 3,
     *     1.4701028926435867D- 1,  1.728843275053138 D- 2,
     *    -3.0266723702418487D- 1                           /
      DATA G0,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14
     * /  -1.07463           D-12,  2.778003          D-11,
     *     2.0080 677        D-10, -4.75788198        D- 9,
     *    -2.6460 77937      D- 8,  5.7071424129      D- 7,
     *     2.32181314997     D- 6, -4.528596911743    D- 5,
     *    -1.2058288258765   D- 4,  2.11709460137864  D- 3,
     *     2.96574486954234  D- 3, -4.716051742840817 D- 2,
     *    -1.918838969353708 D- 2,  3.0074824530274775D- 1,
     *    -4.68282348234583  D- 3                           /
      DATA H0,H1,H2,H3,H4,H5,H6,H7,H8,H9,H10,H11,H12,H13,H14
     * /   2.21931           D-12,  4.11634           D-12,
     *    -4.3816515         D-10, -5.7599652         D-10,
     *     6.230613774       D- 8,  4.885713666       D- 8,
     *    -6.10299268189     D- 6, -1.70664050678     D- 6,
     *     3.7656506609458   D- 4, -6.721328811703    D- 5,
     *    -1.260105985120454 D- 2,  6.82955038740438  D- 3,
     *     1.6844608042350780D- 1, -1.1208094379604528D- 1,
     *    -3.2757913759146522D- 1                           /
      DATA A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10
     * /  -3.0207732         D-10,  2.08062912       D-9,
     *    -7.00035699        D- 9,  1.655768548      D-8,
     *    -3.558232648       D- 8,  8.955040706      D-8,
     *    -3.2722783265      D- 7,  2.05931893593    D-6,
     *    -2.808872167333    D- 5,  1.46097026511633 D-3,
     *     7.9788456080286546D- 1                           /
      DATA B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12
     * /   9.4650569         D-10, -7.05369184       D-9,
     *     2.464612368       D- 8, -5.514232274      D-8,
     *     9.343285420       D- 8, -1.4024484043     D-7,
     *     2.2297669562      D- 7, -4.4752033584     D-7,
     *     1.29573186799     D- 6, -6.06767597238    D-6,
     *     5.407079074732    D- 5, -1.27834898198997 D-3,
     *     2.9920671030107448D- 1                           /
      DATA PI2,PI4
     *   /  0.63661 97723 67581 34D0, 0.78539 81633 97448 31D0 /
C
      AX=DABS(X)
C      IF (AX .LT. 8.797D-77) GO TO 60
      IF (AX .LT. 8.797D-37) GO TO 60
      IF (AX .GE. 1.0D12) GO TO 50
      IF (AX .GE. 8.0D0) GO TO 30
      IF (AX .GE. 6.0D0) GO TO 20
      IF (AX .GE. 4.0D0) GO TO 10
C
C     0 < !X! < 4
C
      Y=X*0.25D0
      Z=Y*Y
      BJ1=((((((((((C0*Z +C1)*Z +C2)*Z +C3)*Z +C4)*Z +C5)*Z +C6)*Z
     *       +C7)*Z +C8)*Z +C9)*Z + 2.0D0) *Y
      BN1=((((((((((((D0*Z +D1)*Z +D2)*Z +D3)*Z +D4)*Z +D5)*Z +D6)*Z
     *       +D7)*Z +D8)*Z +D9)*Z +D10)*Z +D11)*Z -PI2) /X
      BN1=PI2*BJ1*DLOG(AX*0.5D0) +BN1
      GO TO 9999
C
C     4 <= !X! < 6
C
   10 Z=X-5.0D0
      BJ1=(((((((((((((H0*Z +H1)*Z +H2)*Z +H3)*Z +H4)*Z +H5)*Z
     *       +H6)*Z +H7)*Z +H8)*Z +H9)*Z +H10)*Z +H11)*Z +H12)*Z
     *       +H13)*Z +H14
      BN1=((((((((((((((E0*Z +E1)*Z +E2)*Z +E3)*Z +E4)*Z +E5)*Z
     *       +E6)*Z +E7)*Z +E8)*Z +E9)*Z +E10)*Z +E11)*Z +E12)*Z
     *       +E13)*Z +E14)*Z +E15
      GO TO 40
C
C     6 <= !X! < 8
C
   20 Z=X-7.0D0
      BJ1=(((((((((((((G0*Z +G1)*Z +G2)*Z +G3)*Z +G4)*Z +G5)*Z
     *       +G6)*Z +G7)*Z +G8)*Z +G9)*Z +G10)*Z +G11)*Z +G12)*Z
     *       +G13)*Z +G14
      BN1=(((((((((((((F0*Z +F1)*Z +F2)*Z +F3)*Z +F4)*Z +F5)*Z
     *       +F6)*Z +F7)*Z +F8)*Z +F9)*Z +F10)*Z +F11)*Z +F12)*Z
     *       +F13)*Z +F14
      GO TO 40
C
C     8 <= !X! < 1.0D12
C
   30 Z=64.0D0/(X*X)
      A=(((((((((A0*Z +A1)*Z +A2)*Z +A3)*Z +A4)*Z +A5)*Z +A6)*Z +A7)*Z
     *       +A8)*Z +A9)*Z +A10
      B=((((((((((((B0*Z +B1)*Z +B2)*Z +B3)*Z +B4)*Z +B5)*Z +B6)*Z
     *       +B7)*Z +B8)*Z +B9)*Z +B10)*Z +B11)*Z +B12) /AX
      Q=X-PI4
      CC= DCOS(Q)
      CS= DSIN(Q)
      CQ= DSQRT(AX)
      BJ1=(A*CS + B*CC) /CQ
      BN1=(B*CS - A*CC) /CQ
C
   40 IF (AX .LT. 0.0D0) THEN
      BJ1=-BJ1
      BN1=-BN1
      END IF
      GO TO 9999
C
C     PARAMETER ERROR
C
   50 BJ1= 0.D0
      BN1= 0.D0
      GO TO 9999
C
   60 BJ1= 0.0D0
C      BN1=-1.0D75
      BN1=-1.0D35
C
 9999 RETURN
      END
