PROGRAM axist1
C This program solves the AXIsymmetric STraining deformation of
C two spheres. It calculates coefficients for the resistance functions
C
C Q Q
C X (s,lambda) and X (s,lambda)
C 11 12
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.FOR 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
do 690 m=1,maxs
EVEN= M .EQ. 2*(M/2)
NM = (M*(M+1))/2
v(nm)=0d0
if (even) then
do 640 iq=1,m-1
ip=m-iq-1
psum=0d0
do 630 n=1,(iq+2)/2
psum=psum+p( indx2(n,iq-n,ip) )
630 continue
640 v( nm+iq) = 2.5d0*psum
v(nm+m)=0d0
else
do 660 iq=1,m
ip=m-iq
psum=0d0
do 650 n=1,(iq+2)/2
psum=psum+p( indx2(n,iq-n-1,ip) )
650 continue
660 v( nm+iq) = 2.5d0*psum
endif
690 continue
V(0)=-13D0
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,750)
750 FORMAT(1X,
&'RESISTANCE FUNCTIONS from Jeffrey, Morris & Brady, 1993')
WRITE(6,760)
760 FORMAT(1X,'Phys. Fluids A5, 2317')
WRITE (6,800)
800 FORMAT(/1X,
& 'NON-ZERO COEFFICIENTS FOR RXQ11 AND RXQ12 (EQUNS 47)',
& //28X,'AS FLOATING POINT AND AS FRACTION.'
& /1X,'FOR S**( 0)')
WRITE (6,805)
805 FORMAT(10X,'LAMBDA( 0) : 0.')
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=V(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 separately using
C the commented out code below.
C
C WRITE(6,950)
C950 FORMAT(/1X,'RXQ 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+V(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 =<NM = M/2+1. Further, for each fixed M and n,
C the array elements that are non-zero are:
C P(n,m-n+1,0) -----> 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