PROGRAM axitr3 C Comment lines found in axitr1 are not repeated here, and should be read C first. C ------------------------ WARNING --------------------------- C This program uses quadruple precision (REAL*16) which is only supported C by a few compilers. It was the program actually used to obtain the C data files. C C The dimensions given below are for MAXS = 50. REAL*16 P(0:2307277),V(0:2307277),FAC0(151,151),XN,XS REAL*16 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=qext(N) FAC0(N,1)=1q0+XN DO 10 IS=2, (MMAX+1)/2 XS=qext(IS) FAC0(N,IS)=FAC0(N,IS-1)*(XN+XS)/XS 10 CONTINUE C We set the initial conditions. P(INDX1(1,0,0))=1q0 V(INDX1(1,0,0))=1q0 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 Write(*,*)m 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=qext(N) FCTR1=XN*(XN+.5q0)/(XN+1q0) FCTR2=XN*(XN-.5q0)/(XN+1q0) FCTR3=XN*(2q0*XN*XN-.5q0)/(XN+1q0) FCTRV=2q0*XN/((XN+1q0)*(2q0*XN+3q0)) 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=0q0 VSUM=0q0 IF (IQ.EQ.0q0) GO TO 250 JS=(IQ+1)/2 DO 200 IS=1,JS XS=qext(IS) FAC=FAC0(N,IS) PSUM=PSUM+FAC*FCTR1*(2q0*XN*XS-XN-XS+2q0)* & P(INDX1(IS,IQ-IS,IP-N+1))/((XN+XS)*(2q0*XS-1q0)) 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))/(2q0*XS+1q0) 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)=0q0 DO 620 M=1,MAXS NM = (M*(M+1))/2 DO 610 IQ=0, M-1 610 P(NM+IQ)=0.75q0*P( INDX1(2,M-IQ,IQ) ) P(NM+M) =0q0 620 CONTINUE C MMAX=((MAXS+1)*(MAXS+2))/2 -1 OPEN (1,FILE='rxacof.dat') 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 one computer (ETA-10) crashed. 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') 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') C DO 940 M=0,MAXS C VSUM=0q0 C NM=((M+1)*M)/2 C DO 930 IQ=0,M C VSUM=VSUM+V(NM+IQ) C930 CONTINUE C VSUM=VSUM*2q0**(-M) C WRITE(3,700)VSUM C940 CONTINUE C CLOSE(3) C OPEN (4,FILE='rxgeqcof.dat') C DO 990 M=0,MAXS C PSUM=0q0 C NM=((M+1)*M)/2 C DO 970 IQ=0,M C PSUM=PSUM+P(NM+IQ) C970 CONTINUE C PSUM=PSUM*2q0**(-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