PROGRAM AXITR2 C Comment lines found in AXITR1 are not repeated here, and should be read C first. This program modifies AXITR1 only in the output of the C coefficients for the resistance functions RXA and RXG. C The results are written out in a form suitable for the subroutines C RXA, RAX, RXG and RGX. C To calculate large numbers of terms, the memory requirements are: C C MAXS MMAX DIMENSION P OR V MEMORY (PER 8 BYTES) FAC0 C ---- ---- ---------------- ------ (for each array) ---- C 15 16 475 4KB 8,8 C 50 51 12152 97KB 26,26 C 100 101 89927 700KB 51,51 C 150 151 295827 2.4KB 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 C The dimensions given below are for MAXS = 300. DOUBLE PRECISION P(0:2307277),V(0:2307277),FAC0(151,151),XN,XS DOUBLE PRECISION FCTR1,FCTR2,FCTR3,FCTRV,FAC,PSUM,VSUM INTEGER MMAX,MAXS,NM,N,IS,M,NQ,IQ,INDX1,JS,IP COMMON /ARRAYS/P,V,FAC0 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 FAC0 DO 10 N=1,(MMAX+1)/2 XN=DBLE(N) FAC0(N,1)=1D0+XN DO 10 IS=2, (MMAX+1)/2 XS=DBLE(IS) FAC0(N,IS)=FAC0(N,IS-1)*(XN+XS)/XS 10 CONTINUE C We set the initial conditions. P(INDX1(1,0,0))=1D0 V(INDX1(1,0,0))=1D0 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=DBLE(N) FCTR1=XN*(XN+.5D0)/(XN+1D0) FCTR2=XN*(XN-.5D0)/(XN+1D0) FCTR3=XN*(2D0*XN*XN-.5D0)/(XN+1D0) FCTRV=2D0*XN/((XN+1D0)*(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-1=IQ-JS. IP=M-IQ-N+1 PSUM=0D0 VSUM=0D0 IF (IQ.EQ.0D0) GO TO 250 JS=(IQ+1)/2 DO 200 IS=1,JS XS=DBLE(IS) FAC=FAC0(N,IS) PSUM=PSUM+FAC*FCTR1*(2D0*XN*XS-XN-XS+2D0)* & P(INDX1(IS,IQ-IS,IP-N+1))/((XN+XS)*(2D0*XS-1D0)) 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.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 CONTINUE 400 CONTINUE 500 CONTINUE C Collect the coefficients for writing. C DO 600 M=0,MAXS NM=(M*(M+1))/2 DO 600 IQ=0,M V(NM+IQ)=P(INDX1(1,M-IQ,IQ)) 600 CONTINUE P(0)=0D0 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) ) P(NM+M) =0D0 620 CONTINUE C MMAX=((MAXS+1)*(MAXS+2))/2 -1 OPEN (1,FILE='rxacof.dat',STATUS='NEW') C First we write an identification line WRITE(1,650) MAXS 650 FORMAT(' RXA',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='rxgcof.dat',STATUS='NEW') WRITE(2,710) MAXS 710 FORMAT(' RXG',I5) DO 720 M=0,MMAX 720 WRITE (2,700)P(M) CLOSE (2) C Previously the equal spheres case was written out separately using C the commented out code below. C OPEN (3,FILE='rxaeqcof.dat',STATUS='NEW') 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,700)VSUM C940 CONTINUE C CLOSE(3) C OPEN (4,FILE='rxgeqcof.dat',STATUS='NEW') 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,700)PSUM C990 CONTINUE C CLOSE(4) 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