PROGRAM axist1
C This program solves the AXIsymmetric STraining deformation of
C two spheres. It calculates coefficients for the resistance functions
C
C G G
C X (s,lambda) and X (s,lambda)
C 11 12
C
C and the resistance functions
C
C M M
C X (s,lambda) and X (s,lambda)
C 11 12
C
C defined in Jeffrey, 'The calculation of the low Reynolds number
C resistance functions for two unequal spheres', Phys. Fluids (1992).
C
C The program calculates P and V from equations (44)-(46)
C in Jeffrey (1992). To reduce memory requirements, a function INDX2
C is used to map non-zero array elements to a one-dimensional 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, 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+2. Dimension P,V to INDX2(2,0,MAXS). Also dimension
C
C (n+is) /
C the combinatorial factor C0(n,is)= ( ) / (n+1)
C ( n )/
C to ( (MMAX+1)/2, (MMAX+1)/2 ).
C The dimensions given below are for MAXS = 15. In axist2.f there is a
C table of dimensions for different MAXS.
C
DOUBLE PRECISION P(0:558),V(0:558),C0(9,9),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=15
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**2-.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 The term corresponding to IQ = NQ is INDX2(N,N-2,NQ) and
C this is always 0, except when N=2. Moreover, if we plug this combination
C of subscripts into our recurrence formulae, we find that we strike a
C subscript INDX2(IS, NQ-IS, -1), which we have not included in our
C numbering scheme, so we must implement special tests.
C If IQ = NQ we skip the summation. And, since we are on the subject,
C the case IQ=0 is also special.
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 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 INDX2 is legal (an illegal argument to
C INDX2 means that the term can be shown analytically to be zero).
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 The largest value of IS for the next term is given by IS-2 = IQ-IS-2.
C Thus IS = IQ/2 = JS-1. Thus we skip this if IS = JS.
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 Having calculated all the coefficients , I now bring the ones I wish to
C keep to the front of the array, enforcing non-dimensionalizations as I do.
C The force is the sum over m and q of 5*P(1,m-q,q)*t1**(m-q)*t2**q.
C Because we defined the problem using (41) in Jeffrey (1992), I get
C a factor of lambda to watch for.
C Also the G tensor here is the transpose of the one in AXITR1, so we
C must use the symmetry relations to get the output to agree with the
C RXG coefficients in axitr1.
C Thus for RXG I use the mappings:
C for m even 5*P(INDX2(1,Q-1,M-Q+1)) -> V(usual),
C for m odd 5*P(INDX2(1,M-Q,Q)) -> V( usual)
DO 610 M=0,MAXS
EVEN= M .EQ. 2*(M/2)
NM=(M*(M+1))/2
DO 610 IQ=0,M
IF (EVEN) THEN
V(NM+IQ)=5D0*P(INDX2(1,IQ-1,M-IQ+1))
ELSE
V(NM+IQ)=5D0*P(INDX2(1,M-IQ,IQ))
ENDIF
610 CONTINUE
C The mapping for RXM is given in Jeffrey (1992) equations (47a,b).
C For even m, I use the mapping P( INDX2(2,m-q,q) ) -> P(m(m+1)/2+q+1)
C for odd m, I use P( INDX2(2,m-q+1,q-1) ) -> P( m(m+1)/2+q+1 ).
DO 620 M=0,MAXS
EVEN= M .EQ. 2*(M/2)
NM=(M*(M+1))/2
DO 620 IQ=0,M
IF (EVEN) THEN
P(NM+IQ)=P(INDX2(2,M-IQ,IQ))
ELSE
IF (IQ.EQ.0) THEN
P(NM+IQ) = 0D0
ELSE
P(NM+IQ)=P(INDX2(2,M-IQ+1,IQ-1))
ENDIF
ENDIF
620 CONTINUE
C Now that the calculations are complete, I shall change the first elements
C of the output. Later when I sum the series using a standard subroutine,
C I shall use these elements to check that the correct coefficients have
C been read in. The real value of the element will be supplied by the
C subroutine.
C We now list the coefficients as given in the paper. Note that the
C coefficients stored in any file rxgxxx.dat by version 2 of
C this program are the P(1,M-q,q) terms without the 2**M factor.
V(0)=-6D0
P(0)=-9D0
WRITE (6,700)
700 FORMAT (1X,
& 'RESISTANCE FUNCTIONS RXG AND RXM DEFINED IN JEFFREY,',
& ' PHYS. FLUIDS (1992)'//1X,
& 'NON-ZERO COEFFICIENTS FOR RXG11 AND RXG12 (EQUNS 18)',
& //1X,'FOR S**( 0)')
DO 790 M=1,15
WRITE (6,750) M
750 FORMAT (1X,'FOR S**(-',I2,')')
NM=((M+1)*M)/2
DO 780 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)
780 CONTINUE
790 CONTINUE
C We now list the coefficients as given in the paper. Note that the
C coefficients stored in any file rxmxxx.dat by version 2 of
C this program are the P(2,M-q,q) terms without the 2**M factor.
WRITE (6,800)
800 FORMAT(/1X,
& 'NON-ZERO COEFFICIENTS FOR RXM11 AND RXM12 (EQUNS 47)',
& //28X,'AS FLOATING POINT AND AS FRACTION.'
& /1X,'FOR S**( 0)')
WRITE (6,805)
805 FORMAT(10X,'LAMBDA( 0) : 1.')
DO 890 M=1,15
WRITE (6,850) M
850 FORMAT (1X,'FOR S**(-',I2,')')
NM=((M+1)*M)/2
DO 880 IQ=0,M
XN=P(NM+IQ)*DBLE(2**M)
IF(XN.NE.0D0) THEN
IF (M.LT.8) THEN
WRITE (6,870) IQ,XN
ELSE
XS=XN*5D0
WRITE (6,870) IQ,XN,XS
ENDIF
ENDIF
870 FORMAT(10X,'LAMBDA(',I2,') : ',F15.2,10X,F13.0,'/5')
880 CONTINUE
890 CONTINUE
C
C Previously the equal spheres case was written out separatly using the
C commented out code below.
C WRITE(6,900)
C900 FORMAT(/1X,'RXG 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,'RXM 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
INTEGER FUNCTION INDX2(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, 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 INDX2 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-1,0 1,M-2,1 ....................... 1,0,M-1 1,-1,M
C ***** 2,M-2,0 2,M-3,1 ............... 2,0,M-2
C ***** ******* ....................
C ***** ******* n,M-n,0 ... n,n-2,M-2n+2
C Note that the integer arithmetic relies on correct divisibility.
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 is 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