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