PROGRAM axiro1
C This program calculates coefficients for the resistance functions
C
C   C                   C
C  X  (s,lambda)  and  X  (s,lambda) defined in the paper
C   11                  12
C
C  Jeffrey, D.J. and Onishi, Y. 'Calculation of the resistance and
C     mobility functions for two unequal spheres in low-Reynolds-number
C     flow', J. Fluid Mech. Vol. 139 (1984) 261-290.
C
C  The program calculates Q from equations (6.4)-(6.5) in that paper.
C  To reduce memory requirements, a function INDX1 is used to map
C  non-zero array elements to a one-dimensional (dense) array.
C  This mapping has been programmed in a computationally slow form
C  so that the connection between the paper and the program is easy to
C  follow and check. Similarly the output is designed to display the
C  connection with the published results as clearly as possible.
C
C The main counting variable is M=N+P+Q-1, where N, P, and Q are array
C indices. The maximum value of M is MMAX. Dimension Q to INDX1(1,0,MMAX).
C
C                                            (n+s)
C Dimension the combinatorial factor  FAC0 = (   ) to (MMAX/2+1,MMAX/2+1).
C                                            ( n )
C
C The dimensions given below are for MMAX = 15. In AXIRO2.FOR there is a
C table of dimensions for different MMAX.
C
      DOUBLE PRECISION Q(0:387),FAC0(8,8),QSUM,XN,XS
      INTEGER MMAX,NM,N,IS,M,NQ,IQ,INDX1,KS,IP
      COMMON/ARRAYS/Q,FAC0
C This COMMON block reduces disk space requirements with many compilers
C  and helps memory management on some computers.
      MMAX=15
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
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 is calculated.
      DO 500 M=1,MMAX
        NM=M/2+1
        IF(M.EQ.MMAX)NM=1
        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=DBLE(N)
          DO 300 IQ=0,NQ
C Now that M, N, Q are specified, I know P (denoted IP) as well.
C  We obtain KS from the equation KS-1=IQ-KS-1.
            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
C
C We now follow equation (6.7) to list the coefficients as given in the paper
C The coefficients written in any file rxaxxx.dat by version 2 of
C this program are the P(1,M-q,q) terms without the 2**M factor.
C
C
      WRITE (6,700)
700   FORMAT(1X,
     & 'RESISTANCE FUNCTIONS IN JEFFREY AND ONISHI, J.F.M. 1984')
      WRITE (6,710)
710   FORMAT(/1X,'COEFFICIENTS FOR RXC (EQ 6.7)'//1X,'     FOR S**( 0)')
      WRITE (6,720)
720   FORMAT(1X,'LAMBDA( 0) :         1.')
      DO 800 M=1,15
        WRITE (6,750) M
750     FORMAT (1X,'     FOR S**(-',I2,')')
        NM=((M+1)*M)/2
        DO 800 IQ=0,M
          XN=Q(NM+IQ)*DBLE(2**M)
          IF (XN.NE.0D0) THEN
            WRITE (6,770) IQ,XN
770         FORMAT(1X,'LAMBDA(',I2,') : ',F15.2)
          ENDIF
800   CONTINUE
      STOP
      END
      INTEGER FUNCTION INDX1(N,IP,IQ)
C ------------------------------ Arguments -----------------------------
      INTEGER N,IP,IQ
C This function maps 3-dimensional arrays onto 1-dimensional arrays.
C  For each array element P(n,p,q) we define M=n+p+q-1, which is the primary
C  variable. For each value of M, we find that the array elements are
C  non-zero only for 1=< n =<NM = M/2+1. Further, for each fixed M and n,
C  the array elements that are non-zero are:
C                     P(n,m-n+1,0) -----> P(n,n-1,m-2n+2).
C  INDX1 uses these facts to set up a triangular scheme for each M to convert
C   it to a linear array. The array is filled as follows:
C   n,p,q=   1,M,0    1,M-1,1  .......................  1,1,M-1  1,0,M
C                2,M-1,0  2,M-2,1 ...............  2,1,M-2
C                **     ....................
C                **     n,M-n+1,0  ...  n,n-1,M-2n+2
C Note that the integer arithmetic relies on correct divisibility.
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