PROGRAM trnro3
C Comment lines found in trnro1 are not repeated here, and should be read
C first.
C
C ------------------------- WARNING -----------------------------
C This program uses quadruple precision (REAL*16) which is only supported
C by a few compilers. It was the program actually used to obtain the
C data files.
C
C To calculate large numbers of terms, the memory requirements are:
C
C MAXS MMAX DIMENSION P,V OR Q MEMORY (PER 8 BYTES) FAC0
C ---- ---- ---------------- ------ (for each array) ----
C 15 16 475 5KB 8,8
C 50 51 12152 100KB 26,26
C 100 101 89927 700KB 51,51
C 150 151 295827 2.4MB 76,76
C 200 201 692352 5.5MB 101,101
C 250 251 1342002 10.7MB 126,126
C 300 301 2307277 18.5MB 151,151
C
REAL*16 P(0:2307277),V(0:2307277),Q(0:2307277),FAC1(151,151),XN
REAL*16 FCTR1,FCTR2,FCTR3,FCTR4,FCTRQ,FCTRV,FACTOR
REAL*16 PSUM,VSUM,QSUM,FAC,XS
INTEGER MMAX,MAXS,NM,N,IS,M,NQ,IQ,INDX1,JS,IP,KS
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=300
MMAX=MAXS+1
C First tabulate FAC1
DO 10 N=1,(MMAX+1)/2
XN=qext(N)
FAC1(N,1)=1q0
DO 10 IS=2,(MMAX+1)/2
XS=qext(IS)
10 FAC1(N,IS)=FAC1(N,IS-1)*(XN+XS)/(XS-1q0)
C We first set the initial conditions for rotation.
P(INDX1(1,0,0))=0q0
V(INDX1(1,0,0))=0q0
Q(INDX1(1,0,0))=1q0
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=qext(N)
FCTR1=(XN+.5q0)/(XN+1q0)
FCTR2=XN*(XN-.5q0)/(XN+1q0)
FCTR3=XN*(2q0*XN*XN-.5q0)/(XN+1q0)
FCTR4=(8q0*XN*XN-2q0)/(3q0*XN+3q0)
FCTRV=2q0*XN/((XN+1q0)*(2q0*XN+3q0))
FCTRQ=3q0/(2q0*XN*(XN+1q0))
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=0q0
VSUM=0q0
QSUM=0q0
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=qext(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=2q0*(XN*XS+1q0)**2/(XN+XS)
PSUM=PSUM+FAC*FCTR1*(4q0+XN*XS-FACTOR)
& *P(INDX1(IS,IQ-IS,IP-N+1))/(XS*(2q0*XS-1q0))
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+1q0)
ENDIF
ENDIF
IF(IS.LT.JS)THEN
PSUM=PSUM+FAC*FCTR3*V(INDX1(IS,IQ-IS-2,IP-N+1))/(2q0*XS+1q0)
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
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
DO 600 M=0,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.5q0*P(INDX1(1,M-IQ,IQ))
P(0)=0q0
DO 620 M=1,MAXS
NM=(M*(M+1))/2
DO 610 IQ=0,M-1
610 P(NM+IQ)=-0.375q0*P(INDX1(2,M-IQ,IQ))
620 P(NM+M) = 0q0
DO 660 M=1,MAXS,2
NM=(M*(M+1))/2
DO 650 NQ=M,1,-1
Q(NM+NQ)=Q(NM+NQ-1)
650 P(NM+NQ)=P(NM+NQ-1)
Q(NM)=0q0
660 P(NM)=0q0
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
MMAX=((MAXS+1)*(MAXS+2))/2 -1
700 FORMAT(D22.16)
OPEN (1,FILE='rybcof.dat',STATUS='NEW')
WRITE(1,710)MAXS
710 FORMAT(' RYB',I5)
DO 720 M=0,MMAX
720 WRITE (1,700) V(M)
CLOSE (1)
C
OPEN (2,FILE='ryccof.dat',STATUS='NEW')
WRITE(2,740)MAXS
740 FORMAT(' RYC',I5)
DO 750 M=0,MMAX
750 WRITE (2,700) Q(M)
CLOSE (2)
C
OPEN (3,FILE='ryhcof.dat',STATUS='NEW')
WRITE(3,760)MAXS
760 FORMAT(' RYH',I5)
DO 770 M=0,MMAX
770 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='rybeqcof.dat',STATUS='NEW')
C DO 940 M=0,MAXS
C VSUM=0q0
C NM=((M+1)*M)/2
C DO 930 IQ=0,M
C VSUM=VSUM+V(NM+IQ)
C930 CONTINUE
C VSUM=VSUM*2q0**(-M)
C WRITE(4,900)VSUM
C940 CONTINUE
C CLOSE (4)
C OPEN (10,FILE='ryceqcof.dat',STATUS='NEW')
C DO 990 M=0,MAXS
C QSUM=0q0
C NM=((M+1)*M)/2
C DO 970 IQ=0,M
C QSUM=QSUM+Q(NM+IQ)
C970 CONTINUE
C QSUM=QSUM*2q0**(-M)
C WRITE(10,900)QSUM
C990 CONTINUE
C CLOSE (10)
C OPEN (11,FILE='ryheqcof.dat',STATUS='NEW')
C DO 1090 M=0,MAXS
C PSUM=0q0
C NM=((M+1)*M)/2
C DO 1070 IQ=0,M
C PSUM=PSUM+P(NM+IQ)
C1070 CONTINUE
C PSUM=PSUM*2q0**(-M)
C WRITE(11,900)PSUM
C1090 CONTINUE
C CLOSE (11)
STOP
END
C
INTEGER FUNCTION INDX1(N,IP,IQ)
C ------------------------- Arguments ------------------------
INTEGER N,IP,IQ
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