PROGRAM prpst2 C This program solves the PeRPendicular STraining deformation of two spheres. C It calculates coefficients for the resistance functions C M M C and Z (s,lambda) and Z (s,lambda) defined in C 11 12 C C Jeffrey, D.J., 'The calculation of the low Reynolds number resistance C functions for two unequal spheres', Phys. Fluids (1992). C C The program calculates P, V and Q from equations (73)-(76) C in Jeffrey (1992). As in axist1, we reduce memory requirements by using C INDX2 to map non-zero array elements to a one-dimensional (dense) array. C C As in axist1, if we want terms up to S**(-MAXS), we set MMAX=MAXS+2. C The memory requirements are C C MAXS MMAX DIMENSION P, V, Q MEMORY (for one array) FAC1 C ---- ---- ----------------- ------ (per 8 bytes) ---- C 15 17 558 5KB 9,9 C 20 22 1121 9KB 11,11 C 50 52 12856 100KB 26,26 C 100 102 92581 740KB 51,51 C 200 202 702656 5.6MB 101,101 C 300 302 2330231 18.6MB 151,151 C DOUBLE PRECISION P(0:2330231),V(0:2330231),Q(0:2330231) DOUBLE PRECISION C2(151,151),XN,XS DOUBLE PRECISION F1,F2,F3,F4,FV,FP DOUBLE PRECISION PSUM,VSUM,QSUM,FAC,XNS INTEGER MMAX,MAXS,NM,N,IS,M,NQ,IQ,INDX2,JS,IP,KS,ISUB LOGICAL EVEN MAXS=300 MMAX=MAXS+2 C First tabulate C2 DO 10 N=1,(MMAX+1)/2 XN=DBLE(N) C2(N,2)=1D0/(1D0+XN) DO 10 IS=3,(MMAX+1)/2 XS=DBLE(IS) 10 C2(N,IS)=C2(N,IS-1)*(XN+XS)/(XS-2D0) C We first set the initial conditions DO 5 M=0,5 P(M)=0D0 V(M)=0D0 5 Q(M)=0D0 ISUB=INDX2(2,0,0) P(ISUB)=1D0 V(ISUB)=1D0 Q(ISUB)=0D0 C We start at N=2, IQ=0, IP=M and proceed keeping M=IP+IQ+N 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=2 is calculated. DO 500 M=3,MMAX NM=M/2+1 IF(M.EQ.MMAX)NM=2 C Note that N and IS start at 2, because N=1 does not exist for this case DO 400 N=2,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) F1=(XN+.5D0) F2=XN*(XN-.5D0) F3=XN*(2D0*XN*XN-0.5D0) F4=8D0*XN*XN-2D0 FV=XN/(XN+1.5D0) C It was in the next term that a factor 2 was missed in early versions C FQ=2d0/XN 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-1=IQ-JS. C We also need KS defined by KS-1=IQ-KS-1. IP=M-IQ-N JS=IQ/2+1 KS=(IQ+1)/2 PSUM=0D0 VSUM=0D0 QSUM=0D0 IF ( IQ.EQ.0 .OR. IQ.EQ.NQ) GO TO 250 DO 200 IS=2,JS XS=DBLE(IS) FAC=C2(N,IS) C We search the summations in Jeffrey (1992), finding first C all INDX2es that start IS,IQ-IS after which we test the final argument C which can be IP-N+1, IP-N or IP-N-1. c The range of N is chosen to keep IP-N+1 non-negative, C and the range of IS is chosen to keep the combination IS,IQ-IS C legal. Thus the other possibilities require IF tests and jumps. C XNS=XN*XS FP=(XNS+16D0)-2D0*(XNS+4D0)*(XNS+1)/(XN+XS) FP=FP/(XS*(2D0*XS-1D0)) PSUM=PSUM+FAC*F1*FP*P( INDX2(IS,IQ-IS,IP-N+1) ) IF(IP.GE.N) & QSUM=QSUM-FAC*2D0*P( INDX2(IS,IQ-IS,IP-N) )/XNS IF(IP.GT.N) THEN ISUB=INDX2(IS,IQ-IS,IP-N-1) PSUM=PSUM+FAC*F2*P(ISUB) VSUM=VSUM+FAC*FV*P(ISUB) ENDIF IF(IS.LE.KS) THEN PSUM=PSUM-FAC*F4*Q( INDX2(IS,IQ-IS-1,IP-N+1) ) IF(IP.GE.N) & QSUM=QSUM+FAC*XS*Q( INDX2(IS,IQ-IS-1,IP-N) ) ENDIF IF(IS.LT.JS) & PSUM=PSUM+FAC*F3*V( INDX2(IS,IQ-IS-2,IP-N+1) )/(2D0*XS+1D0) 200 CONTINUE 250 ISUB=INDX2(N,IP,IQ) P(ISUB)=PSUM V(ISUB)=P(ISUB)+VSUM 300 Q(ISUB)=QSUM 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 outputting. C I also enforce non-dimensionalizations and other consequences C of the definitions of the resistance functions. C From Eq 77 we see that the odd powers must be shifted one. DO 610 M=0,MAXS EVEN=M.EQ.2*(M/2) NM=(M*(M+1))/2 DO 610 IQ=0,M IF (EVEN) THEN P(NM+IQ)=P( INDX2(2,M-IQ,IQ) ) ELSE P(NM+IQ)=P( INDX2(2,M-IQ+1,IQ-1) ) ENDIF 610 CONTINUE C C Now that the calculations are complete, I shall change the first element C of the output. Later when I sum the series using standard subroutines, C I shall use this element to check that the correct coefficients have C been read in. The real value of the element will be supplied by the C subroutines. C C The signature for RZM is -11. C P(0)=-11D0 MMAX=((MAXS+1)*(MAXS+2))/2 -1 700 FORMAT(D22.16) OPEN (1,FILE='rzmcof.dat',STATUS='NEW') WRITE(1,705) MAXS 705 FORMAT(' RZM',I5) DO 710 M=0,MMAX 710 WRITE (1,700) P(M) C C THE CYBERS AND ETA DO NOT NEED A SPACE IN FORMAT STATEMENTS FOR C DISK FILES. C CLOSE (1) C C Previously the equal spheres case was written out separately using C the commented out code below. C C900 FORMAT(D27.21) C OPEN (11,FILE='rzmeqcof.dat') C DO 1090 M=0,MAXS C PSUM=0D0 C NM=((M+1)*M)/2 C DO 1070 IQ=0,M C PSUM=PSUM+P(NM+IQ) C1070 CONTINUE C PSUM=PSUM*2D0**(-M) C WRITE(11,900)PSUM C1090 CONTINUE C CLOSE (11) STOP END C 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