PROGRAM w2 CALC
C This program calculates the quantity W2 defined in equation 8.22 of
C Jeffrey and Onishi 1984. Starting with the new versions of
C equations 3.20, 3.21, we combine XA(11)+XA(12) and remove the cancelling
C singular terms algebraically (this does not improve the convergence
C of the series which now becomes a conditionally convergent
C alternating series) and then differentiate the resulting expression.
C The terms in the series go down roughly like m^(-2) so the rate of
C convergence should be like (maximum)^(-1). So 2 significant digits
C should be all we can hope for. The values listed in Table 3 of Jeffrey &
C Onishi are likely to be more accurate than those coming from this series,
C because they were calculated using data for the XA functions
C said to be accurate to 6 figures.
C
C This program can be used to obtain new values of W2 by editing the
C  DATA statements so that they contain the desired values of LAMBDA and S.
C  The output is at present set to 'standard output' called * in FORTRAN 77.
C  This can be changed to a unit number for output to a file or printer.
C
      DOUBLE PRECISION DLAMDA(5),S(4),W2,W2DAT
      INTEGER I,J
      DATA DLAMDA/1D0,.5D0,.25D0,.125D0,2D0/
      DATA S/2D0,2.01D0,2.05D0,2.1D0/
      DO 50 J=1,5
        WRITE(*,1000)DLAMDA(J)
 1000   FORMAT(/9H LAMBDA =,F8.3)
        WRITE(*,1100)
 1100   FORMAT(/10X,1HS,10X,' W2')
        DO 50 I=1,4
          W2DAT=W2(S(I),DLAMDA(J))
   50 WRITE(*,1200)S(I),W2DAT
 1200 FORMAT(F13.3,6F16.5)
      STOP
      END
C
      DOUBLE PRECISION FUNCTION W2(S,DLAMDA)
C ----------------------- Arguments of function -------------------
      DOUBLE PRECISION S, DLAMDA
C -----------------------------------------------------------------
C This subprogram calculates:
C
C      d    A                    A               A
C      - ( X  +  ( 1+lambda ) * X   +  lambda * X   )
C      ds   11                   12              22
C
C ------------------------ Local variables ----------------------------
      DOUBLE PRECISION RXACOF(0:45450), G1, G2, G3, G1SUM, G2SUM, G3SUM
      DOUBLE PRECISION ZL, COEF, XL, XLP1, XL13,T,SP2,ZT,DXADS,SIGN
      INTEGER I, ISTEP, IQ, MMIN, MMAX, M, NM, MAXS
      LOGICAL LOADED, GREAT1
      CHARACTER*4 ID
      DATA LOADED/.FALSE./
C -------------------- Start of executable code -------------------
      IF(.NOT.LOADED) THEN
        OPEN(1,FILE='rxa300.dat',STATUS='OLD')
        READ(1,10)  ID,MAXS
10      FORMAT(A4,I5)
        IF(ID.NE.' RXA') THEN
          WRITE(*,20)
20        FORMAT(1X,'WRONG COEFFICIENTS FOR RXA')
          STOP
        ENDIF
        READ(1,30) RXACOF
30      FORMAT(D22.16)
        CLOSE(1)
        LOADED=.TRUE.
      ENDIF
      IF(DLAMDA.LE.1D0) THEN
        GREAT1=.FALSE.
        XL=DLAMDA
      ELSE
        GREAT1=.TRUE.
        XL=1D0/DLAMDA
      ENDIF
C
C  First I calculate:
C      d      A    1+lambda     A
C      -  (  X  +  --------  * X   )
C      ds     11       2        12
C
C  Then I calculate:
C      d      A    1+lambda    A
C      -  (  X  +  -------- * X    )
C      ds     22   2*lambda    12
C
      W2=0D0
      DO 200 I=1,2
        XLP1=XL+1D0
        XL13 = XLP1**3
C  These are the g_i for X^A_11 and X^A_12 (i=1..3 in paper)
        IF(I.EQ.1) THEN
          G1= 2D0*XL**2 / XL13
          G2= (XL+7D0*XL**2+XL**3) / (5D0*XL13)
          G3= (1D0+18D0*XL-29D0*XL**2+18D0*XL**3+XL**4)/(42D0*XL13)
          ISTEP=1
        ELSE
          G1= 2D0*XL / XL13
          G2= (1D0+7D0*XL+XL**2) / (5D0*XL13)
          G3= (1D0+18D0*XL-29D0*XL**2+18D0*XL**3+XL**4)/(42D0*XL13*XL)
          ISTEP=-1
        ENDIF
        T=2D0/S
        SP2=S+2D0
        ZT=2D0/(S*S)
        DXADS= 2D0*G1/(SP2*SP2) + 4D0*G2/(SP2*S)
        DXADS=DXADS - 2D0*G3/S - G3*S*DLOG(SP2/S)+2D0*G3
        MMIN=1
        MMAX=MAXS
        G1SUM=G1
        G2SUM=2D0*G2
        G3SUM=4D0*G3
        SIGN= 1D0
        DO 100 M=MMIN,MMAX
          IF (ISTEP .GT. 0) THEN
            NM= ( M*(M+1) )/2
          ELSE
            NM =  M*(M+1)/2 + M
          ENDIF
          ZL=1D0
          COEF=0D0
          DO 50 IQ=0,M
            COEF=COEF+RXACOF(NM)*ZL
            NM=NM+ISTEP
            ZL=ZL*XL
50        CONTINUE
          COEF=COEF*XLP1**(-M)
          DXADS=DXADS +SIGN*ZT*
     &   ( DBLE(M)*COEF-(DBLE(M)*G1SUM+ G2SUM-G3SUM/DBLE(M+2) ) )
          SIGN=-SIGN
          ZT=ZT*T
100     CONTINUE
        IF(I.EQ.2)DXADS=DXADS*XL
        W2=W2+DXADS
200   CONTINUE
C  Unlike other functions, W2(lambda)=lambda*W2(1/lambda).
      IF(GREAT1) W2=W2*DLAMDA
      RETURN
      END