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 = 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