PROGRAM axitr3
C  Comment lines found in axitr1 are not repeated here, and should be read
C  first. 

C  ------------------------ WARNING ---------------------------
C  This program uses quadruple precision (REAL*16) which is only supported
C  by a few compilers.  It was the program actually used to obtain the 
C  data files.
C

C  The dimensions given below are for MAXS = 50.
      REAL*16 P(0:2307277),V(0:2307277),FAC0(151,151),XN,XS
      REAL*16 FCTR1,FCTR2,FCTR3,FCTRV,FAC,PSUM,VSUM
      INTEGER MMAX,MAXS,NM,N,IS,M,NQ,IQ,INDX1,JS,IP
      COMMON /ARRAYS/P,V,FAC0
C  This COMMON block reduces disk space requirements with many compilers
C   and helps memory management on some computers.
      MAXS=300
      MMAX=MAXS+1
C  First tabulate FAC0
      DO 10 N=1,(MMAX+1)/2
        XN=qext(N)
        FAC0(N,1)=1q0+XN
        DO 10 IS=2, (MMAX+1)/2
          XS=qext(IS)
          FAC0(N,IS)=FAC0(N,IS-1)*(XN+XS)/XS
10    CONTINUE
C  We set the initial conditions.
      P(INDX1(1,0,0))=1q0
      V(INDX1(1,0,0))=1q0
C  We start at N=1, Q=0, P=M and proceed keeping M fixed.
C   NM is the solution of the equation M=N+P+Q-1=NM+NM-1+0-1, except
C   the last time through when only N=1 and 2 are calculated.
      DO 500 M=1,MMAX
       Write(*,*)m
       NM=M/2+1
       IF (M.EQ.MMAX)NM=2
       DO 400 N=1,NM
C  The equation for NQ is M=N+P+Q-1=N+N-1+NQ-1
        NQ=M-2*N+2
        XN=qext(N)
        FCTR1=XN*(XN+.5q0)/(XN+1q0)
        FCTR2=XN*(XN-.5q0)/(XN+1q0)
        FCTR3=XN*(2q0*XN*XN-.5q0)/(XN+1q0)
        FCTRV=2q0*XN/((XN+1q0)*(2q0*XN+3q0))
        DO 300 IQ=0,NQ
C  Now that M, N, Q are specified, I know P (denoted IP) as well.
C   We obtain JS from the equation JS-1=IQ-JS.
         IP=M-IQ-N+1
         PSUM=0q0
         VSUM=0q0
         IF (IQ.EQ.0q0) GO TO 250
         JS=(IQ+1)/2
         DO 200 IS=1,JS
          XS=qext(IS)
          FAC=FAC0(N,IS)
          PSUM=PSUM+FAC*FCTR1*(2q0*XN*XS-XN-XS+2q0)*
     &      P(INDX1(IS,IQ-IS,IP-N+1))/((XN+XS)*(2q0*XS-1q0))
          IF(IP.GT.N) THEN
            PSUM=PSUM-FAC*FCTR2*P(INDX1(IS,IQ-IS,IP-N-1))
            VSUM=VSUM-FAC*FCTRV*P(INDX1(IS,IQ-IS,IP-N-1))
          ENDIF
          IF(IS.LT.JS) THEN
            PSUM=PSUM-FAC*FCTR3*V(INDX1(IS,IQ-IS-2,IP-N+1))/(2q0*XS+1q0)
          ENDIF
200      CONTINUE
250      P(INDX1(N,IP,IQ))=PSUM
         V(INDX1(N,IP,IQ))=PSUM+VSUM
300     CONTINUE
400    CONTINUE
500   CONTINUE
C  Collect the coefficients for writing.
C
      DO 600 M=0,MAXS
        NM=(M*(M+1))/2
        DO 600 IQ=0,M
          V(NM+IQ)=P(INDX1(1,M-IQ,IQ))
600   CONTINUE
      P(0)=0q0
      DO 620 M=1,MAXS
        NM = (M*(M+1))/2
        DO 610 IQ=0, M-1
610     P(NM+IQ)=0.75q0*P( INDX1(2,M-IQ,IQ) )
        P(NM+M) =0q0
620   CONTINUE
C
      MMAX=((MAXS+1)*(MAXS+2))/2 -1
      OPEN (1,FILE='rxacof.dat')
C First we write an identification line
      WRITE(1,650) MAXS
650   FORMAT(' RXA',I5)   
C We do not use an implied DO loop because one computer (ETA-10) crashed.
C  The FORMAT may need a leading space on some compilers.
      DO 690 M=0,MMAX
690   WRITE (1,700)V(M)
700   FORMAT(D22.16)
      CLOSE(1)
C
      OPEN (2,FILE='rxgcof.dat')
      WRITE(2,710) MAXS
710   FORMAT(' RXG',I5)
      DO 720 M=0,MMAX
720   WRITE (2,700)P(M)
      CLOSE (2)
C Previously the equal spheres case was written out separately using
C the commented out code below.

C      OPEN (3,FILE='rxaeqcof.dat')
C      DO 940 M=0,MAXS
C        VSUM=0q0
C        NM=((M+1)*M)/2
C        DO 930 IQ=0,M
C          VSUM=VSUM+V(NM+IQ)
C930     CONTINUE
C        VSUM=VSUM*2q0**(-M)
C        WRITE(3,700)VSUM
C940   CONTINUE
C      CLOSE(3)
C      OPEN (4,FILE='rxgeqcof.dat')
C      DO 990 M=0,MAXS
C        PSUM=0q0
C        NM=((M+1)*M)/2
C        DO 970 IQ=0,M
C          PSUM=PSUM+P(NM+IQ)
C970     CONTINUE
C        PSUM=PSUM*2q0**(-M)
C        WRITE(4,700)PSUM
C990   CONTINUE
C      CLOSE(4)
      STOP
      END
C
      INTEGER FUNCTION INDX1(N,IP,IQ)
C ------------------------------ Arguments -------------------------
      INTEGER N,IP,IQ
C --------------------------- Local Variables ----------------------
      INTEGER M,M2,LM
      M=IP+IQ+N-1
      M2=M/2
      LM=( M2*(M2+1) )/2
      LM=LM*(M2+1)+ ( LM*(M2+2) )/3
C LM is calculated as above to avoid overflow.
C If MMAX=200, then LM=681750
      LM=LM+(M-2*M2)*(M2+1)*(M2+1)
      INDX1=LM +(N-1)*(M-N+3)+IQ
      RETURN
      END