      SUBROUTINE ADT(X,CXREFIN,QABS,QSCA,QEXT,DQDO,
     &               THETAS,MXANG,NTHETAS)
      IMPLICIT NONE

! v2

! arguments:

      INTEGER MXANG,NTHETAS
      DOUBLE PRECISION QABS,QEXT,QSCA,X
      DOUBLE PRECISION 
     &   DQDO(0:MXANG),
     &   THETAS(0:MXANG)
      DOUBLE COMPLEX CXREFIN

! local variables:

      INTEGER NLEVMX
      PARAMETER(NLEVMX=28)
      
      INTEGER JLEV,JT,JTH,JU,NLEV,NTERM
      DOUBLE PRECISION BETA,DU,ERT,FAC,FAC2,FAC3,PI,
     &   RHO_0,RHO_1,RHO_2,U,U0,XSINTH
      DOUBLE COMPLEX CXI,CXIRHO,CXRHO
      DOUBLE COMPLEX
     &   Q(1:NLEVMX),
     &   CXSUM(0:NLEVMX)

! intrinsic functions

      INTRINSIC DIMAG

! external functions

      DOUBLE PRECISION BESJ0
      EXTERNAL BESJ0

!=======================================================================
! Subroutine ADT
! Given:
!   MXANG   = dimensioning information for arrays THETAS and DQDO
!   X       = 2*pi*a/lambda
!   CXREFIN = m = complex refractive index
!                 [we adopt convention Im(CXREFIN) > 0 ]
!   NTHETAS = number of scattering angles - 1
!             NOTE: NTHETAS must not exceed MXANG
!   THETAS(0-NTHETAS) = scattering angles (in radians)
!             NOTE: scattering angles do not need to be uniformly
!                   spaced (or monotonic).
!
! Returns:
!   QABS    = absorption cross section/pi*a^2
!   QSCA    = scattering cross section/pi*a^2
!   QEXT    = extinction cross section/pi*a^2
!   DQDO(0-NTHETAS) = dQ_sca/dOmega for angles JTH=0-NTHETAS
!
! Requires:
!   double precision external function BESJ0(X) to evaluate
!   ordinary Bessel function J_0(x) for calculation of DQDO
!   Below we include function BESJ0 from netlib SPECFUN library.
!
! Argument list: Calling program should have 
!   MXANG,NTHETAS declared as INTEGER
!   X declared as DOUBLE PRECISION
!   CXREFIN declared as DOUBLE COMPLEX
!   THETAS declared as DOUBLE PRECISION THETAS(0:MXANG)
!   DQDO declared as DOUBLE PRECISION DQDO(0:MXANG)
            
! Note: here parameter NLEVMX limits application of ADT to extremely
! large values of X.  Maximum number of step halvings is NLEVMX,
! hence maximum number of integration points is 2**(NLEVMX+1)
! User may try changing value of NLEVMX
!-----------------------------------------------------------------------
!
! Subroutine ADT uses anomalous diffraction theory, valid 
! when |m-1| << 1 and  x >> 1 
!
! Note that |m-1|*x can have an arbitrary value
!
! See Chapter 11 of van de Hulst (1957),
! Light Scattering by Small Particles
!
! B.T. Draine, Princeton University Observatory, 2002.05.07
! history
! 02.05.07 (BTD) first written
! 02.05.09 (BTD) calculation of differential scattering cross section
!                totally rewritten to increase computational
!                efficiency by storing and reusing Bessel function
!                values
! 02.07.25 (BTD) concluded that there was no real gain in computational
!                efficiency.  Abandoned strategy to reuse Bessel
!                values.  Scattering angles no longer need to be 
!                uniformly spaced.
! 02.09.06 (BTD) changed NLEVMX from 18 to 20
! 02.12.18 (BTD) changed NLEVMX from 20 to 24
!                changed ERT from 1.D-3 to 1.D-4
!                changed NLEVMX from 24 to 26
! 15.03.11 (BTD) corrected sign mistake in expression for Qext
!                (went unnoticed until today!)
! 15.06.25 (BTD) v2: removed small-angle approximation
!                replaced Theta by sin(Theta) in argument of J0
!     end history
!=======================================================================
      DATA CXI/(0.D0,1.D0)/,ERT/1.D-4/
      SAVE CXI,ERT

! Following four statements should be enabled if NOT using g77.
! They assume that the compiler supports double complex, since the
! statements DBLE and DIMAG are used.  If double complex is not available
! (see above) you will need to change DIMAG to AIMAG
!
! If using g77, following statments could be commented out, as 
! REALPART and IMAGPART are g77 intrinsic functions
! However, they do not need to be commented out.

      DOUBLE PRECISION IMAGPART,REALPART
      DOUBLE COMPLEX DCXVAR
      REALPART(DCXVAR)=(DBLE(DCXVAR))
      IMAGPART(DCXVAR)=(DIMAG(DCXVAR))
!=======================================================================      

      PI=4.D0*ATAN(1.D0)
      CXRHO=2.D0*X*(CXREFIN-1.D0)
      CXIRHO=CXI*CXRHO
      RHO_0=ABS(CXRHO)
      RHO_1=REALPART(CXRHO)
      RHO_2=IMAGPART(CXRHO)

      BETA=ACOS(RHO_1/RHO_0)

      IF(RHO_0.LT.1.D-3)THEN
         QEXT=(4.D0/3.D0)*RHO_2+0.5D0*(RHO_1**2-RHO_2**2)
         QABS=(4.D0/3.D0)*RHO_2-RHO_2**2
         QSCA=0.5D0*RHO_0**2
      ELSE
         FAC=EXP(-RHO_2)
         FAC2=FAC*FAC
!--------------------------
! 15.03.11 (BTD) correction
!         QEXT=2.D0+4.D0*(COS(2.D0*BETA)-
!     &        FAC*(COS(RHO_1-2.D0*BETA)-RHO_0*SIN(RHO_1-BETA)))/RHO_0**2
         QEXT=2.D0+4.D0*(COS(2.D0*BETA)-
     &        FAC*(COS(RHO_1-2.D0*BETA)+RHO_0*SIN(RHO_1-BETA)))/RHO_0**2
!--------------------------
         QABS=1.D0+(FAC2+0.5D0*(FAC2-1.D0)/RHO_2)/RHO_2
         QSCA=QEXT-QABS
      ENDIF

!-----------------------------------------------------------------------
! calculate differential scattering cross section

      FAC=X*X/PI

      DO JTH=0,NTHETAS

! 2015.06.25 correction:
! prior to 2015.06.25, was using small-angle approx sin(theta)=theta
         
         XSINTH=X*SIN(THETAS(JTH))

         IF(XSINTH.EQ.0.D0)THEN

! if scattering angle is zero, evaluate forward scattering analytically:

            DQDO(JTH)=FAC*ABS(0.5D0+(1.D0+(CXIRHO-1.D0)*EXP(CXIRHO))/
     &                        CXRHO**2)**2
         ELSE

! evaluate using adaptive Simpson's rule, using U=cos(tau) as integration
! variable.  Note that integrand vanishes at endpoints U=0 and U=1.

            U=0.5D0

! 2015.06.25 eq. 14 of Draine & Allaf-Akbari (2006) requires
!            two corrections:
!               Theta -> sin(Theta)
!               exp[-i*rho*sin(u)] -> exp[i*rho*sin(u)]
!            here U=cos(u) and sqrt(1-U**2)=sin(u)

            CXSUM(0)=(1.D0-EXP(CXIRHO*SQRT(1.D0-U*U)))*BESJ0(XSINTH*U)*U
            DO JLEV=1,NLEVMX
               NTERM=2**JLEV
               DU=1.D0/DBLE(NTERM)
               U0=-0.5D0*DU
               CXSUM(JLEV)=0.D0
               DO JU=1,NTERM
                  U=U0+DU*JU
                  FAC2=BESJ0(XSINTH*U)*U
                  FAC3=SQRT(1.D0-U*U)
                  CXSUM(JLEV)=CXSUM(JLEV)+(1.D0-EXP(CXIRHO*FAC3))*FAC2
               ENDDO

! quadrature using Simpson's rule:

               Q(JLEV)=2.D0*CXSUM(JLEV)
               DO JT=0,JLEV-1
                  Q(JLEV)=Q(JLEV)+CXSUM(JT)
               ENDDO
               Q(JLEV)=Q(JLEV)*DU/3.D0
               NLEV=JLEV

! for NLEV>1, check to see whether Simpson's rule integration has converged

               IF(NLEV.GT.1)THEN
                  IF(ABS(Q(NLEV)-Q(NLEV-1))/ABS(Q(NLEV)+Q(NLEV-1)).LT.
     &               ERT)GOTO 3500
               ENDIF
            ENDDO

! if reach here, we failed to meet error tolerance.
! write out warning, but use best estimate for diff. scatt. cross section

            WRITE(0,9200)NLEVMX
 3500       DQDO(JTH)=FAC*ABS(Q(NLEV))**2
         ENDIF
      ENDDO
      
      RETURN

 9200 FORMAT('Warning: adt failed to meet error tolerance in NLEV=',I2,
     &          ' doublings')
      END

      DOUBLE PRECISION FUNCTION BESJ0(X)
      IMPLICIT NONE
!--------------------------------------------------------------------
!
! This subprogram computes approximate values for Bessel functions
!   of the first kind of order zero for arguments  |X| <= XMAX
!   (see comments heading CALJY0).
!
! Obtained from www.netlib.org/specfun/
! specfun library (special functions and accompanying test programs
!                  written by W.J. Cody, Argonne National Laboratory)
!--------------------------------------------------------------------
      INTEGER JINT
      DOUBLE PRECISION X,RESULT
!--------------------------------------------------------------------
      JINT=0
      CALL CALJY0(X,RESULT,JINT)
      BESJ0 = RESULT
      RETURN
!---------- Last line of BESJ0 ----------
      END
      DOUBLE PRECISION FUNCTION BESY0(X)
      IMPLICIT NONE
!--------------------------------------------------------------------
!
! This subprogram computes approximate values for Bessel functions
!   of the second kind of order zero for arguments 0 < X <= XMAX
!   (see comments heading CALJY0).
!
!--------------------------------------------------------------------
      INTEGER JINT
      DOUBLE PRECISION X,RESULT
!--------------------------------------------------------------------
      JINT=1
      CALL CALJY0(X,RESULT,JINT)
      BESY0 = RESULT
      RETURN
!---------- Last line of BESY0 ----------
      END
      SUBROUTINE CALJY0(ARG,RESULT,JINT)
!---------------------------------------------------------------------
!
! This packet computes zero-order Bessel functions of the first and
!   second kind (J0 and Y0), for real arguments X, where 0 < X <= XMAX
!   for Y0, and |X| <= XMAX for J0.  It contains two function-type
!   subprograms,  BESJ0  and  BESY0,  and one subroutine-type
!   subprogram,  CALJY0.  The calling statements for the primary
!   entries are:
!
!           Y = BESJ0(X)
!   and
!           Y = BESY0(X),
!
!   where the entry points correspond to the functions J0(X) and Y0(X),
!   respectively.  The routine  CALJY0  is intended for internal packet
!   use only, all computations within the packet being concentrated in
!   this one routine.  The function subprograms invoke  CALJY0  with
!   the statement
!           CALL CALJY0(ARG,RESULT,JINT),
!   where the parameter usage is as follows:
!
!      Function                  Parameters for CALJY0
!       call              ARG             RESULT          JINT
!
!     BESJ0(ARG)     |ARG| .LE. XMAX       J0(ARG)          0
!     BESY0(ARG)   0 .LT. ARG .LE. XMAX    Y0(ARG)          1
!
!   The main computation uses unpublished minimax rational
!   approximations for X .LE. 8.0, and an approximation from the 
!   book  Computer Approximations  by Hart, et. al., Wiley and Sons, 
!   New York, 1968, for arguments larger than 8.0   Part of this
!   transportable packet is patterned after the machine-dependent
!   FUNPACK program BESJ0(X), but cannot match that version for
!   efficiency or accuracy.  This version uses rational functions
!   that are theoretically accurate to at least 18 significant decimal
!   digits for X <= 8, and at least 18 decimal places for X > 8.  The
!   accuracy achieved depends on the arithmetic system, the compiler,
!   the intrinsic functions, and proper selection of the machine-
!   dependent constants.
!
!*******************************************************************
!
! Explanation of machine-dependent constants
!
!   XINF   = largest positive machine number
!   XMAX   = largest acceptable argument.  The functions AINT, SIN
!            and COS must perform properly for  ABS(X) .LE. XMAX.
!            We recommend that XMAX be a small integer multiple of
!            sqrt(1/eps), where eps is the smallest positive number
!            such that  1+eps > 1. 
!   XSMALL = positive argument such that  1.0-(X/2)**2 = 1.0
!            to machine precision for all  ABS(X) .LE. XSMALL.
!            We recommend that  XSMALL < sqrt(eps)/beta, where beta
!            is the floating-point radix (usually 2 or 16).
!
!     Approximate values for some important machines are
!
!                          eps      XMAX     XSMALL      XINF  
!
!  CDC 7600      (S.P.)  7.11E-15  1.34E+08  2.98E-08  1.26E+322
!  CRAY-1        (S.P.)  7.11E-15  1.34E+08  2.98E-08  5.45E+2465
!  IBM PC (8087) (S.P.)  5.96E-08  8.19E+03  1.22E-04  3.40E+38
!  IBM PC (8087) (D.P.)  1.11D-16  2.68D+08  3.72D-09  1.79D+308
!  IBM 195       (D.P.)  2.22D-16  6.87D+09  9.09D-13  7.23D+75
!  UNIVAC 1108   (D.P.)  1.73D-18  4.30D+09  2.33D-10  8.98D+307
!  VAX 11/780    (D.P.)  1.39D-17  1.07D+09  9.31D-10  1.70D+38
!
!*******************************************************************
!*******************************************************************
!
! Error Returns
!
!  The program returns the value zero for  X .GT. XMAX, and returns
!    -XINF when BESLY0 is called with a negative or zero argument.
!
!
! Intrinsic functions required are:
!
!     ABS, AINT, COS, LOG, SIN, SQRT
!
!
!  Latest modification: June 2, 1989
!
!  Author: W. J. Cody
!          Mathematics and Computer Science Division 
!          Argonne National Laboratory
!          Argonne, IL 60439
!
!--------------------------------------------------------------------
      INTEGER I,JINT
      DOUBLE PRECISION
     1       ARG,AX,CONS,DOWN,EIGHT,FIVE5,FOUR,ONE,ONEOV8,PI2,PJ0,
     2       PJ1,PLG,PROD,PY0,PY1,PY2,P0,P1,P17,QJ0,QJ1,QLG,QY0,QY1,
     3       QY2,Q0,Q1,RESJ,RESULT,R0,R1,SIXTY4,THREE,TWOPI,TWOPI1,
     4       TWOPI2,TWO56,UP,W,WSQ,XDEN,XINF,XMAX,XNUM,XSMALL,XJ0,
     5       XJ1,XJ01,XJ02,XJ11,XJ12,XY,XY0,XY01,XY02,XY1,XY11,XY12,
     6       XY2,XY21,XY22,Z,ZERO,ZSQ
      DIMENSION PJ0(7),PJ1(8),PLG(4),PY0(6),PY1(7),PY2(8),P0(6),P1(6),
     1          QJ0(5),QJ1(7),QLG(4),QY0(5),QY1(6),QY2(7),Q0(5),Q1(5)
!-------------------------------------------------------------------
!  Mathematical constants
!    CONS = ln(.5) + Euler's gamma
!-------------------------------------------------------------------
      DATA ZERO,ONE,THREE,FOUR,EIGHT/0.0D0,1.0D0,3.0D0,4.0D0,8.0D0/,
     1     FIVE5,SIXTY4,ONEOV8,P17/5.5D0,64.0D0,0.125D0,1.716D-1/,
     2     TWO56,CONS/256.0D0,-1.1593151565841244881D-1/,
     3     PI2,TWOPI/6.3661977236758134308D-1,6.2831853071795864769D0/,
     4     TWOPI1,TWOPI2/6.28125D0,1.9353071795864769253D-3/
!-------------------------------------------------------------------
!  Machine-dependent constants
!-------------------------------------------------------------------
      DATA XMAX/1.07D+09/,XSMALL/9.31D-10/,XINF/1.7D+38/
!-------------------------------------------------------------------
!  Zeroes of Bessel functions
!-------------------------------------------------------------------
      DATA XJ0/2.4048255576957727686D+0/,XJ1/5.5200781102863106496D+0/,
     1     XY0/8.9357696627916752158D-1/,XY1/3.9576784193148578684D+0/,
     2     XY2/7.0860510603017726976D+0/,
     3     XJ01/ 616.0D+0/, XJ02/-1.4244423042272313784D-03/,
     4     XJ11/1413.0D+0/, XJ12/ 5.4686028631064959660D-04/,
     5     XY01/ 228.0D+0/, XY02/ 2.9519662791675215849D-03/,
     6     XY11/1013.0D+0/, XY12/ 6.4716931485786837568D-04/,
     7     XY21/1814.0D+0/, XY22/ 1.1356030177269762362D-04/
!-------------------------------------------------------------------
!  Coefficients for rational approximation to ln(x/a)
!--------------------------------------------------------------------
      DATA PLG/-2.4562334077563243311D+01,2.3642701335621505212D+02,
     1         -5.4989956895857911039D+02,3.5687548468071500413D+02/
      DATA QLG/-3.5553900764052419184D+01,1.9400230218539473193D+02,
     1         -3.3442903192607538956D+02,1.7843774234035750207D+02/
!-------------------------------------------------------------------
!  Coefficients for rational approximation of
!  J0(X) / (X**2 - XJ0**2),  XSMALL  <  |X|  <=  4.0
!--------------------------------------------------------------------
      DATA PJ0/6.6302997904833794242D+06,-6.2140700423540120665D+08,
     1         2.7282507878605942706D+10,-4.1298668500990866786D+11,
     2        -1.2117036164593528341D-01, 1.0344222815443188943D+02,
     3        -3.6629814655107086448D+04/
      DATA QJ0/4.5612696224219938200D+05, 1.3985097372263433271D+08,
     1         2.6328198300859648632D+10, 2.3883787996332290397D+12,
     2         9.3614022392337710626D+02/
!-------------------------------------------------------------------
!  Coefficients for rational approximation of
!  J0(X) / (X**2 - XJ1**2),  4.0  <  |X|  <=  8.0
!-------------------------------------------------------------------
      DATA PJ1/4.4176707025325087628D+03, 1.1725046279757103576D+04,
     1         1.0341910641583726701D+04,-7.2879702464464618998D+03,
     2        -1.2254078161378989535D+04,-1.8319397969392084011D+03,
     3         4.8591703355916499363D+01, 7.4321196680624245801D+02/
      DATA QJ1/3.3307310774649071172D+02,-2.9458766545509337327D+03,
     1         1.8680990008359188352D+04,-8.4055062591169562211D+04,
     2         2.4599102262586308984D+05,-3.5783478026152301072D+05,
     3        -2.5258076240801555057D+01/
!-------------------------------------------------------------------
!  Coefficients for rational approximation of
!    (Y0(X) - 2 LN(X/XY0) J0(X)) / (X**2 - XY0**2),
!        XSMALL  <  |X|  <=  3.0
!--------------------------------------------------------------------
      DATA PY0/1.0102532948020907590D+04,-2.1287548474401797963D+06,
     1         2.0422274357376619816D+08,-8.3716255451260504098D+09,
     2         1.0723538782003176831D+11,-1.8402381979244993524D+01/
      DATA QY0/6.6475986689240190091D+02, 2.3889393209447253406D+05,
     1         5.5662956624278251596D+07, 8.1617187777290363573D+09,
     2         5.8873865738997033405D+11/
!-------------------------------------------------------------------
!  Coefficients for rational approximation of
!    (Y0(X) - 2 LN(X/XY1) J0(X)) / (X**2 - XY1**2),
!        3.0  <  |X|  <=  5.5
!--------------------------------------------------------------------
      DATA PY1/-1.4566865832663635920D+04, 4.6905288611678631510D+06,
     1         -6.9590439394619619534D+08, 4.3600098638603061642D+10,
     2         -5.5107435206722644429D+11,-2.2213976967566192242D+13,
     3          1.7427031242901594547D+01/
      DATA QY1/ 8.3030857612070288823D+02, 4.0669982352539552018D+05,
     1          1.3960202770986831075D+08, 3.4015103849971240096D+10,
     2          5.4266824419412347550D+12, 4.3386146580707264428D+14/
!-------------------------------------------------------------------
!  Coefficients for rational approximation of
!    (Y0(X) - 2 LN(X/XY2) J0(X)) / (X**2 - XY2**2),
!        5.5  <  |X|  <=  8.0
!--------------------------------------------------------------------
      DATA PY2/ 2.1363534169313901632D+04,-1.0085539923498211426D+07,
     1          2.1958827170518100757D+09,-1.9363051266772083678D+11,
     2         -1.2829912364088687306D+11, 6.7016641869173237784D+14,
     3         -8.0728726905150210443D+15,-1.7439661319197499338D+01/
      DATA QY2/ 8.7903362168128450017D+02, 5.3924739209768057030D+05,
     1          2.4727219475672302327D+08, 8.6926121104209825246D+10,
     2          2.2598377924042897629D+13, 3.9272425569640309819D+15,
     3          3.4563724628846457519D+17/
!-------------------------------------------------------------------
!  Coefficients for Hart,s approximation,  |X| > 8.0
!-------------------------------------------------------------------
      DATA P0/3.4806486443249270347D+03, 2.1170523380864944322D+04,
     1        4.1345386639580765797D+04, 2.2779090197304684302D+04,
     2        8.8961548424210455236D-01, 1.5376201909008354296D+02/
      DATA Q0/3.5028735138235608207D+03, 2.1215350561880115730D+04,
     1        4.1370412495510416640D+04, 2.2779090197304684318D+04,
     2        1.5711159858080893649D+02/
      DATA P1/-2.2300261666214198472D+01,-1.1183429920482737611D+02,
     1        -1.8591953644342993800D+02,-8.9226600200800094098D+01,
     2        -8.8033303048680751817D-03,-1.2441026745835638459D+00/
      DATA Q1/1.4887231232283756582D+03, 7.2642780169211018836D+03,
     1        1.1951131543434613647D+04, 5.7105024128512061905D+03,
     2        9.0593769594993125859D+01/
!-------------------------------------------------------------------
!  Check for error conditions
!-------------------------------------------------------------------
      AX = ABS(ARG)
      IF ((JINT .EQ. 1) .AND. (ARG .LE. ZERO)) THEN
            RESULT = -XINF
            GO TO 2000
         ELSE IF (AX .GT. XMAX) THEN
            RESULT = ZERO
            GO TO 2000
      END IF
      IF (AX .GT. EIGHT) GO TO 800
      IF (AX .LE. XSMALL) THEN
         IF (JINT .EQ. 0) THEN
               RESULT = ONE
            ELSE
               RESULT = PI2 * (LOG(AX) + CONS)
         END IF
         GO TO 2000
      END IF
!-------------------------------------------------------------------
!  Calculate J0 for appropriate interval, preserving
!     accuracy near the zero of J0
!-------------------------------------------------------------------
      ZSQ = AX * AX
      IF (AX .LE. FOUR) THEN
            XNUM = (PJ0(5) * ZSQ + PJ0(6)) * ZSQ + PJ0(7)
            XDEN = ZSQ + QJ0(5)
            DO 50 I = 1, 4
               XNUM = XNUM * ZSQ + PJ0(I)
               XDEN = XDEN * ZSQ + QJ0(I)
   50       CONTINUE
            PROD = ((AX - XJ01/TWO56) - XJ02) * (AX + XJ0)
         ELSE
            WSQ = ONE - ZSQ / SIXTY4
            XNUM = PJ1(7) * WSQ + PJ1(8)
            XDEN = WSQ + QJ1(7)
            DO 220 I = 1, 6
               XNUM = XNUM * WSQ + PJ1(I)
               XDEN = XDEN * WSQ + QJ1(I)
  220       CONTINUE
            PROD = (AX + XJ1) * ((AX - XJ11/TWO56) - XJ12)
      END IF
      RESULT = PROD * XNUM / XDEN
      IF (JINT .EQ. 0) GO TO 2000
!-------------------------------------------------------------------
!  Calculate Y0.  First find  RESJ = pi/2 ln(x/xn) J0(x),
!    where xn is a zero of Y0
!-------------------------------------------------------------------
      IF (AX .LE. THREE) THEN
            UP = (AX-XY01/TWO56)-XY02
            XY = XY0
         ELSE IF (AX .LE. FIVE5) THEN
            UP = (AX-XY11/TWO56)-XY12
            XY = XY1
         ELSE
            UP = (AX-XY21/TWO56)-XY22
            XY = XY2
      END IF
      DOWN = AX + XY
      IF (ABS(UP) .LT. P17*DOWN) THEN
            W = UP/DOWN
            WSQ = W*W
            XNUM = PLG(1)
            XDEN = WSQ + QLG(1)
            DO 320 I = 2, 4
               XNUM = XNUM*WSQ + PLG(I)
               XDEN = XDEN*WSQ + QLG(I)
  320       CONTINUE
            RESJ = PI2 * RESULT * W * XNUM/XDEN
         ELSE
            RESJ = PI2 * RESULT * LOG(AX/XY)
      END IF
!-------------------------------------------------------------------
!  Now calculate Y0 for appropriate interval, preserving
!     accuracy near the zero of Y0
!-------------------------------------------------------------------
      IF (AX .LE. THREE) THEN
            XNUM = PY0(6) * ZSQ + PY0(1)
            XDEN = ZSQ + QY0(1)
            DO 340 I = 2, 5
               XNUM = XNUM * ZSQ + PY0(I)
               XDEN = XDEN * ZSQ + QY0(I)
  340       CONTINUE
         ELSE IF (AX .LE. FIVE5) THEN
            XNUM = PY1(7) * ZSQ + PY1(1)
            XDEN = ZSQ + QY1(1)
            DO 360 I = 2, 6
               XNUM = XNUM * ZSQ + PY1(I)
               XDEN = XDEN * ZSQ + QY1(I)
  360       CONTINUE
         ELSE
            XNUM = PY2(8) * ZSQ + PY2(1)
            XDEN = ZSQ + QY2(1)
            DO 380 I = 2, 7
               XNUM = XNUM * ZSQ + PY2(I)
               XDEN = XDEN * ZSQ + QY2(I)
  380       CONTINUE
      END IF
      RESULT = RESJ + UP * DOWN * XNUM / XDEN
      GO TO 2000
!-------------------------------------------------------------------
!  Calculate J0 or Y0 for |ARG|  >  8.0
!-------------------------------------------------------------------
  800 Z = EIGHT / AX
      W = AX / TWOPI
      W = AINT(W) + ONEOV8
      W = (AX - W * TWOPI1) - W * TWOPI2
      ZSQ = Z * Z
      XNUM = P0(5) * ZSQ + P0(6)
      XDEN = ZSQ + Q0(5)
      UP = P1(5) * ZSQ + P1(6)
      DOWN = ZSQ + Q1(5)
      DO 850 I = 1, 4
         XNUM = XNUM * ZSQ + P0(I)
         XDEN = XDEN * ZSQ + Q0(I)
         UP = UP * ZSQ + P1(I)
         DOWN = DOWN * ZSQ + Q1(I)
  850 CONTINUE
      R0 = XNUM / XDEN
      R1 = UP / DOWN
      IF (JINT .EQ. 0) THEN
            RESULT = SQRT(PI2/AX) * (R0*COS(W) - Z*R1*SIN(W))
         ELSE
            RESULT = SQRT(PI2/AX) * (R0*SIN(W) + Z*R1*COS(W))
      END IF
 2000 RETURN
!---------- Last line of CALJY0 ----------
      END



