PROGRAM axist4 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=50. 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 690 m=1,maxs EVEN= M .EQ. 2*(M/2) NM = (M*(M+1))/2 v(nm)=0d0 if (even) then do 640 iq=1,m-1 ip=m-iq-1 psum=0d0 do 630 n=1,(iq+2)/2 psum=psum+p( indx2(n,iq-n,ip) ) 630 continue 640 v( nm+iq) = 2.5d0*psum v(nm+m)=0d0 else do 660 iq=1,m ip=m-iq psum=0d0 do 650 n=1,(iq+2)/2 psum=psum+p( indx2(n,iq-n-1,ip) ) 650 continue 660 v( nm+iq) = 2.5d0*psum endif 690 continue C V(0)=-13D0 MMAX=((MAXS+1)*(MAXS+2))/2 -1 OPEN (1,FILE='rxqcof.dat',STATUS='NEW') WRITE(1,710) MAXS 710 FORMAT(' RXQ',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 790 M=0,MMAX 790 WRITE (1,700)V(M) 700 FORMAT(D22.16) CLOSE (1) C Previously the equal spheres case was written out separatly using C the commented out code below. C900 FORMAT(D27.21) C OPEN (3,FILE='rxqeqcof.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 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