program prpst1 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 Dimension P,V,Q to INDX2(2,0,MAXS). Dimension the combinatorial C C (n+is) / ( MMAX+1 MMAX+1 ) C factor C2(N,IS) = ( ) / (n+1) to C2 ( ------ , ------ ) . C (n+2 )/ ( 2 2 ) C C The dimensions given below are for MAXS = 15. In PRPST2.FOR there is a C table of dimensions for different MAXS. C DOUBLE PRECISION P(0:558),V(0:558),Q(0:558),C2(9,9),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=15 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 WRITE (6,700) 700 FORMAT (1X, & 'RESISTANCE FUNCTIONS RZM DEFINED IN JEFFREY,', & ' PHYS. FLUIDS (1992)'//1X, & 'NON-ZERO COEFFICIENTS FOR RZM (EQS 78)', & //1X,'FOR S**( 0)') WRITE (6,705) 705 FORMAT(10X,'LAMBDA( 0) : 1.00') DO 790 M=1,15 WRITE (6,750) M 750 FORMAT (1X,'FOR S**(-',I2,')') NM=((M+1)*M)/2 DO 790 IQ=0,M XN=P(NM+IQ)*2**M IF(XN.NE.0D0) WRITE (6,770) IQ,XN 770 FORMAT(10X,'LAMBDA(',I2,') : ',F12.2) 790 CONTINUE C C Previously the equal spheres case was written out separately using C the commented out code below. C C WRITE(6,1050) C1050 FORMAT(/1X,'RZM FOR EQUAL SPHERES') C DO 1090 M=0,15 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(6,1075)M,PSUM C1075 FORMAT (1X,' FOR (2A/R)**',I2,' =',E21.14) C1090 CONTINUE 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