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 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