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