PROGRAM axist2 C Comment lines found in axist1 are not repeated here, and should be read C first. This program modifies axist1 only in the output of the C coefficients for the resistance functions RXM and RXG. The results C are written out in form suitable for the subroutines RXG, RGX, RXM C and RMX. C To calculate large numbers of terms, the memory requirements are: C C MAXS MMAX DIMENSION P OR V MEMORY (PER 8 BYTES) C0 C ---- ---- ----------------- ------ (for each array) -- C 15 17 558 5KB 9,9 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 C The dimensions given below are for MAXS=300. DOUBLE PRECISION P(0:2330231),V(0:2330231),C0(151,151),XN,XS DOUBLE PRECISION FCTR1,FCTR2,FCTR3,FCTRV,FAC,PSUM,VSUM INTEGER MMAX,MAXS,NM,N,IS,M,NQ,IQ,INDX2,JS,IP LOGICAL EVEN COMMON/ ARRAYS/P,V,C0 C This COMMON block reduces disk space requirements with many compilers C and helps memory management on some computers. MAXS=300 MMAX=MAXS+2 C First tabulate C0 DO 10 N=1,(MMAX+1)/2 XN=DBLE(N) C0(N,1)=1D0 DO 10 IS=2,(MMAX+1)/2 XS=DBLE(IS) C0(N,IS)=C0(N,IS-1)*(XN+XS)/XS 10 CONTINUE C We set the initial conditions. The labelling scheme needs some zero C elements at the start. DO 5 M=0,5 P(M)=0D0 5 V(M)=0D0 P(INDX2(2,0,0))=1D0 V(INDX2(2,0,0))=1D0 C We start at N=1, Q=0, P=M-1 and proceed keeping M fixed. C NM is the solution of the equation M=N+P+Q=NM+NM-2+0, except C the last time through when only N=1 and 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=N+N-2+NQ NQ=M-2*(N-1) XN=DBLE(N) FCTR1=XN*(XN+.5D0) FCTR2=XN*(XN-.5D0) FCTR3=XN*(2D0*XN*XN-.5D0) FCTRV=2D0*XN/(2D0*XN+3D0) 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-2=IQ-JS. IP=M-IQ-N PSUM=0D0 VSUM=0D0 C IF ( IQ.EQ.0 .OR. IQ.EQ.NQ) GO TO 250 JS=IQ/2+1 DO 200 IS=1,JS XS=DBLE(IS) FAC=C0(N,IS) C PSUM=PSUM+FAC*FCTR1*(2D0*XN*XS-XN-XS+2D0)* & P(INDX2(IS,IQ-IS,IP-N+1))/((XN+XS)*(2D0*XS-1D0)) IF(IP.GT.N) THEN PSUM=PSUM-FAC*FCTR2*P(INDX2(IS,IQ-IS,IP-N-1)) VSUM=VSUM-FAC*FCTRV*P(INDX2(IS,IQ-IS,IP-N-1)) ENDIF C 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 P(INDX2(N,IP,IQ))=PSUM V(INDX2(N,IP,IQ))=PSUM+VSUM 300 CONTINUE 400 CONTINUE 500 CONTINUE C Collect the coefficients for writing. DO 610 M=0,MAXS EVEN= M .EQ. 2*(M/2) NM=(M*(M+1))/2 DO 610 IQ=0,M IF (EVEN) THEN V(NM+IQ)=5D0*P(INDX2(1,IQ-1,M-IQ+1)) P(NM+IQ)=P(INDX2(2,M-IQ,IQ)) ELSE V(NM+IQ)=5D0*P(INDX2(1,M-IQ,IQ)) IF (IQ.EQ.0) THEN P(NM+IQ) = 0D0 ELSE P(NM+IQ)=P(INDX2(2,M-IQ+1,IQ-1)) ENDIF ENDIF 610 CONTINUE C V(0)=-6D0 P(0)=-9D0 MMAX=((MAXS+1)*(MAXS+2))/2 -1 OPEN (1,FILE='rxgcof.dat',STATUS='NEW') WRITE(1,650) MAXS 650 FORMAT(' RXG',I5) C We do not use an implied DO loop because the computer (ETA-10) crashes. C The FORMAT may need a leading space on some compilers. DO 690 M=0,MMAX 690 WRITE (1,700)V(M) 700 FORMAT(D22.16) CLOSE (1) C OPEN (2,FILE='rxmcof.dat',STATUS='NEW') WRITE(2,705) MAXS 705 FORMAT(' RXM', I5) DO 710 M=0,MMAX 710 WRITE (2,700)P(M) CLOSE (2) C C Previously the equal spheres case was written out separatly using the C commented out code below. C900 FORMAT(D27.21) C OPEN (3,FILE='rxgeqdof.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(3,900)VSUM C940 CONTINUE C CLOSE (3) C C OPEN (4,FILE='rxmeqcof.dat') C DO 990 M=0,MAXS C PSUM=0D0 C NM=((M+1)*M)/2 C DO 970 IQ=0,M C PSUM=PSUM+P(NM+IQ) C970 CONTINUE C PSUM=PSUM*2D0**(-M) C WRITE(4,900)PSUM C990 CONTINUE C CLOSE (4) STOP END C INTEGER FUNCTION INDX2(N,IP,IQ) C ------------------------ Arguments --------------------- INTEGER N, IP, IQ 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 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