C     Last change:  PGM   8 Nov 2000    1:04 pm
C     PROGRAM SOMNEC(INPUT,OUTPUT,TAPE21)
C
C     PROGRAM TO GENERATE NEC INTERPOLATION GRIDS FOR FIELDS DUE TO
C     GROUND.  FIELD COMPONENTS ARE COMPUTED BY NUMERICAL EVALUATION
C     OF MODIFIED SOMMERFELD INTEGRALS.
C
C     SOMNEC2D IS A DOUBLE PRECISION VERSION OF SOMNEC FOR USE WITH
C     NEC2D.  AN ALTERNATE VERSION (SOMNEC2SD) IS ALSO PROVIDED IN WHICH
C     COMPUTATION IS IN SINGLE PRECISION BUT THE OUTPUT FILE IS WRITTEN
C     IN DOUBLE PRECISION FOR USE WITH NEC2D.  SOMNEC2SD RUNS ABOUT TWIC
C     AS FAST AS THE FULL DOUBLE PRECISION SOMNEC2D.  THE DIFFERENCE
C     BETWEEN NEC2D RESULTS USING A FOR021 FILE FROM THIS CODE RATHER
C     THAN FROM SOMNEC2SD WAS INSIGNFICANT IN THE CASES TESTED.
C
C     Changes made by J Bergervoet, 31-5-95:
C         Parameter 0. --> 0.D0 in calling of routine TEST
C         Status of output files set to 'UNKNOWN'
C***
      IMPLICIT REAL*8(A-H,O-Z)
C***
      COMPLEX*16 CK1,CK1SQ,ERV,EZV,ERH,EPH,CKSM,CT1,CT2,CT3,CL1,CL2,CON,
     1AR1,AR2,AR3,EPSCF
      COMMON /EVLCOM/ CKSM,CT1,CT2,CT3,CK1,CK1SQ,CK2,CK2SQ,TKMAG,TSMAG,C
     1K1R,ZPH,RHO,JH
      COMMON /GGRID/ AR1(11,10,4),AR2(17,5,4),AR3(9,8,4),EPSCF,DXA(3),DY
     1A(3),XSA(3),YSA(3),NXA(3),NYA(3)
      CHARACTER*3  LCOMP(4)
      DATA LCOMP/'ERV','EZV','ERH','EPH'/
      WRITE(*,*) 'SOMNEC2D,  Last changes: May 31 1995,  J. Bergervoet'
      WRITE(*,*)
C
      Write(*,*)
     &'GIVE GROUND PARAMETERS - EPR = RELATIVE DIELECTRIC CONSTANT'
      Write(*,*) '                         SIG = CONDUCTIVITY (MHOS/M)'
      Write(*,*) '                         FMHZ = FREQUENCY (MHZ)'
      Write(*,*)
     &'                         IPT = 1 TO PRINT GRIDS.  =0 OTHERWISE.'
      Write(*,*)
     &'IF SIG .LT. 0. THEN COMPLEX DIELECTRIC CONSTANT = EPR + J*SIG'
      Write(*,*) 'AND FMHZ IS NOT USED.'
C
999   WRITE(*,21)
C             21    FORMAT($,' ENTER EPR,SIG,FMHZ,IPT > ')
      READ(*,*,ERR=999) EPR,SIG,FMHZ,IPT
      WRITE(*,22)
C             22    FORMAT(' STARTING COMPUTATION OF SOMMERFELD INTEGRAL TABLES')
      WRITE(*,*)
      WRITE(*,*)

      WRITE(*,100) EPR
100   FORMAT("  RELATIVE DIELECTRIC CONSTANT (EPR)  = ", D20.5)
      WRITE(*,101) SIG
101   FORMAT("  SIGMA [CONDUCTIVITY IN MHOS/METER]  = ", D20.5)
      WRITE(*,102) FMHZ
102   FORMAT("                     FREQUENCY IN MHZ = ", D20.5)
      IF(IPT == 1) WRITE(*,*) "   GRID FILE [SOM2D.OUT] WILL BE CREATED"
      IF(IPT == 0) WRITE(*,*) "   NO GRID FILE WILL BE CREATED"
      WRITE(*,*)
C***
      IF (SIG.LT.0.) GO TO 1
      WLAM=299.8/FMHZ
      EPSCF=DCMPLX(EPR,-SIG*WLAM*59.96)
      GO TO 2
1     EPSCF=DCMPLX(EPR,SIG)
2     CALL SECOND (TST)
      CK2=6.283185308
      CK2SQ=CK2*CK2
C
C     SOMMERFELD INTEGRAL EVALUATION USES EXP(-JWT), NEC USES EXP(+JWT),
C     HENCE NEED CONJG(EPSCF).  CONJUGATE OF FIELDS OCCURS IN SUBROUTINE
C     EVLUA.
C
      CK1SQ=CK2SQ*DCONJG(EPSCF)
      CK1=SQRT(CK1SQ)
      CK1R=DREAL(CK1)
      TKMAG=100.*ABS(CK1)
      TSMAG=100.*CK1*DCONJG(CK1)
      CKSM=CK2SQ/(CK1SQ+CK2SQ)
      CT1=.5*(CK1SQ-CK2SQ)
      ERV=CK1SQ*CK1SQ
      EZV=CK2SQ*CK2SQ
      CT2=.125*(ERV-EZV)
      ERV=ERV*CK1SQ
      EZV=EZV*CK2SQ
      CT3=.0625*(ERV-EZV)
C
C     LOOP OVER 3 GRID REGIONS
C
      DO 6 K=1,3
      NR=NXA(K)
      NTH=NYA(K)
      DR=DXA(K)
      DTH=DYA(K)
      R=XSA(K)-DR
      IRS=1
      IF (K.EQ.1) R=XSA(K)
      IF (K.EQ.1) IRS=2
C
C     LOOP OVER R.  (R=SQRT(RHO**2 + (Z+H)**2))
C
      DO 6 IR=IRS,NR
      R=R+DR
      THET=YSA(K)-DTH
C
C     LOOP OVER THETA.  (THETA=ATAN((Z+H)/RHO))
C
      DO 6 ITH=1,NTH
      THET=THET+DTH
      RHO=R*COS(THET)
      ZPH=R*SIN(THET)
      IF (RHO.LT.1.E-7) RHO=1.E-8
      IF (ZPH.LT.1.E-7) ZPH=0.
      CALL EVLUA (ERV,EZV,ERH,EPH)
      RK=CK2*R
      CON=-(0.,4.77147)*R/DCMPLX(COS(RK),-SIN(RK))
      GO TO (3,4,5), K
3     AR1(IR,ITH,1)=ERV*CON
      AR1(IR,ITH,2)=EZV*CON
      AR1(IR,ITH,3)=ERH*CON
      AR1(IR,ITH,4)=EPH*CON
      GO TO 6
4     AR2(IR,ITH,1)=ERV*CON
      AR2(IR,ITH,2)=EZV*CON
      AR2(IR,ITH,3)=ERH*CON
      AR2(IR,ITH,4)=EPH*CON
      GO TO 6
5     AR3(IR,ITH,1)=ERV*CON
      AR3(IR,ITH,2)=EZV*CON
      AR3(IR,ITH,3)=ERH*CON
      AR3(IR,ITH,4)=EPH*CON
6     CONTINUE
C
C     FILL GRID 1 FOR R EQUAL TO ZERO.
C
      CL2=-(0.,188.370)*(EPSCF-1.)/(EPSCF+1.)
      CL1=CL2/(EPSCF+1.)
      EZV=EPSCF*CL1
      THET=-DTH
      NTH=NYA(1)
      DO 9 ITH=1,NTH
      THET=THET+DTH
      IF (ITH.EQ.NTH) GO TO 7
      TFAC2=COS(THET)
      TFAC1=(1.-SIN(THET))/TFAC2
      TFAC2=TFAC1/TFAC2
      ERV=EPSCF*CL1*TFAC1
      ERH=CL1*(TFAC2-1.)+CL2
      EPH=CL1*TFAC2-CL2
      GO TO 8
7     ERV=0.
      ERH=CL2-.5*CL1
      EPH=-ERH
8     AR1(1,ITH,1)=ERV
      AR1(1,ITH,2)=EZV
      AR1(1,ITH,3)=ERH
9     AR1(1,ITH,4)=EPH
      CALL SECOND (TIM)
C
C     WRITE GRID ON TAPE21
C
      OPEN(UNIT=21,FILE='SOM2D.NEC',STATUS='UNKNOWN',FORM='UNFORMATTED')
      WRITE (21) AR1,AR2,AR3,EPSCF,DXA,DYA,XSA,YSA,NXA,NYA
      REWIND 21
      IF (IPT.EQ.0) GO TO 14
C
C     PRINT GRID
C
C	DEBUGGING CODE
C      ---------------------------------------------------
      PRINT *,'AR1(1,1,1)= ',AR1(1,1,1)
C      PRINT *,'AR2(1,1,1)= ',AR2(1,1,1)
C      PRINT *,'AR3(1,1,1)= ',AR3(1,1,1)
      PRINT *,'EPSCF= ',EPSCF
      PRINT *,'DXA= ',DXA
      PRINT *,'DYA= ',DYA
      PRINT *,'XSA= ',XSA
      PRINT *,'YSA= ',YSA
      PRINT *,'NXA= ',NXA
      PRINT *,'NYA= ',NYA
      PRINT *,'AR1= ',AR1
      PRINT *,'AR2= ',AR2
      PRINT *,'AR3= ',AR3
      PRINT 444,AR1(1,1,1)
444   FORMAT(11HAR1(1,1,1)=,E12.5)	

C      ---------------------------------------------------

      OPEN (UNIT=3,FILE='SOM2D.OUT',STATUS='NEW',ERR=14)
      WRITE(3,17) EPSCF
      DO 13 K=1,3
      NR=NXA(K)
      NTH=NYA(K)
      WRITE(3,18) K,XSA(K),DXA(K),NR,YSA(K),DYA(K),NTH
      DO 13 L=1,4
      WRITE(3,19) LCOMP(L)
      DO 13 IR=1,NR
      GO TO (10,11,12), K
10    WRITE(3,20) IR,(AR1(IR,ITH,L),ITH=1,NTH)
      GO TO 13
11    WRITE(3,20) IR,(AR2(IR,ITH,L),ITH=1,NTH)
      GO TO 13
12    WRITE(3,20) IR,(AR3(IR,ITH,L),ITH=1,NTH)
13    CONTINUE
14    TIM=TIM-TST
      WRITE(*,16) TIM
      STOP
C
16    FORMAT (6H TIME=,1PE12.5)
17    FORMAT (30H1NEC GROUND INTERPOLATION GRID,/,21H DIELECTRIC CONSTAN
     1T=,1P2E12.5)
18    FORMAT (///,5H GRID,I2,/,4X,5HR(1)=,F7.4,4X,3HDR=,F7.4,4X,3HNR=,I3
     1,/,9H THET(1)=,F7.4,3X,4HDTH=,F7.4,3X,4HNTH=,I3,//)
19    FORMAT (///,1X,A3)
20    FORMAT (4H IR=,I3,/,1X,(10E12.5))
21    FORMAT($,' ENTER EPR,SIG,FMHZ,IPT > ')
22    FORMAT(' STARTING COMPUTATION OF SOMMERFELD INTEGRAL TABLES')
      END
      BLOCK DATA SOMSET
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 AR1,AR2,AR3,EPSCF
      COMMON /GGRID/ AR1(11,10,4),AR2(17,5,4),AR3(9,8,4),EPSCF,DXA(3),DY
     1A(3),XSA(3),YSA(3),NXA(3),NYA(3)
      DATA NXA/11,17,9/,NYA/10,5,8/,XSA/0.,.2,.2/,YSA/0.,0.,.3490658504/
      DATA DXA/.02,.05,.1/,DYA/.1745329252,.0872664626,.1745329252/
      END
      SUBROUTINE BESSEL (Z,J0,J0P)
C
C     BESSEL EVALUATES THE ZERO-ORDER BESSEL FUNCTION AND ITS DERIVATIVE
C     FOR COMPLEX ARGUMENT Z.
C
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      COMPLEX*16 J0,J0P,P0Z,P1Z,Q0Z,Q1Z,Z,ZI,ZI2,ZK,FJ,CZ,SZ,J0X,J0PX
      DIMENSION M(101), A1(25), A2(25), FJX(2)
      EQUIVALENCE (FJ,FJX)
      DATA PI,C3,P10,P20,Q10,Q20/3.141592654,.7978845608,.0703125,.11215
     120996,.125,.0732421875/
      DATA P11,P21,Q11,Q21/.1171875,.1441955566,.375,.1025390625/
      DATA POF,INIT/.7853981635,0/,FJX/0.,1./
      IF (INIT.EQ.0) GO TO 5
1     ZMS=Z*DCONJG(Z)
      IF (ZMS.GT.1.E-12) GO TO 2
      J0=(1.,0.)
      J0P=-.5*Z
      RETURN
2     IB=0
      IF (ZMS.GT.37.21) GO TO 4
      IF (ZMS.GT.36.) IB=1
C     SERIES EXPANSION
      IZ=1.+ZMS
      MIZ=M(IZ)
      J0=(1.,0.)
      J0P=J0
      ZK=J0
      ZI=Z*Z
      DO 3 K=1,MIZ
      ZK=ZK*A1(K)*ZI
      J0=J0+ZK
3     J0P=J0P+A2(K)*ZK
      J0P=-.5*Z*J0P
      IF (IB.EQ.0) RETURN
      J0X=J0
      J0PX=J0P
C     ASYMPTOTIC EXPANSION
4     ZI=1./Z
      ZI2=ZI*ZI
      P0Z=1.+(P20*ZI2-P10)*ZI2
      P1Z=1.+(P11-P21*ZI2)*ZI2
      Q0Z=(Q20*ZI2-Q10)*ZI
      Q1Z=(Q11-Q21*ZI2)*ZI
      ZK=EXP(FJ*(Z-POF))
      ZI2=1./ZK
      CZ=.5*(ZK+ZI2)
      SZ=FJ*.5*(ZI2-ZK)
      ZK=C3*SQRT(ZI)
      J0=ZK*(P0Z*CZ-Q0Z*SZ)
      J0P=-ZK*(P1Z*SZ+Q1Z*CZ)
      IF (IB.EQ.0) RETURN
      ZMS=COS((SQRT(ZMS)-6.)*31.41592654)
      J0=.5*(J0X*(1.+ZMS)+J0*(1.-ZMS))
      J0P=.5*(J0PX*(1.+ZMS)+J0P*(1.-ZMS))
      RETURN
C     INITIALIZATION OF CONSTANTS
5     DO 6 K=1,25
      A1(K)=-.25D0/(K*K)
6     A2(K)=1.D0/(K+1.D0)
      DO 8 I=1,101
      TEST=1.D0
      DO 7 K=1,24
      INIT=K
      TEST=-TEST*I*A1(K)
      IF (TEST.LT.1.D-6) GO TO 8
7     CONTINUE
8     M(I)=INIT
      GO TO 1
      END
      SUBROUTINE EVLUA (ERV,EZV,ERH,EPH)
C
C     EVALUA CONTROLS THE INTEGRATION CONTOUR IN THE COMPLEX LAMBDA
C     PLANE FOR EVALUATION OF THE SOMMERFELD INTEGRALS.
C
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      COMPLEX*16 ERV,EZV,ERH,EPH,A,B,CK1,CK1SQ,BK,SUM,DELTA,ANS,DELTA2,
     1CP1,CP2,CP3,CKSM,CT1,CT2,CT3
      COMMON /CNTOUR/ A,B
      COMMON /EVLCOM/ CKSM,CT1,CT2,CT3,CK1,CK1SQ,CK2,CK2SQ,TKMAG,TSMAG,C
     1K1R,ZPH,RHO,JH
      DIMENSION SUM(6), ANS(6)
      DATA PTP/.6283185308/
      DEL=ZPH
      IF (RHO.GT.DEL) DEL=RHO
      IF (ZPH.LT.2.*RHO) GO TO 4
C
C     BESSEL FUNCTION FORM OF SOMMERFELD INTEGRALS
C
      JH=0
      A=(0.,0.)
      DEL=1./DEL
      IF (DEL.LE.TKMAG) GO TO 2
      B=DCMPLX(.1*TKMAG,-.1*TKMAG)
      CALL ROM1 (6,SUM,2)
      A=B
      B=DCMPLX(DEL,-DEL)
      CALL ROM1 (6,ANS,2)
      DO 1 I=1,6
1     SUM(I)=SUM(I)+ANS(I)
      GO TO 3
2     B=DCMPLX(DEL,-DEL)
      CALL ROM1 (6,SUM,2)
3     DELTA=PTP*DEL
      CALL GSHANK (B,DELTA,ANS,6,SUM,0,B,B)
      GO TO 10
C
C     HANKEL FUNCTION FORM OF SOMMERFELD INTEGRALS
C
4     JH=1
      CP1=DCMPLX(0.D0,.4*CK2)
      CP2=DCMPLX(.6*CK2,-.2*CK2)
      CP3=DCMPLX(1.02*CK2,-.2*CK2)
      A=CP1
      B=CP2
      CALL ROM1 (6,SUM,2)
      A=CP2
      B=CP3
      CALL ROM1 (6,ANS,2)
      DO 5 I=1,6
5     SUM(I)=-(SUM(I)+ANS(I))
C     PATH FROM IMAGINARY AXIS TO -INFINITY
      SLOPE=1000.
      IF (ZPH.GT..001*RHO) SLOPE=RHO/ZPH
      DEL=PTP/DEL
      DELTA=DCMPLX(-1.D0,SLOPE)*DEL/SQRT(1.+SLOPE*SLOPE)
      DELTA2=-DCONJG(DELTA)
      CALL GSHANK (CP1,DELTA,ANS,6,SUM,0,BK,BK)
      RMIS=RHO*(DREAL(CK1)-CK2)
      IF (RMIS.LT.2.*CK2) GO TO 8
      IF (RHO.LT.1.E-10) GO TO 8
      IF (ZPH.LT.1.E-10) GO TO 6
      BK=DCMPLX(-ZPH,RHO)*(CK1-CP3)
      RMIS=-DREAL(BK)/ABS(DIMAG(BK))
      IF(RMIS.GT.4.*RHO/ZPH)GO TO 8
C     INTEGRATE UP BETWEEN BRANCH CUTS, THEN TO + INFINITY
6     CP1=CK1-(.1,.2)
      CP2=CP1+.2
      BK=DCMPLX(0.D0,DEL)
      CALL GSHANK (CP1,BK,SUM,6,ANS,0,BK,BK)
      A=CP1
      B=CP2
      CALL ROM1 (6,ANS,1)
      DO 7 I=1,6
7     ANS(I)=ANS(I)-SUM(I)
      CALL GSHANK (CP3,BK,SUM,6,ANS,0,BK,BK)
      CALL GSHANK (CP2,DELTA2,ANS,6,SUM,0,BK,BK)
      GO TO 10
C     INTEGRATE BELOW BRANCH POINTS, THEN TO + INFINITY
8     DO 9 I=1,6
9     SUM(I)=-ANS(I)
      RMIS=DREAL(CK1)*1.01
      IF (CK2+1..GT.RMIS) RMIS=CK2+1.
      BK=DCMPLX(RMIS,.99*DIMAG(CK1))
      DELTA=BK-CP3
      DELTA=DELTA*DEL/ABS(DELTA)
      CALL GSHANK (CP3,DELTA,ANS,6,SUM,1,BK,DELTA2)
10    ANS(6)=ANS(6)*CK1
C     CONJUGATE SINCE NEC USES EXP(+JWT)
      ERV=DCONJG(CK1SQ*ANS(3))
      EZV=DCONJG(CK1SQ*(ANS(2)+CK2SQ*ANS(5)))
      ERH=DCONJG(CK2SQ*(ANS(1)+ANS(6)))
      EPH=-DCONJG(CK2SQ*(ANS(4)+ANS(6)))
      RETURN
      END
      SUBROUTINE GSHANK (START,DELA,SUM,NANS,SEED,IBK,BK,DELB)
C
C     GSHANK INTEGRATES THE 6 SOMMERFELD INTEGRALS FROM START TO
C     INFINITY (UNTIL CONVERGENCE) IN LAMBDA.  AT THE BREAK POINT, BK,
C     THE STEP INCREMENT MAY BE CHANGED FROM DELA TO DELB.  SHANK S
C     ALGORITHM TO ACCELERATE CONVERGENCE OF A SLOWLY CONVERGING SERIES
C     IS USED
C
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      COMPLEX*16 START,DELA,SUM,SEED,BK,DELB,A,B,Q1,Q2,ANS1,ANS2,A1,A2,
     1AS1,AS2,DEL,AA
      COMMON /CNTOUR/ A,B
      DIMENSION Q1(6,20), Q2(6,20), ANS1(6), ANS2(6), SUM(6), SEED(6)
      DATA CRIT/1.E-4/,MAXH/20/
      RBK=DREAL(BK)
      DEL=DELA
      IBX=0
      IF (IBK.EQ.0) IBX=1
      DO 1 I=1,NANS
1     ANS2(I)=SEED(I)
      B=START
2     DO 20 INT=1,MAXH
      INX=INT
      A=B
      B=B+DEL
      IF (IBX.EQ.0.AND.DREAL(B).GE.RBK) GO TO 5
      CALL ROM1 (NANS,SUM,2)
      DO 3 I=1,NANS
3     ANS1(I)=ANS2(I)+SUM(I)
      A=B
      B=B+DEL
      IF (IBX.EQ.0.AND.DREAL(B).GE.RBK) GO TO 6
      CALL ROM1 (NANS,SUM,2)
      DO 4 I=1,NANS
4     ANS2(I)=ANS1(I)+SUM(I)
      GO TO 11
C     HIT BREAK POINT.  RESET SEED AND START OVER.
5     IBX=1
      GO TO 7
6     IBX=2
7     B=BK
      DEL=DELB
      CALL ROM1 (NANS,SUM,2)
      IF (IBX.EQ.2) GO TO 9
      DO 8 I=1,NANS
8     ANS2(I)=ANS2(I)+SUM(I)
      GO TO 2
9     DO 10 I=1,NANS
10    ANS2(I)=ANS1(I)+SUM(I)
      GO TO 2
11    DEN=0.
      DO 18 I=1,NANS
      AS1=ANS1(I)
      AS2=ANS2(I)
      IF (INT.LT.2) GO TO 17
      DO 16 J=2,INT
      JM=J-1
      AA=Q2(I,JM)
      A1=Q1(I,JM)+AS1-2.*AA
      IF (DREAL(A1).EQ.0..AND.DIMAG(A1).EQ.0.) GO TO 12
      A2=AA-Q1(I,JM)
      A1=Q1(I,JM)-A2*A2/A1
      GO TO 13
12    A1=Q1(I,JM)
13    A2=AA+AS2-2.*AS1
      IF (DREAL(A2).EQ.0..AND.DIMAG(A2).EQ.0.) GO TO 14
      A2=AA-(AS1-AA)*(AS1-AA)/A2
      GO TO 15
14    A2=AA
15    Q1(I,JM)=AS1
      Q2(I,JM)=AS2
      AS1=A1
16    AS2=A2
17    Q1(I,INT)=AS1
      Q2(I,INT)=AS2
      AMG=ABS(DREAL(AS2))+ABS(DIMAG(AS2))
      IF (AMG.GT.DEN) DEN=AMG
18    CONTINUE
      DENM=1.E-3*DEN*CRIT
      JM=INT-3
      IF (JM.LT.1) JM=1
      DO 19 J=JM,INT
      DO 19 I=1,NANS
      A1=Q2(I,J)
      DEN=(ABS(DREAL(A1))+ABS(DIMAG(A1)))*CRIT
      IF (DEN.LT.DENM) DEN=DENM
      A1=Q1(I,J)-A1
      AMG=ABS(DREAL(A1))+ABS(DIMAG(A1))
      IF (AMG.GT.DEN) GO TO 20
19    CONTINUE
      GO TO 22
20    CONTINUE
      WRITE(*,24)
      DO 21 I=1,NANS
21    WRITE(*,25) Q1(I,INX),Q2(I,INX)
22    DO 23 I=1,NANS
23    SUM(I)=.5*(Q1(I,INX)+Q2(I,INX))
      RETURN
C
24    FORMAT (46H **** NO CONVERGENCE IN SUBROUTINE GSHANK ****)
25    FORMAT (1X,1P10E12.5)
      END
      SUBROUTINE HANKEL (Z,H0,H0P)
C
C     HANKEL EVALUATES HANKEL FUNCTION OF THE FIRST KIND, ORDER ZERO,
C     AND ITS DERIVATIVE FOR COMPLEX ARGUMENT Z.
C
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      COMPLEX*16 CLOGZ,H0,H0P,J0,J0P,P0Z,P1Z,Q0Z,Q1Z,Y0,Y0P,Z,ZI,ZI2,ZK,
     1FJ
      DIMENSION M(101), A1(25), A2(25), A3(25), A4(25), FJX(2)
      EQUIVALENCE (FJ,FJX)
      DATA PI,GAMMA,C1,C2,C3,P10,P20/3.141592654,.5772156649,-.024578509
     15,.3674669052,.7978845608,.0703125,.1121520996/
      DATA Q10,Q20,P11,P21,Q11,Q21/.125,.0732421875,.1171875,.1441955566
     1,.375,.1025390625/
      DATA POF,INIT/.7853981635,0/,FJX/0.,1./
      IF (INIT.EQ.0) GO TO 5
1     ZMS=Z*DCONJG(Z)
      IF (ZMS.NE.0.) GO TO 2
      WRITE(*,9)
      STOP
2     IB=0
      IF (ZMS.GT.16.81) GO TO 4
      IF (ZMS.GT.16.) IB=1
C     SERIES EXPANSION
      IZ=1.+ZMS
      MIZ=M(IZ)
      J0=(1.,0.)
      J0P=J0
      Y0=(0.,0.)
      Y0P=Y0
      ZK=J0
      ZI=Z*Z
      DO 3 K=1,MIZ
      ZK=ZK*A1(K)*ZI
      J0=J0+ZK
      J0P=J0P+A2(K)*ZK
      Y0=Y0+A3(K)*ZK
3     Y0P=Y0P+A4(K)*ZK
      J0P=-.5*Z*J0P
      CLOGZ=LOG(.5*Z)
      Y0=(2.*J0*CLOGZ-Y0)/PI+C2
      Y0P=(2./Z+2.*J0P*CLOGZ+.5*Y0P*Z)/PI+C1*Z
      H0=J0+FJ*Y0
      H0P=J0P+FJ*Y0P
      IF (IB.EQ.0) RETURN
      Y0=H0
      Y0P=H0P
C     ASYMPTOTIC EXPANSION
4     ZI=1./Z
      ZI2=ZI*ZI
      P0Z=1.+(P20*ZI2-P10)*ZI2
      P1Z=1.+(P11-P21*ZI2)*ZI2
      Q0Z=(Q20*ZI2-Q10)*ZI
      Q1Z=(Q11-Q21*ZI2)*ZI
      ZK=EXP(FJ*(Z-POF))*SQRT(ZI)*C3
      H0=ZK*(P0Z+FJ*Q0Z)
      H0P=FJ*ZK*(P1Z+FJ*Q1Z)
      IF (IB.EQ.0) RETURN
      ZMS=COS((SQRT(ZMS)-4.)*31.41592654)
      H0=.5*(Y0*(1.+ZMS)+H0*(1.-ZMS))
      H0P=.5*(Y0P*(1.+ZMS)+H0P*(1.-ZMS))
      RETURN
C     INITIALIZATION OF CONSTANTS
5     PSI=-GAMMA
      DO 6 K=1,25
      A1(K)=-.25D0/(K*K)
      A2(K)=1.D0/(K+1.D0)
      PSI=PSI+1.D0/K
      A3(K)=PSI+PSI
6     A4(K)=(PSI+PSI+1.D0/(K+1.D0))/(K+1.D0)
      DO 8 I=1,101
      TEST=1.D0
      DO 7 K=1,24
      INIT=K
      TEST=-TEST*I*A1(K)
      IF (TEST*A3(K).LT.1.D-6) GO TO 8
7     CONTINUE
8     M(I)=INIT
      GO TO 1
C
9     FORMAT (34H ERROR - HANKEL NOT VALID FOR Z=0.)
      END
      SUBROUTINE LAMBDA (T,XLAM,DXLAM)
C
C     COMPUTE INTEGRATION PARAMETER XLAM=LAMBDA FROM PARAMETER T.
C
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      COMPLEX*16 A,B,XLAM,DXLAM
      COMMON /CNTOUR/ A,B
      DXLAM=B-A
      XLAM=A+DXLAM*T
      RETURN
      END
      SUBROUTINE ROM1 (N,SUM,NX)
C
C     ROM1 INTEGRATES THE 6 SOMMERFELD INTEGRALS FROM A TO B IN LAMBDA.
C     THE METHOD OF VARIABLE INTERVAL WIDTH ROMBERG INTEGRATION IS USED.
C
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      COMPLEX*16 A,B,SUM,G1,G2,G3,G4,G5,T00,T01,T10,T02,T11,T20
      COMMON /CNTOUR/ A,B
      DIMENSION SUM(6), G1(6), G2(6), G3(6), G4(6), G5(6), T01(6), T10(6
     1), T20(6)
      DATA NM,NTS,RX/131072,4,1.E-4/
      LSTEP=0
      Z=0.
      ZE=1.
      S=1.
      EP=S/(1.E4*NM)
      ZEND=ZE-EP
      DO 1 I=1,N
1     SUM(I)=(0.,0.)
      NS=NX
      NT=0
      CALL SAOA (Z,G1)
2     DZ=S/NS
      IF (Z+DZ.LE.ZE) GO TO 3
      DZ=ZE-Z
      IF (DZ.LE.EP) GO TO 17
3     DZOT=DZ*.5
      CALL SAOA (Z+DZOT,G3)
      CALL SAOA (Z+DZ,G5)
4     NOGO=0
      DO 5 I=1,N
      T00=(G1(I)+G5(I))*DZOT
      T01(I)=(T00+DZ*G3(I))*.5
      T10(I)=(4.*T01(I)-T00)/3.
C     TEST CONVERGENCE OF 3 POINT ROMBERG RESULT
      CALL TEST (DREAL(T01(I)),DREAL(T10(I)),TR,DIMAG(T01(I)),DIMAG(T10
     1(I)),TI,0.d0)
      IF (TR.GT.RX.OR.TI.GT.RX) NOGO=1
5     CONTINUE
      IF (NOGO.NE.0) GO TO 7
      DO 6 I=1,N
6     SUM(I)=SUM(I)+T10(I)
      NT=NT+2
      GO TO 11
7     CALL SAOA (Z+DZ*.25,G2)
      CALL SAOA (Z+DZ*.75,G4)
      NOGO=0
      DO 8 I=1,N
      T02=(T01(I)+DZOT*(G2(I)+G4(I)))*.5
      T11=(4.*T02-T01(I))/3.
      T20(I)=(16.*T11-T10(I))/15.
C     TEST CONVERGENCE OF 5 POINT ROMBERG RESULT
      CALL TEST (DREAL(T11),DREAL(T20(I)),TR,DIMAG(T11),DIMAG(T20(I)),TI
     1,0.d0)
      IF (TR.GT.RX.OR.TI.GT.RX) NOGO=1
8     CONTINUE
      IF (NOGO.NE.0) GO TO 13
9     DO 10 I=1,N
10    SUM(I)=SUM(I)+T20(I)
      NT=NT+1
11    Z=Z+DZ
      IF (Z.GT.ZEND) GO TO 17
      DO 12 I=1,N
12    G1(I)=G5(I)
      IF (NT.LT.NTS.OR.NS.LE.NX) GO TO 2
      NS=NS/2
      NT=1
      GO TO 2
13    NT=0
      IF (NS.LT.NM) GO TO 15
      IF (LSTEP.EQ.1) GO TO 9
      LSTEP=1
      CALL LAMBDA (Z,T00,T11)
      WRITE(*,18) T00
      WRITE(*,19) Z,DZ,A,B
      DO 14 I=1,N
14    WRITE(*,19) G1(I),G2(I),G3(I),G4(I),G5(I)
      GO TO 9
15    NS=NS*2
      DZ=S/NS
      DZOT=DZ*.5
      DO 16 I=1,N
      G5(I)=G3(I)
16    G3(I)=G2(I)
      GO TO 4
17    CONTINUE
      RETURN
C
18    FORMAT (38H ROM1 -- STEP SIZE LIMITED AT LAMBDA =,1P2E12.5)
19    FORMAT (1X,1P10E12.5)
      END
      SUBROUTINE SAOA (T,ANS)
C
C     SAOA COMPUTES THE INTEGRAND FOR EACH OF THE 6
C     SOMMERFELD INTEGRALS FOR SOURCE AND OBSERVER ABOVE GROUND
C
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      COMPLEX*16 ANS,XL,DXL,CGAM1,CGAM2,B0,B0P,COM,CK1,CK1SQ,CKSM,CT1,
     1CT2,CT3,DGAM,DEN1,DEN2
      COMMON /EVLCOM/ CKSM,CT1,CT2,CT3,CK1,CK1SQ,CK2,CK2SQ,TKMAG,TSMAG,C
     1K1R,ZPH,RHO,JH
      DIMENSION ANS(6)
      CALL LAMBDA (T,XL,DXL)
      IF (JH.GT.0) GO TO 1
C     BESSEL FUNCTION FORM
      CALL BESSEL (XL*RHO,B0,B0P)
      B0=2.*B0
      B0P=2.*B0P
      CGAM1=SQRT(XL*XL-CK1SQ)
      CGAM2=SQRT(XL*XL-CK2SQ)
      IF (DREAL(CGAM1).EQ.0.) CGAM1=DCMPLX(0.D0,-ABS(DIMAG(CGAM1)))
      IF (DREAL(CGAM2).EQ.0.) CGAM2=DCMPLX(0.D0,-ABS(DIMAG(CGAM2)))
      GO TO 2
C     HANKEL FUNCTION FORM
1     CALL HANKEL (XL*RHO,B0,B0P)
      COM=XL-CK1
      CGAM1=SQRT(XL+CK1)*SQRT(COM)
      IF (DREAL(COM).LT.0..AND.DIMAG(COM).GE.0.) CGAM1=-CGAM1
      COM=XL-CK2
      CGAM2=SQRT(XL+CK2)*SQRT(COM)
      IF (DREAL(COM).LT.0..AND.DIMAG(COM).GE.0.) CGAM2=-CGAM2
2     XLR=XL*DCONJG(XL)
      IF (XLR.LT.TSMAG) GO TO 3
      IF (DIMAG(XL).LT.0.) GO TO 4
      XLR=DREAL(XL)
      IF (XLR.LT.CK2) GO TO 5
      IF (XLR.GT.CK1R) GO TO 4
3     DGAM=CGAM2-CGAM1
      GO TO 7
4     SIGN=1.
      GO TO 6
5     SIGN=-1.
6     DGAM=1./(XL*XL)
      DGAM=SIGN*((CT3*DGAM+CT2)*DGAM+CT1)/XL
7     DEN2=CKSM*DGAM/(CGAM2*(CK1SQ*CGAM2+CK2SQ*CGAM1))
      DEN1=1./(CGAM1+CGAM2)-CKSM/CGAM2
      COM=DXL*XL*EXP(-CGAM2*ZPH)
      ANS(6)=COM*B0*DEN1/CK1
      COM=COM*DEN2
      IF (RHO.EQ.0.) GO TO 8
      B0P=B0P/RHO
      ANS(1)=-COM*XL*(B0P+B0*XL)
      ANS(4)=COM*XL*B0P
      GO TO 9
8     ANS(1)=-COM*XL*XL*.5
      ANS(4)=ANS(1)
9     ANS(2)=COM*CGAM2*CGAM2*B0
      ANS(3)=-ANS(4)*CGAM2*RHO
      ANS(5)=COM*B0
      RETURN
      END
      SUBROUTINE TEST (F1R,F2R,TR,F1I,F2I,TI,DMIN)
C
C     TEST FOR CONVERGENCE IN NUMERICAL INTEGRATION
C
      IMPLICIT REAL*8(A-H,O-Z)
      SAVE
      DEN=ABS(F2R)
      TR=ABS(F2I)
      IF (DEN.LT.TR) DEN=TR
      IF (DEN.LT.DMIN) DEN=DMIN
      IF (DEN.LT.1.E-37) GO TO 1
      TR=ABS((F1R-F2R)/DEN)
      TI=ABS((F1I-F2I)/DEN)
      RETURN
1     TR=0.
      TI=0.
      RETURN
      END

      SUBROUTINE SECOND (CPUSECD)
C     Purpose:
C     SECOND returns cpu time in seconds.  Must be customized!!!
      REAL*8 CPUSECD
      integer Iticks

C--   Not customized:
C       Cpusecd = 0.0            ! if we have no clock routine
C--   MACINTOSH:
C       CPUSECD= LONG(362)/60.0
C--   Lahey fortran
C        Call Timer(Iticks)
C        cpusecd = Iticks/100.d0
      END