PROGRAM trnst2
C It calculates coefficients for the resistance functions
C
C G H M
C Y , Y , and Y . 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:12856),V(0:12856),Q(0:12856),FAC1 (26,26)
DOUBLE PRECISION FCTR1,FCTR2,FCTR3,FCTR4,FCTRQ,FCTRV,FACTOR
DOUBLE PRECISION PSUM,VSUM,QSUM,FAC,XN,XS
INTEGER MMAX,MAXS,NM,N,IS,M,NQ,IQ,INDX2,JS,IP,KS,ISUB
LOGICAL EVEN
COMMON/ARRAYS/P,V,Q,FAC1
MAXS=50
MMAX=MAXS+2
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
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)/(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
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=FAC1(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. 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(INDX2(IS,IQ-IS,IP-N+1))/(XS*(2D0*XS-1D0))
IF(IP.GE.N) THEN
QSUM=QSUM-FAC*FCTRQ*P(INDX2(IS,IQ-IS,IP-N))/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))/(XN+1D0)
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
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)=-10D0*Q( INDX2(1,M-IQ,IQ) )/9D0
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)=-10D0*Q( INDX2(1,IQ-2,M-IQ+2) )/9D0
ELSE
Q(NM+IQ)=0D0
ENDIF
P(NM+IQ)=P( INDX2(2,M-IQ+1,IQ-1) )
ENDIF
610 CONTINUE
C
V(0)=-7D0
Q(0)=-8D0
P(0)=-10D0
C
MMAX=((MAXS+1)*(MAXS+2))/2 -1
700 FORMAT(D22.16)
OPEN (1,FILE='rygcof.dat',STATUS='NEW')
WRITE(1,705) MAXS
705 FORMAT(' RYG',I5)
DO 710 M=0,MMAX
710 WRITE (1,700) V(M)
C
C THE CYBERS AND ETA DO NOT NEED A SPACE IN FORMAT STATEMENTS FOR
C DISK FILES.
C
CLOSE (1)
C
OPEN (2,FILE='ryhcof.dat',STATUS='NEW')
WRITE(2,715) MAXS
715 FORMAT(' RYH',I5)
DO 720 M=0,MMAX
720 WRITE (2,700) Q(M)
CLOSE (2)
C
OPEN (3,FILE='rymcof.dat',STATUS='NEW')
WRITE(3,725) MAXS
725 FORMAT(' RYM',I5)
DO 730 M=0,MMAX
730 WRITE (3,700) P(M)
CLOSE (3)
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 (4,FILE='rygeqcof.dat')
C DO 940 M=0,MAXS
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(4,900)VSUM
C940 CONTINUE
C CLOSE (4)
C OPEN (10,FILE='ryheqcof.dat')
C DO 990 M=0,MAXS
C QSUM=0D0
C NM=((M+1)*M)/2
C DO 970 IQ=0,M
C QSUM=QSUM+Q(NM+IQ)
C970 CONTINUE
C QSUM=QSUM*2D0**(-M)
C WRITE(10,900)QSUM
C990 CONTINUE
C CLOSE (10)
C
C OPEN (11,FILE='rymeqcof.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)
C
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 =<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