PROGRAM AXITR1 C This program solves the AXIsymmetric TRanslation of two spheres. C It calculates coefficients for the resistance functions C C A A C X (s,lambda) and X (s,lambda) C 11 12 C C defined in the paper Jeffrey, D.J. and Onishi, Y. 'Calculation of the C resistance and mobility functions for two unequal spheres in C low-Reynolds-number flow', J. Fluid Mech. Vol. 139 (1984) 261-290. C C It also calculates coefficients for the resistance functions C C G G C X (s,lambda) and X (s,lambda) C 11 12 C C defined in the paper Jeffrey,'The calculation of the low Reynolds number C resistance functions for two unequal spheres', Phys. Fluids (1992). C The program calculates P and V from equations (3.6)-(3.9) in Jeffrey and C Onishi (1984). To reduce memory requirements, a function INDX1 is used to C map non-zero array elements to a one-dimensional (dense) array. C The coding and output of this program are intended to display the C connection with published results as clearly as possible. C C The main counting variable is M=N+P+Q-1, where N, P, Q are array indices. C If the maximum power we wish to find is S**(- MAXS), then the maximum value C of M is MMAX=MAXS+1. Dimension P,V to INDX1(2,1,MMAX-2). Also dimension C C (n+is) C the combinatorial factor FAC0(n,is)= ( ) to ( (MMAX+1)/2, (MMAX+1)/2 ). C ( n ) C C The dimensions given below are for MAXS = 15. In AXITR2.FOR there is a C table of dimensions for different MAXS. C DOUBLE PRECISION P(0:475),V(0:475),FAC0(8,8),XN,XS DOUBLE PRECISION 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 C This COMMON block reduces disk space requirements with many compilers C and helps memory management on some computers. MAXS=15 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) FAC0(N,IS)=FAC0(N,IS-1)*(XN+XS)/XS 10 CONTINUE C We set the initial conditions. P(INDX1(1,0,0))=1D0 V(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 and 2 are calculated. DO 500 M=1,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-1=N+N-1+NQ-1 NQ=M-2*N+2 XN=DBLE(N) FCTR1=XN*(XN+.5D0)/(XN+1D0) FCTR2=XN*(XN-.5D0)/(XN+1D0) FCTR3=XN*(2D0*XN*XN-.5D0)/(XN+1D0) FCTRV=2D0*XN/((XN+1D0)*(2D0*XN+3D0)) 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. DO 300 IQ=0,NQ IP=M-IQ-N+1 PSUM=0D0 VSUM=0D0 C If Q=0, the sum does not need to be done. This would happen anyway C in FORTRAN 77 because JS=0 when Q=0. However, in FORTRAN IV a test C was needed. IF (IQ.EQ.0D0) GO TO 250 JS=(IQ+1)/2 DO 200 IS=1,JS XS=DBLE(IS) FAC=FAC0(N,IS) C The way the loops and the index mapping are set up, the first term C in (3.9) will always be legal, but the other terms in (3.8)-(3.9) must be C tested to make certain that the INDX1 is legal (an illegal argument to C INDX1 means that the term can be shown analytically to be zero). C PSUM=PSUM+FAC*FCTR1*(2D0*XN*XS-XN-XS+2D0)* & P(INDX1(IS,IQ-IS,IP-N+1))/((XN+XS)*(2D0*XS-1D0)) 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 C The largest value of IS for the next term is given by IS-1=IQ-IS-2 C Thus IS=(IQ-1)/2=(IQ+1)/2-1=JS-1. Thus we skip this if IS=JS. IF(IS.LT.JS) PSUM=PSUM-FAC*FCTR3* & V(INDX1(IS,IQ-IS-2,IP-N+1))/(2D0*XS+1D0) 200 CONTINUE 250 P(INDX1(N,IP,IQ))=PSUM V(INDX1(N,IP,IQ))=PSUM+VSUM 300 CONTINUE 400 CONTINUE 500 CONTINUE C C The force is the sum over m and q of P(1,m-q,q)*t1**(m-q)*t2**q . C I extract these coefficients and store them in the front of array V, C which is no longer needed. Notice the array V is used before P. C V( m(m+1)/2+q ) = P(INDX1(1,m-q,q)) 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 C C The stresslet is the sum over m and q of (3/4)P(2,m-q,q)*t1**(m-q)*t2**q. C I extract the coefficients and store them in the front of array P. C P( m(m+1)/2+q ) = (3/4) P(INDX1(2,m-q,q)) for 0.lt.q.lt.m C P( m(m+1)/2+m ) = 0. C Remember P(2,0,m-1) was not stored, so in particular P(2,0,0) was not. C P(0)=0D0 DO 620 M=1,MAXS NM = (M*(M+1))/2 DO 610 IQ=0, M-1 610 P(NM+IQ)=0.75D0*P( INDX1(2,M-IQ,IQ) ) P(NM+M) =0D0 620 CONTINUE C C We now follow equation (3.15) 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 WRITE (6,700) 700 FORMAT(1X, &'RESISTANCE FUNCTIONS RXA FROM JEFFREY AND ONISHI, J.F.M. 1984') WRITE (6,705) 705 FORMAT(1X,'NON-ZERO COEFFICIENTS FOR USE IN EQS (3.13-3.14)'/) DO 790 M=0,15 WRITE (6,750) M 750 FORMAT (1X,'FOR S**(-',I2,')') NM=((M+1)*M)/2 DO 790 IQ=0,M XN=V(NM+IQ)*DBLE(2**M) IF (XN.NE.0D0) WRITE (6,770) IQ,XN 770 FORMAT(10X,'LAMBDA(',I2,') = ',F15.2) 790 CONTINUE C C We now list the coefficients to give in the paper. C Note that the coefficients stored in any file RXGxxx.DAT are the C P(2,M-q,q) terms without the 2**M factor. C WRITE (6,800) 800 FORMAT(/1X, & 'RESISTANCE FUNCTIONS RXG FROM JEFFREY, PHYS. FLUIDS 1992') WRITE (6,805) 805 FORMAT(1X,'NON-ZERO COEFFICIENTS FOR USE IN EQS (18a,18b)'/) DO 860 M=0,15 WRITE (6,830) M 830 FORMAT (1X,'FOR S**(-',I2,')') NM=((M+1)*M)/2 DO 860 IQ=0,M XN=P(NM+IQ)*DBLE(2**M) IF (XN.NE.0D0) WRITE (6,840) IQ,XN 840 FORMAT(10X,'LAMBDA(',I2,') : ',F15.2) 860 CONTINUE C Previously the equal spheres case was written out separately using the C commented out code below. C C WRITE(6,900) C900 FORMAT(/1X,'RXA FOR EQUAL SPHERES: COEFFICIENTS FOR (2A/R)**N') C DO 940 M=0,15 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(6,935)M,VSUM C935 FORMAT (1X,' FOR (2A/R)**',I2,' =',E21.14) C940 CONTINUE C WRITE(6,950) C950 FORMAT(/1X,'RXG FOR EQUAL SPHERES: COEFFICIENTS FOR (2A/R)**N') C DO 990 M=0,15 C PSUM=0D0 C NM=((M+1)*M)/2 C DO 970 IQ=0,M C PSUM=PSUM+P(NM+IQ) C970 CONTINUE C PSUM=PSUM*2D0**(-M) C WRITE(6,975)M,PSUM C975 FORMAT (1X,' FOR (2A/R)**',I2,' =',E21.14) C990 CONTINUE STOP END C INTEGER FUNCTION INDX1(N,IP,IQ) C C------------------------------ARGUMENTS------------------------------ INTEGER N,IP,IQ C--------------------------------------------------------------------- C 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 = 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 M,M2,LM C--------------------------------------------------------------------- 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