PROGRAM trntr2 C Comment lines found in trntr1 are not repeated here, and should be read C first. This program modifies trntr1 only in the output of the C coefficients for the resistance functions RYA, RYB, and RYG. C The results are written out in a form suitable for the subroutines C RYA, RAY, RYG, RGY, RBY, RYB. C To calculate large numbers of terms, the memory requirements are: C C MAXS MMAX DIMENSION P,Q or V MEMORY (for one array) FAC1 C ---- ---- ------------------- ------ (per 8 bytes) ---- 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.5MB 76,76 C 200 201 692352 5.5MB 101,101 C 250 251 1342002 10.75MB 126,126 C 300 301 2307277 18.5MB 151,151 C DOUBLE PRECISION P(0:2307277),V(0:2307277),Q(0:2307277) DOUBLE PRECISION FAC1(151,151),XN DOUBLE PRECISION FCTR1,FCTR2,FCTR3,FCTR4,FCTRQ,FCTRV,FACTOR DOUBLE PRECISION 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 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 for translation. P(INDX1(1,0,0))=1D0 V(INDX1(1,0,0))=1D0 Q(INDX1(1,0,0))=0D0 DO 500 M=1,MMAX NM=M/2+1 IF(M.EQ.MMAX)NM=2 DO 400 N=1,NM 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 IP=M-IQ-N+1 JS=(IQ+1)/2 KS=IQ/2 PSUM=0D0 VSUM=0D0 QSUM=0D0 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=DBLE(IS) FAC=FAC1(N,IS) C FACTOR=2D0*(XN*XS+1D0)**2/(XN+XS) PSUM=PSUM+FAC*FCTR1*(4D0+XN*XS-FACTOR) & *P(INDX1(IS,IQ-IS,IP-N+1))/(XS*(2D0*XS-1D0)) 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)) 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+1D0) ENDIF ENDIF IF(IS.LT.JS) THEN PSUM=PSUM+FAC*FCTR3*V(INDX1(IS,IQ-IS-2,IP-N+1))/(2D0*XS+1D0) 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 C I map the RYA coefficients using: P(INDX1(1,m-q,q)) -> V( 1+m(m+1)/2+q ). C I map the RYB ones using: 2Q(INDX1(1,m-q,q)) -> Q( 1+m(m+1)/2+q ). C DO 600 M=1,MAXS NM=(M*(M+1))/2 DO 600 IQ=0,M Q(NM+IQ)=2D0*Q(INDX1(1,M-IQ,IQ)) 600 V(NM+IQ)=P(INDX1(1,M-IQ,IQ)) C C I map the RYG ones using: (3/4) P(INDX1(2,m-q,q)) -> P( 1+m(m+1)/2+q ). C The term P(INDX1(2,0,1)) was not stored and is set separately. C DO 620 M=1,MAXS NM = (M*(M+1))/2 DO 610 IQ=0,(M-1) 610 P(NM+IQ)=0.75D0*P( INDX1(2,M-IQ,IQ) ) 620 P(NM+M) = 0D0 C C Now that the calculations are complete, I shall change the first element C of the output to a signature value. Later when I sum the series C using standard subroutines, I shall use this element to check that C the correct coefficients have been read in. The real value of the C element will be supplied by the subroutines. C C The signature of RYA is -2, of RYB -3, and of RYG -7. C P(0)=-7D0 Q(0)=-3D0 V(0)=-2D0 C MMAX=((MAXS+1)*(MAXS+2))/2 -1 700 FORMAT(D22.16) OPEN (1,FILE='ryacof.dat',STATUS='NEW') WRITE(1,705) MAXS 705 FORMAT(' RYA',I5) DO 710 M=0,MMAX 710 WRITE (1,700) V(M) CLOSE (1) C OPEN (2,FILE='rybcof.dat',STATUS='NEW') WRITE(2,715) MAXS 715 FORMAT(' RYB',I5) DO 720 M=0,MMAX 720 WRITE (2,700) Q(M) CLOSE (2) C OPEN (3,FILE='rygcof.dat',STATUS='NEW') WRITE(3,725) MAXS 725 FORMAT(' RYG',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. C900 FORMAT(D27.21) C OPEN (4,FILE='ryaeqcof.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) CC C OPEN (10,FILE='rybeqcof.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) CC C OPEN (11,FILE='rygeqcof.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) 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