PROGRAM axist4
C Comment lines found in axist1 are not repeated here, and should be read
C  first. This program modifies axist1 only in the output of the
C  coefficients for the resistance functions RXM and RXG. The results
C  are written out in form suitable for the subroutines RXG, RGX, RXM
C  and RMX.
C  To calculate large numbers of terms, the memory requirements are:
C
C   MAXS    MMAX     DIMENSION P OR V    MEMORY (PER 8 BYTES)     C0
C   ----    ----     -----------------   ------ (for each array)  --
C     15      17           558              5KB                   9,9
C     50      52         12856            100KB                  26,26
C    100     102         92581            740KB                  51,51
C    200     202        702656            5.6MB                 101,101
C    300     302       2330231           18.6MB                 151,151
C
C  The dimensions given below are for MAXS=50.
      DOUBLE PRECISION P(0:2330231),V(0:2330231),C0(151,151),XN,XS
      DOUBLE PRECISION FCTR1,FCTR2,FCTR3,FCTRV,FAC,PSUM,VSUM
      INTEGER MMAX,MAXS,NM,N,IS,M,NQ,IQ,INDX2,JS,IP
      LOGICAL EVEN
      COMMON/ ARRAYS/P,V,C0
C  This COMMON block reduces disk space requirements with many compilers
C   and helps memory management on some computers.
      MAXS=300
      MMAX=MAXS+2
C  First tabulate C0
      DO 10 N=1,(MMAX+1)/2
        XN=DBLE(N)
        C0(N,1)=1D0
        DO 10 IS=2,(MMAX+1)/2
          XS=DBLE(IS)
          C0(N,IS)=C0(N,IS-1)*(XN+XS)/XS
10    CONTINUE
C We set the initial conditions. The labelling scheme needs some zero
C  elements at the start.
      DO 5 M=0,5
        P(M)=0D0
5     V(M)=0D0
      P(INDX2(2,0,0))=1D0
      V(INDX2(2,0,0))=1D0
C We start at N=1, Q=0, P=M-1 and proceed keeping M fixed.
C  NM is the solution of the equation M=N+P+Q=NM+NM-2+0, except
C  the last time through when only N=1 and 2 are calculated.
      DO 500 M=3,MMAX
        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=N+N-2+NQ
          NQ=M-2*(N-1)
          XN=DBLE(N)
          FCTR1=XN*(XN+.5D0)
          FCTR2=XN*(XN-.5D0)
          FCTR3=XN*(2D0*XN*XN-.5D0)
          FCTRV=2D0*XN/(2D0*XN+3D0)
          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-2=IQ-JS.
            IP=M-IQ-N
            PSUM=0D0
            VSUM=0D0
C
            IF ( IQ.EQ.0 .OR. IQ.EQ.NQ) GO TO 250
            JS=IQ/2+1
            DO 200 IS=1,JS
              XS=DBLE(IS)
              FAC=C0(N,IS)
C
              PSUM=PSUM+FAC*FCTR1*(2D0*XN*XS-XN-XS+2D0)*
     &         P(INDX2(IS,IQ-IS,IP-N+1))/((XN+XS)*(2D0*XS-1D0))
              IF(IP.GT.N) THEN
                PSUM=PSUM-FAC*FCTR2*P(INDX2(IS,IQ-IS,IP-N-1))
                VSUM=VSUM-FAC*FCTRV*P(INDX2(IS,IQ-IS,IP-N-1))
              ENDIF
C
              IF(IS.LT.JS) THEN
                PSUM=PSUM-FAC*FCTR3*V(INDX2(IS,IQ-IS-2,IP-N+1))/
     &          (2D0*XS+1D0)
              ENDIF
  200       CONTINUE
  250       P(INDX2(N,IP,IQ))=PSUM
            V(INDX2(N,IP,IQ))=PSUM+VSUM
  300     CONTINUE
  400   CONTINUE
  500 CONTINUE
C  Collect the coefficients for writing.

      do 690 m=1,maxs
        EVEN= M .EQ. 2*(M/2)
        NM = (M*(M+1))/2
        v(nm)=0d0
        if (even) then
          do 640 iq=1,m-1
          ip=m-iq-1
          psum=0d0
          do 630 n=1,(iq+2)/2
            psum=psum+p( indx2(n,iq-n,ip) )
630         continue
640       v( nm+iq) = 2.5d0*psum
         v(nm+m)=0d0
        else
          do 660 iq=1,m
          ip=m-iq
          psum=0d0
          do 650 n=1,(iq+2)/2
            psum=psum+p( indx2(n,iq-n-1,ip) )
650         continue
660       v( nm+iq) = 2.5d0*psum
        endif
690   continue
C
      V(0)=-13D0
      MMAX=((MAXS+1)*(MAXS+2))/2 -1
      OPEN (1,FILE='rxqcof.dat',STATUS='NEW')
      WRITE(1,710) MAXS
 710  FORMAT(' RXQ',I5)
C  We do not use an implied DO loop because the computer (ETA-10) crashes.
C   The FORMAT may need a leading space on some compilers.
      DO 790 M=0,MMAX
790   WRITE (1,700)V(M)
700   FORMAT(D22.16)
      CLOSE (1)
C  Previously the equal spheres case was written out separatly using 
C  the commented out code below.
C900   FORMAT(D27.21)
C      OPEN (3,FILE='rxqeqcof.dat')
C      DO 940 M=0,MAXS
C        VSUM=0D0
C        NM=((M+1)*M)/2
C        DO 930 IQ=0,M
C          VSUM=VSUM+V(NM+IQ)
C930     CONTINUE
C        VSUM=VSUM*2D0**(-M)
C        WRITE(3,900)VSUM
C940   CONTINUE
C      CLOSE (3)
C
      STOP
      END
C
      INTEGER FUNCTION INDX2(N,IP,IQ)
C --------------------------- Arguments -----------------------
      INTEGER N, IP, IQ
C ------------------------- Local Variables -------------------
      INTEGER M,M2,LM
      M=IP+IQ+N
      M2=M/2
      LM=( M2*(M2+1) )/2
      LM=LM*(M2+1)+ ( LM*(M2+2) )/3
C LM calculated as above to avoid overflow.
C  If MMAX=200, then LM = 681750
      LM=LM+(M-2*M2)*(M2+1)*(M2+1)
      INDX2=LM+(N-1)*(M-N+3)+IQ
      RETURN
      END