PROGRAM axiro2
C  Comment lines found in axiro1.f are not repeated here, and should be
C   read first. This program modifies axiro1 only in the output of the
C   coefficients for the resistance function RXC.
C   The results are written out in a form suitable for the subroutines RXC,
C   and RCX.
C   To calculate large numbers the memory requirements are:
C
C   MMAX    DIMENSION Q    MEMORY  (PER 8 BYTES)      FAC0
C   ----    -----------    ------  (for each array)  -------
C     15           387        3KB                      8,8
C     20           845        7KB                     11,11
C     50         11425       92KB                     26,26
C    100         87225      698KB                     51,51
C    200        681950      5.5MB                    101,101
C    300       2284175     18.3MB                    151,151
C
C
      DOUBLE PRECISION Q(0:2284175),FAC0(151,151),QSUM,XN,XS
      INTEGER MAXS, MMAX,NM,IS,M,NQ,INDX1,KS,N,IP,IQ
      COMMON/ARRAYS/Q,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=DBLE(N)
        FAC0(N,1)=1D0+XN
        DO 10 IS=2,(MMAX+1)/2
          XS=DBLE(IS)
   10 FAC0(N,IS)=FAC0(N,IS-1)*(XN+XS)/XS
C We set the initial conditions for rotation.
      Q(INDX1(1,0,0))=1D0
      DO 500 M=1,MMAX
        NM=M/2+1
        IF(M.EQ.MMAX)NM=1
        DO 400 N=1,NM
          NQ=M-2*N+2
          XN=DBLE(N)
          DO 300 IQ=0,NQ
            IP=M-IQ-N+1
            KS=IQ/2
            QSUM=0D0
            IF(KS.LT.1.OR.IP.LT.N) GOTO 300
            DO 200 IS=1,KS
              XS=DBLE(IS)
  200       QSUM=QSUM+FAC0(N,IS)*XS*Q(INDX1(IS,IQ-IS-1,IP-N))
  300     Q(INDX1(N,IP,IQ))=QSUM/(XN+1D0)
  400   CONTINUE
  500 CONTINUE
C
C The full set of coefficients having now been calculated, I collect the
C  ones I wish to keep to the front of the array ready for writing.
C  I also enforce non-dimensionalizations and other consequences
C  of the definition of RXC. I map the coefficients in 2 steps. Firstly
C  Q(INDX1(1,m-q,q)) -> Q( m(m+1)/2+q ).
C  Then for M odd, Q( m(m+1)/2+q ) -> Q( m(m+1)/2+q+1 ) for 0<q<=m-1
C
      DO 600 M=0,MMAX
        NM=(M*(M+1))/2
        DO 600 IQ=0,M
  600 Q(NM+IQ)=Q(INDX1(1,M-IQ,IQ))
      DO 630 M=1,MMAX,2
        NM=(M*(M+1))/2
        DO 620 NQ=M,1,-1
  620   Q(NM+NQ)=Q(NM+NQ-1)
  630 Q(NM)=0D0

      OPEN(1,FILE='rxccof.dat',STATUS='NEW')
      WRITE(1,650) MAXS
 650  FORMAT(' RXC',I5)
      MMAX=((MAXS+1)*(MAXS+2))/2 -1
      DO 690 M=0,MMAX
      WRITE (1,700)Q(M)
700   FORMAT(D22.16)
690   CONTINUE
      CLOSE(1)
      STOP
      END
C
      INTEGER FUNCTION INDX1(N,IP,IQ)
C ----------------------------- Arguments -----------------------
      INTEGER N,IP,IQ
C -------------------------- Local variables --------------------
      INTEGER M2,LM,M
      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