PROGRAM trnst1 C This program solves the TRaNsverse STraining deformation of two spheres. C It calculates coefficients for the resistance functions C C G G C Y (s,lambda) and Y (s,lambda) C 11 12 C C H H C Y (s,lambda) and Y (s,lambda) C 11 12 C C M M C and Y (s,lambda) and Y (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 (58)-(61) in C 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 C1(N,IS) = ( ) / (n+1) to C1 ( ------ , ------ ) . C (n+1 )/ ( 2 2 ) C C The dimensions given below are for MAXS = 15. In trnst2.f there is a C table of dimensions for different MAXS. C DOUBLE PRECISION P(0:558),V(0:558),Q(0:558),C1(9,9),XN,XS DOUBLE PRECISION FCTR1,FCTR2,FCTR3,FCTR4,FCTRV,FACTOR DOUBLE PRECISION PSUM,VSUM,QSUM,FAC INTEGER MMAX,MAXS,NM,N,IS,M,NQ,IQ,INDX2,JS,IP,KS,ISUB,J LOGICAL EVEN COMMON/ARRAYS/P,V,Q,C1 C MAXS=15 MMAX=MAXS+2 C First tabulate C1 DO 10 N=1,(MMAX+1)/2 XN=DBLE(N) C1(N,1)=1D0/(XN+1D0) DO 10 IS=2,(MMAX+1)/2 XS=DBLE(IS) 10 C1(N,IS)=C1(N,IS-1)*(XN+XS)/(XS-1D0) C We first set the initial conditions DO 5 M=0,5 P(M)=0D0 V(M)=0D0 5 Q(M)=0D0 P( INDX2(2,0,0) )=1D0 V( INDX2(2,0,0) )=1D0 Q( INDX2(2,0,0) )=0D0 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,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-1=N+N-1+NQ-1 NQ=M-2*N+2 XN=DBLE(N) FCTR1=(XN+.5D0) FCTR2=XN*(XN-.5D0) FCTR3=XN*(2D0*XN*XN-.5D0) FCTR4=(4D0*XN*XN-1D0) 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-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=1,JS XS=DBLE(IS) FAC=C1(N,IS) C We search the summations in Jeffrey (1984), 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 FACTOR=2D0*(XN*XS+1D0)**2/(XN+XS) PSUM=PSUM+FAC*FCTR1*(4D0+XN*XS-FACTOR) & *P(INDX2(IS,IQ-IS,IP-N+1))/(XS*(2D0*XS-1D0)) IF(IP.GE.N) THEN QSUM=QSUM-FAC*P(INDX2(IS,IQ-IS,IP-N))/(XN*XS) ENDIF IF(IP.GT.N) THEN ISUB=INDX2(IS,IQ-IS,IP-N-1) PSUM=PSUM+FAC*FCTR2*P(ISUB) VSUM=VSUM+FAC*FCTRV*P(ISUB) ENDIF IF(IS.LE.KS) THEN PSUM=PSUM-FAC*FCTR4*Q( INDX2(IS,IQ-IS-1,IP-N+1)) IF(IP.GE.N) THEN QSUM=QSUM+FAC*XS*Q( INDX2(IS,IQ-IS-1,IP-N) ) ENDIF ENDIF 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 ISUB=INDX2(N,IP,IQ) P(ISUB)=PSUM V(ISUB)=PSUM+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 C The mapping for RYM is given in equation (62) in Jeffrey(1992). C Because of the lambda dividing RYM12, we map the coefficients using: C m even P(INDX2(2,m-q,q)) -> P(m(m+1)/2+q+1). C m odd P(INDX2(2,m-q+1,q-1)) -> P(m(m+1)/2+q+1). C The RYG coefficients are mapped : C m even (5/3)*P(INDX2(1,q-1,m-q+1)->V(usual) C m odd (5/3)*P(INDX2(1,m-q,q)->V(usual) C The RYH coefficients are for m even -5*Q(1,m-q,q)/3 c for m odd 5*Q(1,q-2,m-q+2)/3 DO 610 M=0,MAXS EVEN=M.EQ.2*(M/2) NM=(M*(M+1))/2 DO 610 IQ=0,M IF (EVEN) THEN Q(NM+IQ)=-5D0*Q( INDX2(1,M-IQ,IQ) )/3D0 IF(IQ.EQ.0) THEN V(NM+IQ)=0D0 ELSE V(NM+IQ)=5D0*P(INDX2(1,IQ-1,M-IQ+1))/3D0 ENDIF P(NM+IQ)=P( INDX2(2,M-IQ,IQ) ) ELSE V(NM+IQ)=5D0*P( INDX2(1,M-IQ,IQ) )/3D0 IF(IQ.GE.2) THEN Q(NM+IQ)=-5D0*Q( INDX2(1,IQ-2,M-IQ+2) )/3D0 ELSE Q(NM+IQ)=0D0 ENDIF 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 V(0)=-7D0 Q(0)=-8D0 P(0)=-10D0 WRITE (6,700) 700 FORMAT (1X, & 'RESISTANCE FUNCTIONS RYG, RYH AND RYM DEFINED IN JEFFREY,', & ' PHYS. FLUIDS (1992)'//1X, & 'NON-ZERO COEFFICIENTS FOR RYG (EQS 28)', & //28X,'AS FLOATING POINT AND AS FRACTION.', & /1X,'FOR S**( 0)') C zero anyway DO 750 M=1,15 WRITE (6,730) M 730 FORMAT (1X,'FOR S**(-',I2,')') NM=((M+1)*M)/2 DO 750 IQ=0,M XN=V(NM+IQ)*DBLE(2**M) IF(XN.NE.0D0) THEN IF (M.GT.6) THEN J=2**(M-6) XS=XN*DBLE(J) WRITE (6,735) IQ,XN,XS,J ELSE WRITE (6,735) IQ,XN ENDIF ENDIF 735 FORMAT(10X,'LAMBDA(',I2,') : ',F15.5,10X,F15.2,'/',I4) 750 CONTINUE C WRITE (6,760) 760 FORMAT(/1X,'NON-ZERO COEFFICIENTS FOR RYH (EQS 36)' & //1X,'FOR S**( 0)') DO 790 M=1,15 WRITE (6,770) M 770 FORMAT (1X,'FOR S**(-',I2,')') NM=((M+1)*M)/2 DO 790 IQ=0,M XN=Q(NM+IQ)*2**M IF (XN.NE.0D0) THEN IF (M.GT.9) THEN J=2**(M-9) XS=XN*DBLE(J) WRITE (6,775) IQ,XN,XS,J ELSE WRITE (6,775) IQ,XN ENDIF ENDIF 775 FORMAT(10X,'LAMBDA(',I2,') : ',F15.5,10X,F15.2,'/',I4) 790 CONTINUE C WRITE (6,800) 800 FORMAT(/1X, & 'NON-ZERO COEFFICIENTS FOR RYM (EQS 63)', & //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 890 IQ=0,M XN=P(NM+IQ)*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,') : ',F17.4,10X,F15.2,'/5') 890 CONTINUE C C Previously the equal spheres case was written out separately using C the commented out code below. C C WRITE(6,900) C900 FORMAT(/1X,'RYG FOR EQUAL SPHERES') 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,1050) C1050 FORMAT(/1X,'RYH 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,'RYM 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 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