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