PROGRAM trnro1 C This program solves the TRaNsverse ROtation of two spheres. C It calculates coefficients for the resistance functions C C B B C Y (s,lambda) and Y (s,lambda) and C 11 12 C C C C C Y (s,lambda) and Y (s,lambda) defined in the paper C 11 12 C C Jeffrey, D.J. and Onishi, Y. 'Calculation of the resistance and C mobility functions for two unequal spheres in low-Reynolds-number C flow', J. Fluid Mech. Vol. 139 (1984) 261-290, C C H H C and Y (s,lambda) and Y (s,lambda) defined in the paper C 11 12 C C Jeffrey, D.J. 'The calculation of the low Reynolds number C resistance functions for two unequal spheres', Phys. Fluids (1992). C C The program calculates P, V and Q from equations (4.6)-(4.11) in Jeffrey C and Onishi (1984). As in axitr1, we reduce memory requirements by using C INDX1 to map non-zero array elements to a one-dimensional (dense) array. C Earlier versions which reduced storage requirements by not saving V C have been discarded, since memory is no longer a constraint. C C As in axitr1, if we want terms up to S**(-MAXS), we set MMAX=MAXS+1. C Dimension P,V,Q to INDX1(2,1,MMAX-2). Dimension the combinatorial C C (n+is) ( MMAX+1 MMAX+1 ) C factor FAC1(N,IS) = ( ) to FAC1 ( ------ , ------ ) . C (n+1 ) ( 2 2 ) C C The dimensions given below are for MAXS = 15. In trnro2.f there is a C table of dimensions for different MAXS. C DOUBLE PRECISION P(0:475),V(0:475),Q(0:475),FAC1(8,8),XN DOUBLE PRECISION FCTR1,FCTR2,FCTR3,FCTR4,FCTRQ,FCTRV,FACTOR,XS DOUBLE PRECISION PSUM,VSUM,QSUM,FAC INTEGER MMAX,MAXS,NM,N,IS,M,NQ,IQ,INDX1,JS,IP,KS,J COMMON/ARRAYS/P,V,Q,FAC1 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 FAC1 DO 10 N=1,(MMAX+1)/2 XN=DBLE(N) FAC1(N,1)=1D0 DO 10 IS=2,(MMAX+1)/2 XS=DBLE(IS) 10 FAC1(N,IS)=FAC1(N,IS-1)*(XN+XS)/(XS-1D0) C We first set the initial conditions for rotation. P(INDX1(1,0,0))=0D0 V(INDX1(1,0,0))=0D0 Q(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+.5D0)/(XN+1D0) FCTR2=XN*(XN-.5D0)/(XN+1D0) FCTR3=XN*(2D0*XN*XN-.5D0)/(XN+1D0) FCTR4=(8D0*XN*XN-2D0)/(3D0*XN+3D0) FCTRV=2D0*XN/((XN+1D0)*(2D0*XN+3D0)) FCTRQ=3D0/(2D0*XN*(XN+1D0)) 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+1 JS=(IQ+1)/2 KS=IQ/2 PSUM=0D0 VSUM=0D0 QSUM=0D0 C If JS is less than 1, we skip the loop on IS and set the elements to 0 IF(JS.LT.1)GO TO 250 DO 200 IS=1,JS XS=DBLE(IS) FAC=FAC1(N,IS) C We search the summations in Jeffrey & Onishi (1984), finding first C all INDX1es 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. The range of N is chosen to keep IP-N+1 C non-negative, 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 FACTOR=2D0*(XN*XS+1D0)**2/(XN+XS) PSUM=PSUM+FAC*FCTR1*(4D0+XN*XS-FACTOR) & *P(INDX1(IS,IQ-IS,IP-N+1))/(XS*(2D0*XS-1D0)) IF(IP.GE.N)THEN QSUM=QSUM-FAC*FCTRQ*P(INDX1(IS,IQ-IS,IP-N))/XS ENDIF IF(IP.GT.N)THEN PSUM=PSUM+FAC*FCTR2*P(INDX1(IS,IQ-IS,IP-N-1)) C A typographical error in equation (4.9) of Jeffrey & Onishi (1984) C has been corrected in the following statement. VSUM=VSUM+FAC*FCTRV*P(INDX1(IS,IQ-IS,IP-N-1)) ENDIF IF(IS.LE.KS)THEN PSUM=PSUM-FAC*FCTR4*Q(INDX1(IS,IQ-IS-1,IP-N+1)) IF(IP.GE.N)THEN QSUM=QSUM+FAC*XS*Q(INDX1(IS,IQ-IS-1,IP-N))/(XN+1D0) ENDIF ENDIF IF(IS.LT.JS)THEN PSUM=PSUM+FAC*FCTR3*V(INDX1(IS,IQ-IS-2,IP-N+1))/(2D0*XS+1D0) ENDIF 200 CONTINUE 250 P(INDX1(N,IP,IQ))=PSUM V(INDX1(N,IP,IQ))=PSUM+VSUM 300 Q(INDX1(N,IP,IQ))=QSUM 400 CONTINUE 500 CONTINUE C The full set of coefficients having now been calculated, I collect the C ones I wish to keep in the fronts of the arrays ready for writing. C I also enforce non-dimensionalizations and other consequences C of the definitions of the resistance functions. C C The definition of RYB means that this program calculates C RYB(1,1) and RYB(2,1) C whereas programs TRNTRx calculate RYB(1,1) and RYB(1,2). To show that C they are equivalent I map the RYB coefficients using a two-step process. C First (3/2)P(INDX1(1,m-q,q)) -> V( m(m+1)/2+q ). C Then for even M, the change to RYB(1,2) requires C V( m(m+1)/2+q ) <--> V( m(m+1)/2+m-q+1) for q< (m+1)/2. C I map the RYC coefficients also in 2 steps. Firstly C Q(INDX1(1,m-q,q)) -> Q( m(m+1)/2+q ). C Then for M odd, Q(m(m+1)/2+q-1) -> Q(m(m+1)/2+q) while 1.le.q.le.m C The M-loop starts at 1 because V(0) and Q(0) are set separately. C DO 600 M=1,MAXS NM=(M*(M+1))/2 DO 600 IQ=0,M Q(NM+IQ)=Q(INDX1(1,M-IQ,IQ)) 600 V(NM+IQ)=1.5D0*P(INDX1(1,M-IQ,IQ)) C The stresslet gives us YH11 and YH12, shifted, like YC12, one power C of lambda. We therefore extract the coefficients to P. The element C P(2,0,M-1) is not calculated and is set to zero separately. DO 620 M=1,MAXS NM=(M*(M+1))/2 DO 610 IQ=0,M-1 610 P(NM+IQ)=-0.375D0*P(INDX1(2,M-IQ,IQ)) 620 P(NM+M) = 0D0 C Because of the definition of the problem, both YC and YH must be shifted C one power when M is odd. DO 660 M=1,MAXS,2 NM=(M*(M+1))/2 C We must start at the top. DO 650 NQ=M,1,-1 Q(NM+NQ)=Q(NM+NQ-1) 650 P(NM+NQ)=P(NM+NQ-1) Q(NM)=0D0 660 P(NM)=0D0 C This loop reflects the coefficients of RYB when M is even. DO 670 M=2,MAXS,2 NM=(M*(M+1))/2 N=M/2 DO 670 IQ=1,N PSUM=V(NM+IQ) V(NM+IQ)=V(NM+M-IQ+1) 670 V(NM+M-IQ+1)=PSUM C C Now that the calculations are complete, I shall change the first element C of the output to a signature value. Later when I sum the series C using standard subroutines, I shall use this element to check that C the correct coefficients have been read in. The real value of the C element will be supplied by the subroutines. C C The signature of RYB is -3, of RYC -4, and of RYH -8. C P(0)=-8D0 Q(0)=-4D0 V(0)=-3D0 C C All coefficients to M=15 are printed here for comparison with the C published results. The fractional versions are not always in C simplest form. C C We follow equations (5.3), (5.4) to list the coefficients given in the C paper for YB. Note that the coefficients stored in any file rybxxx.cof C by version 3 of this program are the Q(1,k-q,q) terms before they are C multiplied by the 2**k factor. C WRITE (6,700) 700 FORMAT (1X, &'RESISTANCE FUNCTIONS DEFINED IN JEFFREY & ONISHI, J.F.M. 1984.') WRITE (6,710) 710 FORMAT(/1X, & 'NON-ZERO COEFFICIENTS FOR RYB11 AND RYB12 (EQS 5.3,5.4)', & //28X,'AS FLOATING POINT AND AS FRACTION.', & /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) THEN IF (M.LE.3) THEN WRITE (6,770) IQ,XN ELSE J=2**(M-3) XS=XN*DBLE(J) WRITE (6,770) IQ,XN,XS,J ENDIF ENDIF 770 FORMAT(10X,'LAMBDA(',I2,') : ',F15.5,10X,F13.0,'/',I4) 780 CONTINUE 790 CONTINUE C C We now follow equations (7.7), (7.8) to list the coefficients given in C paper for YC. Note that the coefficients stored in any file rycxxx.cof C by version 3 of this program are the Q(1,k-q,q) terms without the 2**k C factor. C WRITE (6,800) 800 FORMAT(//1X,'NON-ZERO COEFFICIENTS FOR RYC (EQS 7.7,7.8)', & //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=Q(NM+IQ)*DBLE(2**M) IF (XN.NE.0D0) THEN IF (M.LE.6) THEN WRITE (6,870) IQ,XN ELSE J=2**(M-6) XS=XN*DBLE(J) WRITE (6,870) IQ,XN,XS,J ENDIF ENDIF 870 FORMAT(10X,'LAMBDA(',I2,') : ',F15.5,10X,F13.0,'/',I4) 880 CONTINUE 890 CONTINUE C C We now print the coefficients of YH. Note that the coefficients stored C in any file ryhxxx.cof are the P(2,k-q,q) terms without the 2**k factor. C WRITE(6,895) 895 FORMAT(//1X, &'RESISTANCE FUNCTIONS defined in Jeffrey, Phys. Fluids, 1992.') WRITE (6,900) 900 FORMAT(/1X,'NON-ZERO COEFFICIENTS FOR RYH (EQS 34)', & //28X,'AS FLOATING POINT AND AS FRACTION.', & /1X,'FOR S**( 0)') DO 990 M=1,15 WRITE (6,950) M 950 FORMAT (1X,'FOR S**(-',I2,')') NM=((M+1)*M)/2 DO 980 IQ=0,M XN=P(NM+IQ)*DBLE(2**M) IF (XN.NE.0D0) THEN IF (M.GT.9) THEN J=2**(M-9) XS=XN*DBLE(J) WRITE (6,970) IQ,XN,XS,J ELSE WRITE (6,970) IQ,XN ENDIF ENDIF 970 FORMAT(10X,'LAMBDA(',I2,') : ',F15.5,10X,F13.0,'/',I4) 980 CONTINUE 990 CONTINUE C C Previously the equal spheres case was written out separately using C the commented out code below. C WRITE(6,1000) C1000 FORMAT(/1X,'RYB FOR EQUAL SPHERES') C DO 1040 M=0,15 C VSUM=0D0 C NM=((M+1)*M)/2 C DO 1030 IQ=0,M C VSUM=VSUM+V(NM+IQ) C1030 CONTINUE CC VSUM=VSUM*2D0**(-M) C WRITE(6,1035)M,VSUM C1035 FORMAT (1X,' FOR (2A/R)**',I2,' =',E21.14) C1040 CONTINUE C WRITE(6,1050) C1050 FORMAT(/1X,'RYC FOR EQUAL SPHERES') C DO 1090 M=0,15 C QSUM=0D0 C NM=((M+1)*M)/2 C DO 1070 IQ=0,M C QSUM=QSUM+Q(NM+IQ) C1070 CONTINUE C QSUM=QSUM*2D0**(-M) C WRITE(6,1075)M,QSUM C1075 FORMAT (1X,' FOR (2A/R)**',I2,' =',E21.14) C1090 CONTINUE C WRITE(6,1150) C1150 FORMAT(/1X,'RYH FOR EQUAL SPHERES') C DO 1190 M=0,15 C PSUM=0D0 C NM=((M+1)*M)/2 C DO 1170 IQ=0,M C PSUM=PSUM+P(NM+IQ) C1170 CONTINUE C PSUM=PSUM*2D0**(-M) C WRITE(6,1175)M,PSUM C1175 FORMAT (1X,' FOR (2A/R)**',I2,' =',E21.14) C1190 CONTINUE STOP END INTEGER FUNCTION INDX1(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-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 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