############################################################################## # # # Maple V Programs for Two-Center STO Electron Repulsion Integrals # # # # Frank E. Harris # # # # Department of Physics, University of Utah and Quantum Theory Project, # # University of Florida # # # # Adjunct to paper for International Journal of Quantum Chemistry entitled # # "Analytical Evaluation of Two-Center STO Electron Repulsion Integrals # # via Ellipsoidal Expansion", submitted September, 2001. # # # # (Last modified 23 September 2001) # # # # Author's e-mail address: harris@qtp.ufl.edu # # # ############################################################################## # These programs have been developed under Maple V, Release 5, and will # probably also run on newer releases. They are not intended for # production use, as no attempt has been made to optimize them either for # speed or for accuracy. Their main purposes are to validate the methods # described in the paper identified above, to enable studies of the # numerical stability of these and other methods, and to provide tools for # checking production programs. Equation numbers referred to in the # various procedures are those of the IJQC paper. # # This file contains programs that can be used to compute individual two- # electron integrals and their constituent parts by various methods. It # is in a format permitting them to be processed (either under command-line # Maple or a Maple worksheet) by cut-and-paste methods or by use of the # following Maple command in either environment: # # read "name_you_have_given_the_file_(inside_double-quotes)"; # # In some PC environments it will necessary to give a full path # specification, replacing the PC-world's "\" characters by "/". # #----------------------------------------------------------------------------- # # Program 1: Reference closed-form computation of a two-electron STO # integral, by formula in Eq. (10). # # This program is painfully slow and is included mainly to define everything. # It should be used (if at all) with an appropriate setting of the Maple # precision variable "Digits". More efficient computation can be achieved # using other programs of this package. # # Evaluates integral psi1*(1)psi2(1) (1/r12) psi3*(2)psi4(2) dtau1 dtau2. # R is intercenter distance. Each orbital is specified by its values of # zeta, n, l, and m. The signs of the parameters zeta identify the center # on which the STO's are situated (positive zeta --> centered on A, # negative zeta --> centered on B). Calculations are, of course, for # |zeta|. The parameter mumax is the largest value retained in the mu # summation (that resulting from the Neumann expansion). The positive z # axis for each orbital is in the direction toward the other center, the # phi dependence of the orbitals is of the form exp(i m phi), with the phi # coordinate the same for both centers (i.e. one has a right-handed # coordinate system, while the other is left-handed), the associated # Legendre functions are as defined in the Handbook of Mathematical # Functions, and the STO's are normalized. Notice that the ordering of # the orbitals corresponds to Mulliken notation, so that the output is, in # that notation, the two-electron integral [12|34]. Hartree atomic units # are used: zeta is in inverse bohr, R is in bohr, [12|34] is in hartree. # # Usage: Two_el_0(R,zeta1,n1,l1,m1,zeta2,n2,l2,m2,zeta3,n3,l3,m3,zeta4,n4, # l4,m4,mumax); # #----------------------------------------------------------------------------- # # Program 2. Closed-form computation of two-electron STO integral, using # precalculated arrays # # A suitable setting of the Maple precision variable "Digits" is # recommended. # # Usage: Two_el_1(R,zeta1,n1,l1,m1,zeta2,n2,l2,m2,zeta3,n3,l3,m3,zeta4,n4, # l4,m4,mumax); # #----------------------------------------------------------------------------- # # Program 3. Computation of two-electron STO integral for small alpha # (uses precalculated arrays), calculating L and W using # Eqs. (41), (44), and (54). # # Usage: Two_el_2(R,zeta1,n1,l1,m1,zeta2,n2,l2,m2,zeta3,n3,l3,m3,zeta4,n4, # l4,m4,mumax,tmax,jmax); # # The new argument tmax is the maximum t value retained in the # infinite sum of Eq. (44); jmax is the maximum j retained in Eq. # (54). # #============================================================================= # # Overall program test. Computes five integrals in each of three ways, at a # very low level of precision. The output includes # text indicating the expected values and relevant # comments. This is preimarily included to provide a # way of verifying that the program suite works. # This test is supplied as a separate program. # # Usage: read "your-name-for-the-test-program"; # # #============================================================================= # # Test Programs--nonrecursive processes # #----------------------------------------------------------------------------- # # Test 1. Calculate an individual L function, Eq. (28), directly from # the definition by numerical integration, for parameters mu, # sigma, p, and argument alpha. Caution: SLOW! # # Usage: LLL(mu,sigma,p,alpha); # #----------------------------------------------------------------------------- # # Test 2. Calculate array of functions L, Eq. (28), by closed formula. # Array LX[mu,p] contains a set of L for mu from sigma to mumax # and p from 0 to N1 for a given value of sigma and a given alpha. # # Usage: Build_LX(LX,mumax,sigma,N1,alpha); # #----------------------------------------------------------------------------- # # Test 3. Calculate array of L by formula for small alpha. Array LY[mu,p] # contains a set of L. Arguments as in Test 1, plus tmax (the # maximum t in Eq. (44)). # # Usage: Build_LY(LY,mumax,sigma,pmax,alpha,tmax); # #----------------------------------------------------------------------------- # # Test 4. Calculate array of W, Eq. (12), by closed formula. Output # array WX[mu,p1,p2] contains a set of W for mu from sigma to mumax, # p1 from 0 to N1, and p2 from 0 to N2, for given sigma, alpha1, # and alpha2. # # Usage: Build_WX(WX,mumax,sigma,N1,N2,alpha1,alpha2); # #----------------------------------------------------------------------------- # # Test 5. Calculate array of W by formula for small alpha, Eq. (54). Output # WY[mu,p1,p2] contains a set of W. Arguments as in Test3, plus # tmax (cutoff in Eq. (44)) and jmax (cutoff in Eq. (54)). # # Usage: Build_WY(WY,mumax,sigma,N1,N2,alpha1,alpha2,tmax,jmax); # #----------------------------------------------------------------------------- # # Test 6. Calculate array of W by alternate formula for small alpha, Eq. # (55). Only for sigma=p1=p2=0. Output WYY[mu] contains a set of # W for mu from 0 to mumax. Cutoffs tmax and jmax as in Test 4, # the p sum in Eq. (55) is cut off at pmax. # # Usage: Build_WYY(WW1,mumax,0,0,0,alpha1,alpha2,tmax,jmax,pmax); # #----------------------------------------------------------------------------- # # Test 7. Comparison of LX and LY at current precision (cf. Tests 2 and 3) # with exact values for input values of sigma, p, and alpha, for # mu from sigma to mumax. The exact values are generated by closed # formula at precision 60S higher than the current value. The # output is a table, sent to the current standard output device (the # terminal if alternate arrangements have not been made). Argument # tmax is the summation cutoff for LY. # # Usage: Test_L(mumax,sigma,p,alpha,tmax); # #----------------------------------------------------------------------------- # # Test 8. Comparison of WX and WY at current precision (cf. Tests 4 and 5) # with exact values for input values of sigma, p1, p2, alpha1, and # alpha2, for mu from sigma to mumax. Arguments tmax and jmax are # the summation cutoffs for WY. Tabular output. # # Usage: Test_W(mumax,sigma,p1,p2,alpha1,alpha2,tmax,jmax); # #----------------------------------------------------------------------------- # # Test 9. Comparison of WYY at current precision (cf. Test 6) with exact # values for input sigma=p1=p2=0, alpha1, and alpha2, for mu from 0 # to mumax. Arguments tmax, jmax, and pmax are summation cutoffs # for WYY. Tabular output. # # Usage: Test_WYY(mumax,alpha1,alpha2,tmax,jmax,pmax); # #============================================================================= # # Test Programs--recursive processes # #----------------------------------------------------------------------------- # # Test 10. Tests recursion in mu for L, comparing upward recursion using # Eq. (69) with downward recursion using that equation and with # downward recursion followed by the adjustment described by # Eqs. (74) and (75), for sigma=p=0, and argument alpha. Range # of mu is 0 to mumax. Tabular output. # # Usage: Test_L_recur(mumax,alpha); # #----------------------------------------------------------------------------- # # Test 11. Tests recursion in p for L, starting from L[mu,p] with sigma=0 # for the necessary range of mu values. Compares "vertical" # recurrence, using Eqs. (65)-(68) for mu=0 and Eqs. (70)-(73) # for nonzero mu, with "oblique" recurrence, using Eq. (60) for # all mu values except 0 and mumax, the maximum value considered. # Runs recursion in p from zero to pmax for argument alpha, # tabulates output for mu from 0 to mumax and for p=pmax only. # # Usage: Test_Lp_recur(mumax,pmax,alpha); # #----------------------------------------------------------------------------- # # Test 12. Tests recursion in sigma for L, starting from L[mu,p] for the # input p value, for sigma=0, and for the necessary range of mu # values. Runs recursion for argument alpha until the input # sigma value is reached, tabulates output for this p and sigma # and for mu from sigma to mumax. # # Usage: Test_Ls_recur(mumax,sigma,p,alpha); # #----------------------------------------------------------------------------- # # Test 13. Tests recursion in mu for W, comparing Kotani's method, Eqs. (81)- # (86), with exact values and with the results of upward and # downward recursion using Eqs. (96) and (97), and with downward # recursion adjusted (repeatedly if necessary) using Eqs. (109) and # (110). Obtains and tabulates results for mu from 0 to mumax, # with sigma=p1=p2=0. # # Usage: Test_W_recur(mumax,alpha1,alpha2); # #----------------------------------------------------------------------------- # # Test 14. Tests recursion in p1 and p2 for W, using Eqs. (103), (106), and # (108). Uses adjusted downward recursion for p1=p2=0 and sigma=0, # then recurs in p through the input p1 and p2. Obtains and # tabulates results for mu from 0 to mumax. # # Usage: Test_Recur_Wp(mumax,p1,p2,alpha1,alpha2); # #----------------------------------------------------------------------------- # # Test 15. Tests recursion in sigma for W, using Eq. (77). Uses result of # recursion in p1 and p2 as starting values for the sigma recursion. # Obtains and tabulates results for mu from sigma to mumax for the # input values of sigma, p1, and p2. # # Usage: Test_Recur_Ws(mumax,sigma,p1,p2,alpha1,alpha2); # #----------------------------------------------------------------------------- # # Test 16. Compares accurate values of auxiliary function T with those # evaluated according to Eqs. (157)-(160). # # Usage: Check_TX(mumax,alpha); # #============================================================================= # # Additional tests can be carried out by using the utility programs, which # are self-documented. # #============================================================================= ###################################### # Reference closed-form procedures # ###################################### #----------------------------------------------------------------------------- # Master subroutine for two-electron integral, Eq. (10): # Integral psi1*(1)psi2(1) (1/r12) psi3*(2)psi4(2) dtau1 dtau2: # R is intercenter distance, negative zeta identify orbitals on Center B. # mumax is maximum of mu summation. # KK is coefficient defined after Eq. (7); C1 and C2 (produced by Make_C) # are coefficients of powers of xi and eta in wavefunction products, as # defined in Eqs. (7) and (8); Mu_Term is procedure generating the mu term # of the Neumann expansion Two_el_0 := proc(R,zeta1,n1,l1,m1,zeta2,n2,l2,m2,zeta3,n3,l3,m3,zeta4,n4, l4,m4,mumax) local N1,N2,mu,sigma,K1,K2,alpha1,alpha2,beta1,beta2,C1,C2; sigma := m2-m1; if m3-m4 <> sigma then RETURN(0) fi; sigma := abs(sigma); K1 := KK(R,abs(zeta1),n1,l1,m1,abs(zeta2),n2,l2,m2); K2 := KK(R,abs(zeta3),n3,l3,m3,abs(zeta4),n4,l4,m4); alpha1 := (R/2)*(abs(zeta1)+abs(zeta2)); alpha2 := (R/2)*(abs(zeta3)+abs(zeta4)); beta1 := (R/2)*(zeta1+zeta2); beta2 := (R/2)*(zeta3+zeta4); N1 := n1+n2-sigma; N2 := n3+n4-sigma; C1 := array(0..N1,0..N1); C2 := array(0..N2,0..N2); Make_C(C1,signum(zeta1),n1,l1,m1,signum(zeta2),n2,l2,m2); Make_C(C2,signum(zeta3),n3,l3,m3,signum(zeta4),n4,l4,m4); (4./R)*K1*K2*sum('Mu_Term(mu,sigma,N1,N2,C1,C2,alpha1,alpha2,beta1,beta2)', mu=sigma..mumax); end: printf("Two_el_0 loaded--original\n"); #----------------------------------------------------------------------------- # Normalization factor--coefficient defined after Eq. (7): KK := proc(R,zeta1,n1,l1,m1,zeta2,n2,l2,m2); evalf((1/2)*(R*zeta1)^(n1+1/2)*(R*zeta2)^(n2+1/2)* ((2*l1+1)*(2*l2+1)*(l1-m1)!*(l2-m2)! /((2*n1)!*(2*n2)!*(l1+m1)!*(l2+m2)!))^(1/2)); end: printf("Normalization factor KK loaded\n"); #----------------------------------------------------------------------------- # Make mu term of angular sum in Neumann formula. # II is generalized Bessel function, Eq. (11). # WW is w of Eq. (13). Mu_Term := proc(mu,sigma,N1,N2,C1,C2,alpha1,alpha2,beta1,beta2) local S,CX1,CX2,p1,p2,q; CX1 := array(0..N1); for p1 from 0 to N1 do CX1[p1] := sum('C1[p1,q]*II(mu,sigma,q,beta1)',q=0..N1) od; CX2 := array(0..N2); for p2 from 0 to N2 do CX2[p2] := sum('C2[p2,q]*II(mu,sigma,q,beta2)',q=0..N2) od; S := 0.; for p1 from 0 to N1 do if CX1[p1]=0 then next fi; for p2 from 0 to N2 do if CX2[p2]=0 then next fi; S := S + CX1[p1]*CX2[p2]* (WW(mu,sigma,p1,p2,alpha1,alpha2)+WW(mu,sigma,p2,p1,alpha2,alpha1)) od od; (-1)^sigma*2*(2*mu+1)*S end: printf("Mu_Term loaded--original\n"); #----------------------------------------------------------------------------- # Make coefficients of powers of xi and eta for wavefunction products. # CC[p,q] contains coeff of xi^p eta^q, as defined in Eqs. (7) and (8). # Note that orthopoly package defines P(n,x) as Legendre polynomial. # Load_Coeffs moves the coefficients of the expression in S1 into the # array CC. s1 or s2 are +1 for orbital on A, -1 for orbital on B. Make_C := proc(CC,s1,n1,l1,m1,s2,n2,l2,m2) local N1,S1,S2,x,xi,eta; with(orthopoly); if m1=0 then S1 := P(l1,x) else S1:=(-1)^abs(m1)*diff(P(l1,x),x$abs(m1)) fi; if m2=0 then S2 := P(l2,x) else S2:=(-1)^abs(m2)*diff(P(l2,x),x$abs(m2)) fi; S1 := (xi+s1*eta)^(n1-abs(m1)-1)*subs(x=(s1+xi*eta)/(xi+s1*eta),S1) *(xi+s2*eta)^(n2-abs(m2)-1)*subs(x=(s2+xi*eta)/(xi+s2*eta),S2) *(xi^2-eta^2)*((1-eta^2)*(xi^2-1))^((abs(m1)+abs(m2)-abs(m1-m2))/2); N1 := n1+n2-abs(m1-m2); Load_Coeffs(CC,S1,xi,eta,N1); S1; end: printf("Make_C loaded\n"); #----------------------------------------------------------------------------- # Place coefficients of bivariate polynomial in S1 into array CC. # S is in variables xi and eta, of maximum degree nmax in each variable # CC (a previously defined array of dimensions 0..nmax by 0..nmax) is # loaded with elements CC[p,q] which are coefficients of xi^p eta^q in S. Load_Coeffs := proc(CC,S,xi,eta,nmax) local SS,i,j; SS := array(0..nmax); SS[0] := subs(xi=0,S); for i from 1 to nmax do SS[i] := subs(xi=0,diff(S,xi$i))/i! od; for i from 0 to nmax do CC[i,0] := subs(eta=0,SS[i]); for j from 1 to nmax do CC[i,j] := subs(eta=0,diff(SS[i],eta$j))/j! od od end: printf("Load_Coeffs loaded\n"); #----------------------------------------------------------------------------- # Function w, Eq. (13), calculated by closed formula, Eq. (39). # AC(mu,sigma,s) is coefficient A of Eq. (23). # LX(mu,sigma,p,alpha) is L of Eq. (28). # KX(mu,sigma,p,alpha) is k of Eq. (29). WW := proc(mu,sigma,p1,p2,alpha1,alpha2) local S,j,k,s; S := 0.; for j from 0 to (mu+sigma)/2 do s := mu+sigma-2*j; S := S - sum('AC(mu,sigma,s)*((p2+s)!/(k!*alpha2^(p2+s-k+1))) *LX(mu,sigma,p1+k,alpha1+alpha2)',k=0..p2+s) od; ((mu+sigma)!/(mu-sigma)!)^2 *(S + LX(mu,sigma,p1,alpha1)*KX(mu,sigma,p2,alpha2)) end: printf("WW, closed formula, loaded\n"); #----------------------------------------------------------------------------- # Make function L, Eq. (28), using Eq. (33). # LX0 is L with mu=sigma=0. # AC(mu,sigma,s) is coefficient A of Eq. (26). # BC(mu,sigma,s) is coefficient B of Eq. (27). # AA(n,alpha) is A sub n integral, Eq. (34). LX := proc(mu,sigma,p,alpha) local S,j,smax; if mu+sigma >= 1 then smax := mu+sigma-1; S := sum('BC(mu,sigma,smax-2*j)*AA(p+smax-2*j,alpha)',j=0..smax/2) else S := 0. fi; smax := mu+sigma; S + sum('AC(mu,sigma,smax-2*j)*LX0(p+smax-2*j,alpha)',j=0..smax/2) end: printf("LX, closed formula, loaded\n"); #----------------------------------------------------------------------------- # Function L, Eq. (28) for mu=sigma=0, calculated using Eq. (36). # L_Sum is the summation from that equation. # AA(n,alpha) is A sub n integral, Eq. (34), for argument alpha. LX0 := proc(p,alpha); evalf((1/2)*((-1)^(p+1)*AA(p,-alpha)*Ei(1,2*alpha) + (gamma+log(2*alpha))*AA(p,alpha)+L_Sum(p,alpha))) end: printf("LX0, closed formula, loaded\n"); #----------------------------------------------------------------------------- # Summation in Eq. (36). # LC(p,t) is coefficient L defined in Eq. (37). L_Sum := proc(p,alpha) local t; if p<2 then RETURN(0) fi; (p!*exp(-alpha)/alpha^(p+1))*sum('LC(p,t)*alpha^t/t!',t=1..p-1) end: printf("L_Sum loaded--original\n"); #============================================================================= ######################################## # Array-based closed-form procedures # ######################################## # For additional documentation see the reference forms of these programs. #----------------------------------------------------------------------------- Two_el_1 := proc(R,zeta1,n1,l1,m1,zeta2,n2,l2,m2,zeta3,n3,l3,m3,zeta4,n4, l4,m4,mumax) local N1,N2,mu,sigma,K1,K2,alpha1,alpha2,beta1,beta2,C1,C2, II1,II2,WW1; sigma := m2-m1; if m3-m4 <> sigma then RETURN(0) fi; sigma := abs(sigma); K1 := KK(R,abs(zeta1),n1,l1,m1,abs(zeta2),n2,l2,m2); K2 := KK(R,abs(zeta3),n3,l3,m3,abs(zeta4),n4,l4,m4); alpha1 := (R/2)*(abs(zeta1)+abs(zeta2)); alpha2 := (R/2)*(abs(zeta3)+abs(zeta4)); beta1 := (R/2)*(zeta1+zeta2); beta2 := (R/2)*(zeta3+zeta4); N1 := n1+n2-sigma; N2 := n3+n4-sigma; C1 := array(0..N1,0..N1); C2 := array(0..N2,0..N2); Make_C(C1,signum(zeta1),n1,l1,m1,signum(zeta2),n2,l2,m2); Make_C(C2,signum(zeta3),n3,l3,m3,signum(zeta4),n4,l4,m4); Build_II(II1,mumax,sigma,N1,beta1); Build_II(II2,mumax,sigma,N2,beta2); Build_WX(WW1,mumax,sigma,N1,N2,alpha1,alpha2); (4./R)*K1*K2* sum('Mu_Term1(mu,sigma,N1,N2,C1,C2,II1,II2,WW1,alpha1,alpha2,beta1,beta2)', mu=sigma..mumax); end: printf("Two_el_1 loaded--array version\n"); #----------------------------------------------------------------------------- Mu_Term1 := proc(mu,sigma,N1,N2,C1,C2,II1,II2,WW1,alpha1,alpha2,beta1,beta2) local S,CX1,CX2,p1,p2,q; CX1 := array(0..N1); for p1 from 0 to N1 do CX1[p1] := sum('C1[p1,q]*II1[mu,q]',q=0..N1) od; CX2 := array(0..N2); for p2 from 0 to N2 do CX2[p2] := sum('C2[p2,q]*II2[mu,q]',q=0..N2) od; S := 0.; for p1 from 0 to N1 do if CX1[p1]=0 then next fi; for p2 from 0 to N2 do if CX2[p2]=0 then next fi; S := S + CX1[p1]*CX2[p2]*WW1[mu,p1,p2] od od; (-1)^sigma*2*(2*mu+1)*S end: printf("Mu_Term1 loaded--array version\n"); #----------------------------------------------------------------------------- # Version using preloaded array LC1 L_SumX := proc(p,alpha) local t; global LC1; if p<2 then RETURN(0) fi; (p!*exp(-alpha)/alpha^(p+1))*sum('LC1[p,t]*alpha^t/t!',t=1..p-1) end: printf("L_SumX loaded--array version\n"); #----------------------------------------------------------------------------- # Builds array of W functions--capital W, Eq. (12), using Eq. (39). Build_WX := proc(WW1,mumax,sigma,N1,N2,alpha1,alpha2) local KX1,KX2,LX1,LX2,LX12,mu,p1,p2,j,k,s,S; global AC1; WW1 := array(sigma..mumax,0..N1,0..N2); Build_AC(mumax,sigma); Build_KX(KX1,mumax,sigma,N1,alpha1); Build_KX(KX2,mumax,sigma,N2,alpha2); Build_LX(LX1,mumax,sigma,N1,alpha1); Build_LX(LX2,mumax,sigma,N2,alpha2); Build_LX(LX12,mumax,sigma,N1+N2+2*(mumax+sigma),alpha1+alpha2); for mu from sigma to mumax do for p1 from 0 to N1 do for p2 from 0 to N2 do S := 0.; for k from 0 to (mu+sigma)/2 do s := mu+sigma-2*k; S := S - AC1[mu,sigma,s] *((p2+s)!*sum('LX12[mu,p1+j]/(j!*alpha2^(p2+s-j+1))',j=0..p2+s) +(p1+s)!*sum('LX12[mu,p2+j]/(j!*alpha1^(p1+s-j+1))',j=0..p1+s)) od; WW1[mu,p1,p2] := ((mu+sigma)!/(mu-sigma)!)^2 *(S+LX1[mu,p1]*KX2[mu,p2]+KX1[mu,p1]*LX2[mu,p2]) od od od end: printf("Build_WX loaded\n"); #----------------------------------------------------------------------------- # Builds array of L functions, using Eq. (33) Build_LX := proc(LX1,mumax,sigma,N,alpha) local AA1,LX10,mu,p,k,S; global AC1,BC1; LX1 := array(sigma..mumax,0..N,[]); Build_AC(mumax,sigma); Build_BC(mumax,sigma); Build_AA(AA1,mumax+sigma+N,alpha); Build_LX0(LX10,AA1,mumax+sigma+N,alpha); for mu from sigma to mumax do for p from 0 to N do if mu+sigma>0 then S:=sum('BC1[mu,sigma,mu+sigma-2*k-1]*AA1[p+mu+sigma-2*k-1]', k=0..(mu+sigma-1)/2) else S := 0. fi; LX1[mu,p] := S + sum('AC1[mu,sigma,mu+sigma-2*k]*LX10[p+mu+sigma-2*k]', k=0..(mu+sigma)/2) od od end: printf("Build_LX loaded\n"); #----------------------------------------------------------------------------- # Builds array of L functions for mu=sigma=0, using Eq. (36) Build_LX0 := proc(LX10,AA1,pmax,alpha) local AA2,E1,G1,p; Build_LC(pmax); LX10 := array(0..pmax,[]); Build_AA(AA2,pmax,-alpha); E1 := evalf(Ei(1,2*alpha)); G1 := evalf(gamma+log(2*alpha)); for p from 0 to pmax do LX10[p] := (1/2)*((-1)^(p+1)*E1*AA2[p]+G1*AA1[p]+L_SumX(p,alpha)) od end: printf("Build_LX0 loaded\n"); #============================================================================= ########################### # Small-alpha procedure # ########################### #----------------------------------------------------------------------------- # tmax is summation index in Eq. (44), jmax is that in Eq. (54). Two_el_2 := proc(R,zeta1,n1,l1,m1,zeta2,n2,l2,m2,zeta3,n3,l3,m3,zeta4,n4,l4,m4, mumax,tmax,jmax) local N1,N2,mu,sigma,K1,K2,alpha1,alpha2,beta1,beta2,C1,C2, II1,II2,WW1; sigma := m2-m1; if m3-m4 <> sigma then RETURN(0) fi; sigma := abs(sigma); K1 := KK(R,abs(zeta1),n1,l1,m1,abs(zeta2),n2,l2,m2); K2 := KK(R,abs(zeta3),n3,l3,m3,abs(zeta4),n4,l4,m4); alpha1 := (R/2)*(abs(zeta1)+abs(zeta2)); alpha2 := (R/2)*(abs(zeta3)+abs(zeta4)); beta1 := (R/2)*(zeta1+zeta2); beta2 := (R/2)*(zeta3+zeta4); N1 := n1+n2-sigma; N2 := n3+n4-sigma; C1 := array(0..N1,0..N1); C2 := array(0..N2,0..N2); Make_C(C1,signum(zeta1),n1,l1,m1,signum(zeta2),n2,l2,m2); Make_C(C2,signum(zeta3),n3,l3,m3,signum(zeta4),n4,l4,m4); Build_II(II1,mumax,sigma,N1,beta1); Build_II(II2,mumax,sigma,N2,beta2); Build_WY(WW1,mumax,sigma,N1,N2,alpha1,alpha2,tmax,jmax); (4./R)*K1*K2* sum('Mu_Term1(mu,sigma,N1,N2,C1,C2,II1,II2,WW1,alpha1,alpha2,beta1,beta2)', mu=sigma..mumax); end: printf("Two_el_2 loaded--small alpha version\n"); #----------------------------------------------------------------------------- Build_WY := proc(WW1,mumax,sigma,N1,N2,alpha1,alpha2,tmax,jmax) local Ibar1,Ibar2,LY1,LY2,LY12,mu,p1,p2,s,pps,S,k,j; global AC1; WW1 := array(sigma..mumax,0..N1,0..N2); Build_AC(mumax,sigma); Build_Ibar(Ibar1,mumax,sigma,N1,alpha1); Build_Ibar(Ibar2,mumax,sigma,N2,alpha2); Build_LY(LY1,mumax,sigma,N1,alpha1,tmax); Build_LY(LY2,mumax,sigma,N2,alpha2,tmax); Build_LY(LY12,mumax,sigma,mumax+sigma+N1+N2+jmax+1,alpha1+alpha2,tmax); for mu from sigma to mumax do for p1 from 0 to N1 do for p2 from 0 to N2 do S := 0.; for k from 0 to (mu+sigma)/2 do s := mu + sigma - 2*k; pps:= p1 + p2 + s + 1; S := S + AC1[mu,sigma,s] *((p2+s)!*sum('LY12[mu,pps+j]*alpha2^j/(p2+s+j+1)!',j=0..jmax) +(p1+s)!*sum('LY12[mu,pps+j]*alpha1^j/(p1+s+j+1)!',j=0..jmax)) od; WW1[mu,p1,p2] := ((mu+sigma)!/(mu-sigma)!)^2 *(S-LY1[mu,p1]*Ibar2[mu,p2]-Ibar1[mu,p1]*LY2[mu,p2]) od od od end: printf("Build_WY loaded\n"); #----------------------------------------------------------------------------- # Make function L, Eq. (28), by method for small alpha, using Eqs. # (41), (44), and (45), in LY1[mu,p] for mu through mumax, p through pmax, # for input sigma. MC1 is precalculated M coefficient array. # S0 is the finite summation in Eq. (41). # S2 is the summation in Eq. (44), retained through t=tmax. Build_LY := proc(LY1,mumax,sigma,pmax,alpha,tmax) local SS0,SS2,II1,E1,mu,p; LY1 := array(sigma..mumax, 0..pmax); Build_II(II1,mumax,sigma,pmax,alpha); Build_S0(SS0,mumax,sigma,pmax,alpha); Build_S2(SS2,mumax,sigma,pmax,alpha,tmax); E1 := evalf(Ei(1,alpha)); for mu from sigma to mumax do for p from 0 to pmax do LY1[mu,p] := SS0[mu,p] + E1*II1[mu,p]*(-1)^mu + SS2[mu,p] od od end: printf("LY (L for small alpha) loaded\n"); #----------------------------------------------------------------------------- # Evaluates finite summation, Eq. (41), in L evaluation for small alpha Build_S0 := proc(SS0,mumax,sigma,pmax,alpha) local AA1,mu,p,k; SS0 := array(sigma..mumax,0..pmax); Build_AA(AA1,pmax-1,alpha); for mu from sigma to mumax do for p from 0 to pmax do if p-mu+sigma-1 < 0 then SS0[mu,p] := 0 else SS0[mu,p] := (-1)^sigma*sum('(2*k+mu-sigma)! *AA1[p-2*k-mu+sigma-1]/(DF(2*k+2*mu+1)*DF(2*k))', k=0..(p-mu+sigma-1)/2) fi od od end: printf("S0 (term of LY) loaded\n"); #----------------------------------------------------------------------------- # Infinite series, Eq. (44), retained through t=tmax, in L evaluation # for small alpha Build_S2 := proc(SS2,mumax,sigma,pmax,alpha,tmax) local mu,p,t; global MC1; Build_MC(mumax,sigma,pmax,tmax); SS2 := array(sigma..mumax,0..pmax); for mu from sigma to mumax do for p from 0 to pmax do SS2[mu,p] := exp(-alpha)*sum('MC1[mu,sigma,p,t]*alpha^t', t=0..tmax) od od end: printf("S2 (term of LY) loaded\n"); #============================================================================= ################################# # Alternate Method, small alpha # ################################# #----------------------------------------------------------------------------- # Alternate calculation for small alpha, using Eq. (55). Build_KIA generates # the summation that is the first term of that equation; Build_GH makes the # following two sums. The final (unsummed) term is added in this procedure. # Additional terms (not calculated) are needed for nonzero sigma. Thus # calculations for nonzero sigma (also nonzero p1,p2) are blocked. Build_WYY := proc(WW1,mumax,sigma,p1,p2,alpha1,alpha2,tmax,jmax,pmax) local mu,alpha12,Ibar1,Ibar2,LY1,LY2,KIA1,KIA2,GH1; if sigma <> 0 or p1<>0 or p2<>0 then RETURN(0) fi; # Temporary WW1 := array(0..mumax); alpha12 := alpha1 + alpha2; Build_Ibar(Ibar1,mumax,sigma,0,alpha1); Build_Ibar(Ibar2,mumax,sigma,0,alpha2); Build_LY(LY1,mumax,sigma,0,alpha1,tmax); Build_LY(LY2,mumax,sigma,0,alpha2,tmax); Build_KIA(KIA1,mumax,sigma,alpha1,alpha12,pmax); Build_KIA(KIA2,mumax,sigma,alpha2,alpha12,pmax); Build_GH(GH1,mumax,sigma,alpha1,alpha2,tmax,jmax); for mu from sigma to mumax do WW1[mu] := (KIA1[mu]+KIA2[mu]+GH1[mu] -(LY1[mu,0]*Ibar2[mu,0]+Ibar1[mu,0]*LY2[mu,0])) *((mu+sigma)!/(mu-sigma)!)^2 od end: printf("Build_WYY--Eq. (55) version loaded\n"); #----------------------------------------------------------------------------- # IY1 is indexed [mu,q], where q=-p-1 in Eq. (55) Build_KIA := proc(KIA1,mumax,sigma,alpha,alpha12,pmax) local Kbar1,IY1,AA1,p,mu; Build_Kbar(Kbar1,mumax,sigma,0,alpha); Build_IY(IY1,mumax,sigma,pmax+1,alpha); Build_AA(AA1,pmax,alpha12); for mu from 0 to mumax do KIA1[mu] := (-1)^sigma *Kbar1[mu,0]*sum('(-1)^(p+1)*IY1[mu,-p-1]*AA1[p]',p=0..pmax) od end: printf("Build_KIA loaded\n"); #----------------------------------------------------------------------------- # Makes G and H terms according to Eqs. (57) and (58). Build_GH := proc(GH1,mumax,sigma,alpha1,alpha2,tmax,jmax) local SS2,II1,T,E1,X1,s,mu,j,k,alpha12; global AC1; alpha12 := alpha1 + alpha2; Build_AC(mumax,sigma); Build_S2X(SS2,mumax,sigma,mumax+sigma+jmax+1,alpha12,tmax); Build_II(II1,mumax,sigma,mumax+sigma+jmax+1,alpha12); E1 := evalf(Ei(1,alpha12)); X1 := exp(-alpha12); for mu from sigma to mumax do GH1[mu] := 0; for k from 0 to (mu+sigma)/2 do T := 0; s := mu+sigma-2*k; for j from 0 to jmax do T := T + ((alpha1^j+alpha2^j)/(s+j+1)!) *((-1)^mu*E1*II1[mu,s+j+1] + X1*SS2[mu,s+j+1]) od; GH1[mu] := GH1[mu] + s!*AC1[mu,sigma,s]*T od od end: printf("Build_GH loaded\n"); #----------------------------------------------------------------------------- # Infinite series, Eq. (44), retained through t=tmax, in L evaluation # for small alpha (not including exponential)--used for Eqs. (57)-(58). Build_S2X := proc(SS2,mumax,sigma,pmax,alpha,tmax) local mu,p,t; global MC1; Build_MC(mumax,sigma,pmax,tmax); SS2 := array(sigma..mumax,0..pmax); for mu from sigma to mumax do for p from 0 to pmax do SS2[mu,p] := sum('MC1[mu,sigma,p,t]*alpha^t', t=0..tmax) od od end: printf("S2X (term of GH) loaded\n"); #============================================================================= ########################## # Recursive procedures # ########################## #----------------------------------------------------------------------------- # L by upward recursion using Eq. (69) with exact starting values L_recur_up := proc(LL1,mumax,alpha) local A0,LL0,mu; Digits := Digits+20; Build_LX(LL0,2,0,0,alpha); Digits := Digits-20; A0 := exp(-alpha)/alpha; LL1[0] := LL0[0,0]; LL1[1] := LL0[1,0]; for mu from 2 to mumax do LL1[mu] := LL1[mu-2]+((2*mu-1)/alpha)*LL1[mu-1]-(2*mu-1)*A0/(mu*(mu-1)) od end: printf("L_recur_up loaded\n"); #----------------------------------------------------------------------------- # L by downward recursion using Eq. (69) with exact starting values # "opt" is a switch; if it is non-numeric or <> 0, adjust L # values according to Eqs. (74) and (75). L_recur_dn := proc(LL1,mumax,alpha,opt) local A0,LL0,II1,F,mu; Digits := Digits+60; Build_LX(LL0,mumax,0,0,alpha); Digits := Digits-60; A0 := exp(-alpha)/alpha; LL1[mumax] := LL0[mumax,0]; LL1[mumax-1] := LL0[mumax-1,0]; for mu from mumax-2 to 0 by -1 do LL1[mu] := LL1[mu+2]-((2*mu+3)/alpha)*LL1[mu+1] +(2*mu+3)*A0/((mu+1)*(mu+2)) od; if (not type(opt,numeric)) or opt <> 0 then Build_II(II1,mumax,0,0,alpha); F := (LL0[0,0]-LL1[0])/II1[0,0]; for mu from 0 to mumax do LL1[mu] := LL1[mu] + (-1)^mu*F*II1[mu,0] od fi end: printf("L_recur_dn loaded\n"); #----------------------------------------------------------------------------- # Recursively build L[mu,p] from L[mu,0], "vertical" process using # Eqs. (65)-(68) for mu=0 and Eqs. (70)-(73) for nonzero mu. # Recursion started with exact values of L[mu,0]. Lp_recur_vert := proc(LL1,mumax,pmax,alpha) local LL0,AA1,g,mu,p; Digits := Digits+80; Build_LX(LL0,mumax,0,0,alpha); Digits := Digits-80; for mu from 0 to mumax do LL1[mu,0] := LL0[mu,0] od; Build_AA(AA1,mumax-2,alpha); g[0] := alpha*LL1[0,0]; g[1] := g[0] - exp(alpha)*Ei(1,2*alpha); LL1[0,1] := (LL1[0,0]+g[1])/alpha; for p from 2 to pmax do g[p] := g[p-2] - AA1[p-2]; LL1[0,p] := (p*LL1[0,p-1] + g[p])/alpha od; for mu from 1 to mumax do g[0] := alpha*LL1[mu,0]; g[1] := alpha*LL1[mu-1,0] + mu*LL1[mu,0] - exp(-alpha)/mu; LL1[mu,1] := (LL1[mu,0]+g[1])/alpha; for p from 2 to pmax do g[p] := g[p-2] + mu*(LL1[mu,p-1]-LL1[mu-1,p-2]); LL1[mu,p] := (p*LL1[mu,p-1]+g[p])/alpha od od end: printf("Lp_recur_vert loaded\n"); #----------------------------------------------------------------------------- # Recursively build L[mu,p] from L[mu,0], "oblique" process using # Eqs. (65-68) for mu=0, Eqs. (70-73) for mu=mumax, and Eq. (60) for # other mu values. Recursion started with exact values of L[mu,0]. Lp_recur_obl := proc(LL1,mumax,pmax,alpha) local LL0,g,mu,p,AA1; Digits := Digits+80; Build_LX(LL0,mumax,0,0,alpha); Digits := Digits-80; for mu from 0 to mumax do LL1[mu,0] := LL0[mu,0] od; Build_AA(AA1,mumax-2,alpha); g[0] := alpha*LL1[0,0]; g[1] := g[0] - exp(alpha)*Ei(1,2*alpha); LL1[0,1] := (LL1[0,0]+g[1])/alpha; for p from 2 to pmax do g[p] := g[p-2] - AA1[p-2]; LL1[0,p] := (p*LL1[0,p-1] + g[p])/alpha od; g[0] := alpha*LL1[mumax,0]; g[1] := alpha*LL1[mumax-1,0] + mumax*LL1[mumax,0]-exp(-alpha)/mumax; LL1[mumax,1] := (LL1[mumax,0]+g[1])/alpha; for mu from 1 to mumax-1 do LL1[mu,1] := ((mu+1)*LL1[mu+1,0] + mu*LL1[mu-1,0])/(2*mu+1) od; for p from 2 to pmax do g[p] := g[p-2] + mumax*(LL1[mumax,p-1]-LL1[mumax-1,p-2]); LL1[mumax,p] := (p*LL1[mumax,p-1]+g[p])/alpha; for mu from 1 to mumax-1 do LL1[mu,p] := ((mu+1)*LL1[mu+1,p-1] + mu*LL1[mu-1,p-1])/(2*mu+1) od od end: printf("Lp_recur_obl loaded\n"); #----------------------------------------------------------------------------- # Recursively build L of nonzero sigma from a set of L[mu,p] of zero sigma, # using Eq. (61). Ls_recur := proc(LL1,mumax,sigma,p,alpha) local LL0,mu,s; Digits := Digits+80; Build_LX(LL0,mumax+sigma,0,p,alpha); Digits := Digits-80; for mu from 0 to mumax+sigma do LL1[mu,0] := LL0[mu,p] od; for s from 1 to sigma do for mu from s to mumax+sigma-s do LL1[mu,s] := (LL1[mu+1,s-1] - LL1[mu-1,s-1])/(2*mu+1) od od end: printf("Ls_recur loaded\n"); #----------------------------------------------------------------------------- # Make W[mu,p1,p2] array using Kotani method--First makes W[0,p1,p2] # for p1 and p2 through nmax, then increases mu, making W for p1,p2 # values for nmax-mu, so final output includes W[mu,0,0] through nmax. # Array ZZ is auxiliary function Z, Eq. (80). # Array SS1 is auxiliary function S, Eq. (83), p1=p2=0, alpha1 and alpha2 # as shown there. SS2 is with alpha1 and alpha2 reversed. Build_W_Kot := proc(WK1,mumax,alpha1,alpha2) local p1,p2,k,SS1,SS2,ZZ; Build_SS(SS1,mumax,alpha1,alpha2); Build_SS(SS2,mumax,alpha2,alpha1); ZZ := array(0..mumax,0..mumax,0..mumax); WK1 := array(0..mumax,0..mumax,0..mumax); for p1 from 0 to mumax do for p2 from 0 to mumax do WK1[0,p1,p2] := SS1[0,p1,p2] + SS2[0,p2,p1] od od; for p1 from 0 to mumax-1 do for p2 from 0 to mumax-1 do WK1[1,p1,p2] := SS1[1,p1,p2+1] + SS2[1,p2,p1+1] od od; for p1 from 0 to mumax-2 do for p2 from 0 to mumax-2 do ZZ[0,p1,p2] := SS1[0,p1,p2+2] + SS2[0,p2,p1+2] + SS1[1,p1+1,p2] + SS2[1,p2+1,p1] - 3*WK1[1,p1+1,p2+1]; WK1[2,p1,p2] := (WK1[0,p1,p2] - 3*ZZ[0,p1,p2])/4 od od; for p1 from 0 to mumax-3 do for p2 from 0 to mumax-3 do ZZ[1,p1,p2] := 3*WK1[1,p1,p2+2] + 3*WK1[1,p1+2,p2] - WK1[0,p1+1,p2+1] - WK1[1,p1,p2] - 5*WK1[2,p1+1,p2+1]; WK1[3,p1,p2] := (4*WK1[1,p1,p2] - 5*ZZ[1,p1,p2])/9 od od; for k from 2 to mumax-2 do for p1 from 0 to mumax-k-2 do for p2 from 0 to mumax-k-2 do ZZ[k,p1,p2] := ZZ[k-2,p1,p2] + (2*k+1)*(WK1[k,p1+2,p2] + WK1[k,p1,p2+2]) - (2*k-1)*WK1[k-1,p1+1,p2+1] - (2*k+3)*WK1[k+1,p1+1,p2+1]; WK1[k+2,p1,p2] := ((k+1)^2*WK1[k,p1,p2] - (2*k+3)*ZZ[k,p1,p2])/(k+2)^2 od od od; end: printf("Build_W_Kot loaded\n"); #----------------------------------------------------------------------------- # Auxiliary function S, Eq. (83), for Kotani's recurrence method for W. # Evaluation according to Eq. (87). Precision increased temporarily so # starting values will be certainly accurate. Build_SS := proc(SS1,mumax,alpha1,alpha2) local p1,p2,mu,LX1,LX12; Digits := Digits+40; SS1 := array(0..1,0..mumax,0..mumax); Build_LX(LX1,1,0,mumax,alpha1); Build_LX(LX12,1,0,2*mumax,alpha1+alpha2); Digits := Digits-40; for mu from 0 to 1 do for p1 from 0 to mumax do SS1[mu,p1,0] := (exp(-alpha2)*LX1[mu,p1] - LX12[mu,p1])/alpha2; for p2 from 1 to mumax do SS1[mu,p1,p2] := (p2*SS1[mu,p1,p2-1] + exp(-alpha2)*LX1[mu,p1] - LX12[mu,p1+p2])/alpha2 od od od end: printf("Build_SS loaded\n"); #----------------------------------------------------------------------------- # Makes W[[nu,mu] for arguments alpha1, alpha2, from inhomogeneous formulas # to maximum mu and mu = mumax for p=sigma=0. # Result in array WX1. First index of WX1 is for Q, second is for P. # Upward recurrence, Eqs. (96) and (97), with exact starting values of W. Recur_W_up := proc(WX1,mumax,alpha1,alpha2) local i,alpha12,WW1,TX1; alpha12 := alpha1+alpha2; WX1 := array(0..mumax,0..mumax); Digits := Digits+20; Build_WX(WW1,1,0,1,1,alpha1,alpha2); Digits := Digits-20; WX1[0,0] := WW1[0,0,0]; WX1[1,0] := WW1[0,1,0] - evalf(exp(-alpha12)/(alpha1*alpha12)); WX1[0,1] := WW1[0,0,1] - evalf(exp(-alpha12)/(alpha2*alpha12)); WX1[1,1] := WW1[1,0,0]; Build_TX(TX1,2*mumax,alpha12); for i from 0 to mumax-2 do WX1[i+2,i] := WX1[i,i]+((2*i+3)/alpha1)*WX1[i+1,i]+(TX1[i+2,i] -TX1[i,i+2])/alpha1; WX1[i+2,i+1]:=WX1[i,i+1]+((2*i+3)/alpha1)*WX1[i+1,i+1]+(TX1[i+1,i] -TX1[i+1,i+2]+TX1[i+2,i+1]-TX1[i,i+1])/alpha1; WX1[i+1,i+2]:=WX1[i+1,i]+((2*i+3)/alpha2)*WX1[i+1,i+1]+(TX1[i+1,i] -TX1[i,i+1]-TX1[i+1,i+2]+TX1[i+2,i+1])/alpha2; WX1[i+2,i+2] := WX1[i+2,i]+((2*i+3)/alpha2)*WX1[i+2,i+1] +(TX1[i+2,i]-TX1[i,i+2])/alpha2; od; end: printf("Recur_W_up loaded\n"); #----------------------------------------------------------------------------- # Makes W[[nu,mu] for arguments alpha1, alpha2, from inhomogeneous formulas # to maximum mu and mu = mumax for p=sigma=0. # Result in array WX1. First index of WX1 is for Q, second is for P. # Downward recurrence, Eqs. (96) and (97), with exact starting values of W. # If opt numeric and nonzero, adjust values of W according to Eqs. (109) # and (110). Only fills out off-diagonal part of W array if adjusted # values are requested. Recur_W_dn := proc(WY1,mumax,alpha1,alpha2,opt) local alpha12,TX1,IX1,IX2,W0,F,T,mm,mu,nu,mm1; alpha12 := alpha1 + alpha2; Digits := Digits+60; Build_WXXX(WY1,W0,mumax,alpha1,alpha2); Digits := Digits-60; Build_TX(TX1,2*mumax,alpha12); Recur_Dn(WY1,TX1,mumax,alpha1,alpha2); if type(opt,numeric) and opt=0 then RETURN(0) fi; Build_II(IX1,mumax,0,0,alpha1); Build_II(IX2,mumax,0,0,alpha2); mm := mumax; while mm > 0 do F := (W0-WY1[0,0])/(IX1[0,0]*IX2[0,0]); mm1 := mm; for mu from mm to 0 by -1 do T := WY1[mu,mu] + F*IX1[mu,0]*IX2[mu,0]; if (abs(WY1[mu,mu]/T) < 10^5) then for nu from 0 to mu-1 do WY1[mu,nu] := WY1[mu,nu] + F*(-1)^(mu+nu)*IX1[mu,0]*IX2[nu,0]; WY1[nu,mu] := WY1[nu,mu] + F*(-1)^(mu+nu)*IX1[nu,0]*IX2[mu,0] od; WY1[mu,mu]:=T; mm := mu; else Recur_Dn(WY1,TX1,mu+2,alpha1,alpha2); break; fi; od od end: printf("Recur_W_dn loaded\n"); #----------------------------------------------------------------------------- Recur_Dn := proc(WY1,TX1,mumax,alpha1,alpha2) local i,j,j1,j2,T; for i from mumax to 0 by -1 do j1 := min(mumax-2,i); j2 := min(mumax-2,i-1); for j from j1 to 0 by -1 do T := TX1[i,j+2]-TX1[i,j]-TX1[j+2,i]+TX1[j,i]; WY1[j,i] := WY1[j+2,i] - ((2*j+3)*WY1[j+1,i]-T)/alpha1 od; for j from j2 to 0 by -1 do T := TX1[i,j+2]-TX1[i,j]-TX1[j+2,i]+TX1[j,i]; WY1[i,j] := WY1[i,j+2] - ((2*j+3)*WY1[i,j+1]-T)/alpha2 od od end: printf("Recur_Dn loaded\n"); #----------------------------------------------------------------------------- # Explicit computation of W(nu,mu,alpha1,alpha2) by a generalization of # Eq. (39), for nu,mu = mumax or mumax-1. Build_WXXX := proc(WXX,W0,mumax,alpha1,alpha2) local KX1,KX2,LX1,LX2,LX12,alpha12,s,k,j,S1,S2,S3,S4,S5,S6,S7,S8; global AC1; alpha12 := alpha1 + alpha2; Build_AC(mumax,0); Build_KX(KX1,mumax,0,0,alpha1); Build_KX(KX2,mumax,0,0,alpha2); Build_LX(LX1,mumax,0,0,alpha1); Build_LX(LX2,mumax,0,0,alpha2); Build_LX(LX12,mumax,0,mumax,alpha12); S1 := 0; S2 := 0; S3 := 0; S4 := 0; for k from 0 to mumax/2 do s := mumax-2*k; S1:=S1 -AC1[mumax,0,s]*s!*sum('LX12[mumax,j]/(j!*alpha1^(s-j+1))',j=0..s); S2:=S2 -AC1[mumax,0,s]*s!*sum('LX12[mumax,j]/(j!*alpha2^(s-j+1))',j=0..s); S3:=S3 -AC1[mumax,0,s]*s!*sum('LX12[mumax-1,j]/(j!*alpha1^(s-j+1))',j=0..s); S4:=S4 -AC1[mumax,0,s]*s!*sum('LX12[mumax-1,j]/(j!*alpha2^(s-j+1))',j=0..s); od; S5 := 0; S6 := 0; S7 := 0; S8 := 0; for k from 0 to (mumax-1)/2 do s := mumax-1-2*k; S5:=S5 -AC1[mumax-1,0,s]*s!*sum('LX12[mumax,j]/(j!*alpha1^(s-j+1))',j=0..s); S6:=S6 -AC1[mumax-1,0,s]*s!*sum('LX12[mumax,j]/(j!*alpha2^(s-j+1))',j=0..s); S7:=S7 -AC1[mumax-1,0,s]*s!*sum('LX12[mumax-1,j]/(j!*alpha1^(s-j+1))',j=0..s); S8:=S8 -AC1[mumax-1,0,s]*s!*sum('LX12[mumax-1,j]/(j!*alpha2^(s-j+1))',j=0..s); od; WXX[mumax,mumax] := S1 + S2 + LX1[mumax,0]*KX2[mumax,0] + KX1[mumax,0]*LX2[mumax,0]; WXX[mumax-1,mumax-1] := S7 + S8 + LX1[mumax-1,0]*KX2[mumax-1,0] + KX1[mumax-1,0]*LX2[mumax-1,0]; WXX[mumax,mumax-1] := S3 + S6 + LX1[mumax,0]*KX2[mumax-1,0] + KX1[mumax,0]*LX2[mumax-1,0]; WXX[mumax-1,mumax] := S4 + S5 + LX1[mumax-1,0]*KX2[mumax,0] + KX1[mumax-1,0]*LX2[mumax,0]; W0 := LX1[0,0]*KX2[0,0]+KX1[0,0]*LX2[0,0]-LX12[0,0]*(1/alpha1+1/alpha2); end: printf("Build_WXXX loaded\n"); #----------------------------------------------------------------------------- # Make W for nonzero p1 and p2, using Eqs. (103) and (108), through p1max # and p2max. Recur_Wp := proc(WYP,mumax,p1max,p2max,alpha1,alpha2) local WY1,WY2,TT,p1,p2,mu,nu,pmax; pmax := max(p1max,p2max); Recur_W_dn(WY1,mumax+pmax,alpha1,alpha2,1); Build_Sbar(SS1,mumax+p2max,p1max,alpha1,alpha2); Build_Sbar(SS2,mumax+p1max,p2max,alpha2,alpha1); for p2 from 0 to pmax do # Compute array for p1=0 and current p2, place in WY1 if p2 <> 0 then for nu from 0 to mumax+p1max do TT[0] := WY1[nu,1] + SS2[nu,p2-1,0]; for mu from 1 to mumax+p2max-p2 do TT[mu] := ((mu+1)*WY1[nu,mu+1]+mu*WY1[nu,mu-1])/(2*mu+1) od; for mu from 0 to mumax+p2max-p2 do WY1[nu,mu] := TT[mu] od od fi; # Save W for p1=0 and current p2 in WYP for mu from 0 to mumax do WYP[mu,0,p2] := WY1[mu,mu] od; # Copy WY1 to WY2 for nu from 0 to mumax+p1max do for mu from 0 to mumax+p2max-p2 do WY2[nu,mu] := WY1[nu,mu] od od; # Advance p1 for the current p2 for p1 from 1 to p1max do for mu from 0 to mumax+p2max-p2 do TT[0] := WY2[1,mu] + SS1[mu,p1-1,p2]; for nu from 1 to mumax+p1max-p1 do TT[nu] := ((nu+1)*WY2[nu+1,mu]+nu*WY2[nu-1,mu])/(2*nu+1) od; for nu from 0 to mumax+p1max-p1 do WY2[nu,mu] := TT[nu] od od; # Save W for current p1 and p2 in WYP for mu from 0 to mumax do WYP[mu,p1,p2] := WY2[mu,mu] od od od end: printf("Recur_Wp loaded\n"); #----------------------------------------------------------------------------- # Sbar, Eq. (105), using Eq. (106). Build_Sbar := proc(SS1,mumax,pmax,alpha1,alpha2) local mu,p1,p2,KX1; Build_KX(KX1,mumax,0,2*pmax,alpha1+alpha2); for mu from 0 to mumax do for p2 from 0 to pmax do SS1[mu,0,p2] := KX1[mu,p2]/alpha1; for p1 from 1 to pmax do SS1[mu,p1,p2] := (p1*SS1[mu,p1-1,p2] + KX1[mu,p1+p2])/alpha1 od od od end: printf("Build_Sbar loaded\n"); #----------------------------------------------------------------------------- # Makes T_nu,mu(alpha) = int(Q_nu(t) P_mu(t) exp(-alpha xi),xi=1..infinity) # defined in Eq. (98) in array TX1[nu,mu] for indices to mumax # High precision for starting values, current precision for iteration. # Procedure in Eqs. (157)-(160). Build_TX := proc(TX1,mumax,alpha) local LX1,KX1,nu,mu; Digits := Digits+60; Build_LX(LX1,mumax,0,0,alpha); Build_KX(KX1,mumax,0,0,alpha); Digits := Digits-60; TX1 := array(-1..mumax,-1..mumax,[]); for nu from 0 to mumax do TX1[nu,-1] := 0.; TX1[-1,nu] := KX1[nu,0]; # For mu = nu TX1[nu,0] := LX1[nu,0] od: for mu from 1 to mumax do TX1[0,mu]:=((TX1[-1,mu-1]+TX1[1,mu-1])*(2*mu-1)-(mu-1)*TX1[0,mu-2])/mu; for nu from 1 to mumax-mu do TX1[nu,mu]:=((nu*TX1[nu-1,mu-1]+(nu+1)*TX1[nu+1,mu-1])*(2*mu-1)/(2*nu+1) - (mu-1)*TX1[nu,mu-2])/mu od od end: printf("Build_TX loaded\n"); #----------------------------------------------------------------------------- # From an array of W[mu,0,p1,p2] made with procedure Recur_Wp, uses Eq. (77) # to advance sigma to the input value.Puts W[mu,sigma,p1,p2] in WY1[mu,p1,p2] # for the input sigma value, mu through mumax, p1,p2 through N1,N2. Recur_Ws := proc(WY1,mumax,sigma,N1,N2,alpha1,alpha2) local p1,p2,mu,s; Recur_Wp(WY1,mumax+sigma,N1+sigma,N2+sigma,alpha1,alpha2); # Initially stores results shifted down in mu by 1 each step in sigma for s from 1 to sigma do for p1 from 0 to N1+sigma-s do for p2 from 0 to N2+sigma-s do for mu from s to mumax+sigma-s do WY1[mu-s,p1,p2] := ((mu-s+1)*(mu-s+2)^2*WY1[mu-s+2,p1,p2] + (mu+s)*(mu+s-1)^2*WY1[mu-s,p1,p2])/(2*mu+1) - (mu-s+1)*(mu+s)*WY1[mu-s+1,p1+1,p2+1] od od od od; for p1 from 0 to N1 do for p2 from 0 to N2 do for mu from mumax by -1 to sigma do WY1[mu,p1,p2] := WY1[mu-sigma,p1,p2] od od od end: printf("Recur_Ws loaded\n"); #============================================================================= ################################################### # Utilities, Special Functions,Coefficient sets # ################################################### #----------------------------------------------------------------------------- # Generalized spherical Bessel function i, Eq. (11), for general q >= 0. # Not designed to be efficient; recursive construction # starting from function with q=0, from procedure II0. II := proc(mu,sigma,q,beta) local S; if q=0 then II0(mu,sigma,beta) else if L=M then S := 0 else S := (mu-sigma)*II(mu-1,sigma,q-1,beta) fi; -((mu+sigma+1)*II(mu+1,sigma,q-1,beta)+S)/(2*mu+1) fi end: printf("II, Generalized spherical Bessel function i, loaded\n"); #----------------------------------------------------------------------------- # Generalized spherical Bessel function i, Eq. (11), for q=0. # Not designed to be efficient; uses Maple procedure for modified # Bessel function BesselI. II0 := proc(mu,sigma,beta); if sigma=0 then if beta=0 then if mu=0 then 1 else 0 fi elif beta<0 then (-1)^mu*II0(mu,0,-beta) else evalf((Pi/(2*beta))^(1/2)*BesselI(mu+1/2,beta)) fi else evalf((1./(2*mu+1))*(II0(mu+1,sigma-1,beta)-II0(mu-1,sigma-1,beta))) fi end: printf("II0, Generalized spherical Bessel function i (q=0) loaded\n"); #----------------------------------------------------------------------------- # Make generalized Bessel function i for negative q, defined by Eq. (22). # Evaluation using Eqs. (17) and (161)--(163). # IY(mu,sigma,q,beta) = i(mu,sigma,-q-1,beta) is evaluated only for q > 0. IY := proc(mu,sigma,q,beta); if sigma>0 then (IY(mu-1,sigma-1,q,beta) - IY(mu+1,sigma-1,q,beta))/(2*mu+1) elif mu>1 and q<-1 then -((mu-1)*IY(mu-2,0,q,beta)+(2*mu-1)*IY(mu-1,0,q+1,beta))/mu elif mu=1 and q<-1 then (-1)^(q+1)*beta^(-q-1)/(-q-1)! - IY(0,0,q+1,beta) elif mu>1 then -((mu-1)*IY(mu-2,0,-1,beta)+(2*mu-1)*II(mu-1,0,0,beta))/mu elif mu=1 then 1-II0(0,0,beta) elif q<-3 then (beta*(IY(0,0,q+1,beta)-IY(0,0,q+3,beta)) +(q+3)*IY(0,0,q+2,beta)+(-1)^(q+1)*beta^(-q-2)/(-q-2)!)/(q+1) elif q=-3 then (beta*(IY(0,0,-2,beta)-II0(0,0,beta)) +beta)/(-2) elif q=-2 then (beta*(IY(0,0,-1,beta)-II(0,0,1,beta))+II0(0,0,beta) -1)/(-1) else -Shi(beta) fi end: printf("IY, Generalized spherical Bessel function i, (q<0), loaded\n"); #----------------------------------------------------------------------------- # Make generalized spherical Bessel function k, Eq. (29), # using Eqs. (31) and (32): KX := proc(mu,sigma,p,alpha); if sigma=0 then if p=0 then KX0(mu,alpha) elif mu=0 then KX(1,0,p-1,alpha) else ((mu+1)*KX(mu+1,0,p-1,alpha) + mu*KX(mu-1,0,p-1,alpha))/(2*mu+1) fi else (KX(mu+1,sigma-1,p,alpha) - KX(mu-1,sigma-1,p,alpha))/(2*mu+1) fi end: printf("KX, Generalized spherical Bessel function k, loaded\n"); #----------------------------------------------------------------------------- # Make modified spherical Bessel function k_n(alpha), using Eq. (30), # at non-standard scaling such that k0(alpha) = exp(-alpha)/alpha KX0 := proc(n,alpha); if n=-1 or n=0 then exp(-alpha)/alpha else KX0(n-2,alpha) + (2*n-1)*KX0(n-1,alpha)/alpha fi end: printf("KX0, modified spherical Bessel function k, loaded\n"); #----------------------------------------------------------------------------- # Complementary i function of Eq. (53), for sigma = 0. # Evaluated using Eq. (49) with xi = 1. I_bar := proc(mu,sigma,p,alpha) local k; sum('AC(mu,sigma,mu+sigma-2*k)*Aa(p+mu+sigma-2*k,alpha)', k=0..(mu+sigma)/2) end: printf("I_bar, complementary i function, loaded\n"); #----------------------------------------------------------------------------- # AA is the regular A sub n integral for argument x, Eq. (34), # calculated using Eq. (35). Method works for both positive and # negative x. AA := proc(n,x) local j,A0,S; A0 := exp(-x)/x; S := A0; for j from 0 to n do S := (j/x)*S + A0 od; evalf(S) end: printf("AA, A sub n integral, loaded\n"); #----------------------------------------------------------------------------- # Complementary A integral of Eq. (50), evaluated according to Eq. (52): # To be used for small enough alpha that expansion converges after alpha^16. Aa := proc(n,alpha) local j; if alpha<>0 then n!*exp(-alpha)*sum('alpha^j/(n+j+1)!',j=0..16) else 1./(n+1) fi end: printf("Aa, complementary A sub n, loaded\n"); #----------------------------------------------------------------------------- # Sum (k=0..infinity) 1/(2k+i)(2k+j), i,j>0, does not include case i=j. # Formulas in Eq. (46). T_sum := proc(i,j) local k,i1,j1,S; options remember; if i mod 2 = j mod 2 then i1:=min(i,j); j1:=max(i,j); S := (1/(j1-i1))*sum('1/(2*k+i1)',k=0..(j1-i1-2)/2) else if i mod 2 = 0 then i1:=i; j1:=j else i1:=j; j1:=i fi; if j1 > 1 then S := sum('1/(2*k+1)',k=0..(j1-3)/2) else S := 0 fi; S := (1/(i1-j1))*(log(2)+sum('1/(2*k)',k=1..(i1-2)/2)-S) fi; end: printf("T_sum loaded\n"); #----------------------------------------------------------------------------- # Coefficient A, Eq. (26). # DF(n) is n!! AC := proc(mu,sigma,s); (-1)^((s-mu-sigma)/2)*DF(mu-sigma+s-1)/(s!*DF(mu+sigma-s)) end: printf("A coefficient loaded\n"); #----------------------------------------------------------------------------- # Coefficients B, Eq. (27). BC := proc(mu,sigma,s) local k; -sum('(-1)^k*DF(2*mu-2*k-1)/ ((mu+sigma-s-2*k)*(mu+sigma-2*k)!*DF(2*k))',k=0..(mu+sigma-s-1)/2) end: printf("B coefficient loaded\n"); #----------------------------------------------------------------------------- # Coefficient LC(p,t), Eq. (37), used in procedure L_Sum: LC := proc(p,t) local S,k,m,m1,m2; S := -sum('1/k',k=1..p-t); for k from 0 to p-1 do m1 := max(0,t-p+k+1); m2 := min(t,k); S := S + sum('(-1)^m*2^(t-m)*t!/(m!*(t-m)!)',m=m1..m2)/(p-k) od; S end: printf("LC, L coefficient loaded\n"); #----------------------------------------------------------------------------- # Coefficient M(mu,sigma,p,t), Eq. (45). # T_sum is a summation evaluated in Eq. (46). # twok1 is 2 k_1. MC := proc(mu,sigma,p,t) local S,j,l,twok1; S := 0; twok1 := max(t+p-mu+sigma+1,0); twok1 := twok1 + (twok1 mod 2); for j from 0 to (mu+sigma)/2 do S := S + sum('(-1)^(l+j)*DF(2*mu-2*j-1) *T_sum(twok1+2*mu-2*j+1,twok1+mu-sigma-p-l) /((mu+sigma-2*j)!*DF(2*j)*l!*(t-l)!)',l=0..t) od end: printf("MC, M coefficient loaded\n"); #----------------------------------------------------------------------------- # Make n!! (including (-1)!!=1): DF := proc(n) options remember; if n=-1 then RETURN(1) fi; if n mod 2 = 0 then 2^(n/2)*(n/2)! else n!/(2^((n-1)/2)*((n-1)/2)!) fi; end: printf("DF, makes n!!, loaded\n"); #----------------------------------------------------------------------------- # Make (t^2-1)^sigma Q mu sigma(t), where Q is an associated Legendre # function of the second kind. Slow, not simplified. QQ := proc(mu,sigma,t) local T,x,j; T := (1/2)*log((x+1)/(x-1)); if mu<>0 then with(orthopoly); T := P(mu,x)*T - sum('P(j-1,x)*P(mu-j,x)/j',j=1..mu); if sigma <> 0 then T:=(x^2-1)^sigma*diff(T,x$sigma) else T fi; fi; eval(subs(x=t,T)) end: printf("QQ loaded\n"); #============================================================================= ############################### # Array-building procedures # ############################### # If maximum index values are determined in advance, it is more efficient # to make arrays of quantities that will be used more than once, and to # build sets of functions recursively in a coordinated fashion. The # procedures in this section are for this purpose, and include coefficients # and functions for which a single evaluation method is sufficient. The # procedures in later sections are written to use these arrays. #----------------------------------------------------------------------------- # Make coefficients AC, put into array AC1[mu,sigma,s] for mu values # through maxmu, sigma values through maxsigma, all s. Build_AC := proc(maxmu,maxsigma) local mu,sigma,s; global AC1,AC_mu,AC_sigma; if not type(AC_mu,integer) then AC_mu:= -1 fi; if not type(AC_sigma,integer) then AC_sigma:= -1 fi; for mu from 0 to AC_mu do for sigma from AC_sigma+1 to min(mu,maxsigma) do for s from mu+sigma by -2 to 0 do AC1[mu,sigma,s] := AC(mu,sigma,s) od od od; for mu from AC_mu+1 to maxmu do for sigma from 0 to min(mu,maxsigma) do for s from mu+sigma by -2 to 0 do AC1[mu,sigma,s] := AC(mu,sigma,s) od od od; AC_mu := maxmu; AC_sigma := maxsigma; end: printf("Build_AC loaded\n"); #----------------------------------------------------------------------------- # Make coefficients BC, put into array BC1[mu,sigma,s] for mu values # through maxmu, sigma values through maxsigma, all s. Build_BC := proc(maxmu,maxsigma) local mu,sigma,s; global BC1,BC_mu,BC_sigma; if not type(BC_mu,integer) then BC_mu:= -1 fi; if not type(BC_sigma,integer) then BC_sigma:= -1 fi; for mu from 0 to BC_mu do for sigma from BC_sigma+1 to min(mu,maxsigma) do for s from mu+sigma-1 by -2 to 0 do BC1[mu,sigma,s] := BC(mu,sigma,s) od od od; for mu from BC_mu+1 to maxmu do for sigma from 0 to min(mu,maxsigma) do for s from mu+sigma-1 by -2 to 0 do BC1[mu,sigma,s] := BC(mu,sigma,s) od od od; BC_mu := maxmu; BC_sigma := maxsigma; end: printf("Build_BC loaded\n"); #----------------------------------------------------------------------------- # Make coefficients LC, put into array LC1[p,t] for p values # through maxp, all t. Based on Eq. (38). Build_LC := proc(maxp) local p,S,t,j,k; global LC1,LC_p; if not type(LC_p,integer) then LC_p := -1 fi; if maxp <= LC_p then RETURN(0) fi; for t from 1 to maxp-1 do S[t] := 2^t; for j from t-1 by -1 to 1 do S[j] := S[j+1] + (-1)^(t-j)*2^j*t!/(j!*(t-j)!) od; for p from max(t+1,LC_p+1) to maxp do LC1[p,t] := sum('S[k]*(1/(k+p-t) - 1/k)',k=1..t) od od; LC_p := maxp end: printf("Build_LC loaded\n"); #----------------------------------------------------------------------------- # Make coefficients M, Eq. (45), for mu through maxmu, sigma through # maxsigma, p through maxp, t through maxt. Build_MC := proc(maxmu,maxsigma,maxp,maxt) local mu,sigma,p,t,p0,t0,t1; global MC_mu,MC_sigma,MC_p,MC_t,MC1; if not type(MC_mu,integer) then MC_mu := -1 fi; if not type(MC_sigma,integer) then MC_sigma := -1 fi; if not type(MC_p,integer) then MC_p := -1 fi; if not type (MC_t,integer) then MC_t := -1 fi; if maxmu <= MC_mu and maxsigma <= MC_sigma and maxp <= MC_p and maxt <= MC_t then RETURN(0) fi; for mu from 0 to maxmu do for sigma from 0 to min(mu,maxsigma) do if mu > MC_mu or sigma > MC_sigma then p0:=0; t0:=0 else p0:=MC_p+1; t0:=MC_t+1 fi; for p from 0 to maxp do if p >= p0 then t1 := 0 else t1 := t0 fi; for t from t1 to maxt do MC1[mu,sigma,p,t] := MC(mu,sigma,p,t) od od od od; MC_mu:=maxmu; MC_sigma:=maxsigma; MC_p:=maxp; MC_t:=maxt end: printf("Build_MC loaded\n"); #----------------------------------------------------------------------------- # Build functions A sub n of alpha, put into array AA1[n] for n values # through maxn. Build_AA := proc(AA1,maxn,alpha) local n; if maxn < 0 then RETURN(0) fi; AA1 := array(0..maxn,[]); AA1[0] := exp(-alpha)/alpha; for n from 1 to maxn do AA1[n] := (n/alpha)*AA1[n-1]+AA1[0] od end: printf("Build_AA (A sub n) loaded\n"); #----------------------------------------------------------------------------- # Build functions Aa sub n of alpha, put into array Aa1[n] for n values # through maxn. This procedure is only used for small alpha. For large # alpha, it would be preferable to calculate as n!/alpha^(n+1)-A_n(alpha). Build_Aa := proc(Aa1,maxn,alpha) local NN,n,S; Aa1 := array(0..maxn,[]); Aa1[maxn] := 1.; NN := trunc(maxn + 10*(abs(alpha)+2)); for n from NN to maxn by -1 do Aa1[maxn] := n/(n+1+alpha*(1-Aa1[maxn])) od; for n from maxn-1 to 1 by -1 do Aa1[n] := n/(n+1+alpha*(1-Aa1[n+1])) od; if alpha <> 0 then Aa1[0] := (1-exp(-alpha))/alpha else Aa1[0] := 1 fi; for n from 1 to maxn do Aa1[n] := Aa1[n-1]*Aa1[n] od end: printf("Build_Aa loaded\n"); #----------------------------------------------------------------------------- # Make functions II for argument beta, put into array II1[mu,q] # for mu through maxmu, input sigma, q through maxq. Build_II := proc(II1,maxmu,sigma,maxq,beta) local mu,j,q,S; II1 := array(0..maxmu+maxq,0..maxq,[]); if abs(beta) < 10^(-10) then II1[sigma,0] := 1.; for j from 1 to sigma do II1[sigma,0] := II1[sigma,0]/(2*j+1) od; for mu from sigma+1 to maxmu+maxq do II1[mu,0] := 0 od else for mu from sigma to maxmu+maxq do II1[mu,0] := evalf(sqrt(Pi/(2*beta))*BesselI(mu+1/2,beta)/beta^sigma) od fi; for q from 1 to maxq do for mu from sigma to maxmu+maxq-q do if mu=sigma then S:=0 else S := (mu-sigma)*II1[mu-1,q-1] fi; II1[mu,q] := -(S+(mu+sigma+1)*II1[mu+1,q-1])/(2*mu+1) od od end: printf("Build_II (generalized i) loaded\n"); #----------------------------------------------------------------------------- # Make functions I, Eq. (22), for negative q for mu through maxmu, # sigma through maxsigma, q through negative values through -maxmq. # Algorithm not entirely stable. Temporary solution: Increase # precision at start of procedure, restore original precision at end. Build_IY := proc(II1,maxmu,sigma,maxmq,beta) local II2,II3,mu,ss,mq; if maxmq <= 0 then RETURN() fi; Digits := Digits+10; Build_II(II2,max(maxmu+sigma-1,1),0,0,beta); for mu from 0 to max(maxmu+sigma-1,1) do II3[mu,0,0] := II2[mu,0] od; II3[0,0,-1] := -Shi(beta); if maxmq > 1 then II3[0,0,-2] := 1-beta*(II3[1,0,0]+II3[0,0,-1])-II3[0,0,0] fi; for mq from 3 to maxmq do II3[0,0,-mq] := (beta*(II3[0,0,1-mq]-II3[0,0,3-mq]) +(3-mq)*II3[0,0,2-mq]-(-1)^mq*beta^(mq-2)/(mq-2)!)/(1-mq) od; if maxmu+sigma > 0 then for mq from 1 to maxmq do II3[1,0,-mq] := -II3[0,0,1-mq] - (-1)^mq*beta^(mq-1)/(mq-1)!; for mu from 2 to maxmu+sigma do II3[mu,0,-mq]:=-((mu-1)*II3[mu-2,0,-mq]+(2*mu-1)*II3[mu-1,0,1-mq])/mu od; for ss from 1 to sigma do for mu from ss to maxmu+sigma-ss do II3[mu,ss,-mq]:=(II3[mu-1,ss-1,-mq]-II3[mu+1,ss-1,-mq])/(2*mu+1) od od od fi; for mu from sigma to maxmu do for mq from 1 to maxmq do II1[mu,-mq] := II3[mu,sigma,-mq] od od; Digits := Digits-10; end: printf("Build_IY (generalized i, negative q) loaded\n"); #----------------------------------------------------------------------------- # Make functions KX for argument alpha, put into array KX1[mu,p] # for mu through maxmu, input sigma, p through maxp. Build_KX := proc(KX1,maxmu,sigma,maxp,alpha) local KX2,mu,k,p; KX1 := array(sigma..maxmu+maxp,0..maxp,[]); KX2[0,0] := exp(-alpha)/alpha; KX2[-1,0] := KX2[0,0]; for mu from 1 to maxmu+sigma+maxp do KX2[mu,0] := KX2[mu-2,0] + (2*mu-1)*KX2[mu-1,0]/alpha od; for k from 1 to sigma do for mu from k to maxmu+sigma+maxp-k do KX2[mu,k] := (KX2[mu+1,k-1]-KX2[mu-1,k-1])/(2*mu+1) od od; for mu from sigma to maxmu+maxp do KX1[mu,0] := KX2[mu,sigma] od; for p from 1 to maxp do KX1[sigma,p] := KX1[sigma+1,p-1]; for mu from sigma+1 to maxmu+maxp-p do KX1[mu,p] := ((mu+sigma+1)*KX1[mu+1,p-1] + (mu-sigma)*KX1[mu-1,p-1])/(2*mu+1) od od end: printf("Build_KX (generalized k) loaded\n"); #----------------------------------------------------------------------------- # Make functions ibar, Eq. (53), for argument alpha, put in array Ibar1[mu,p] # for mu through maxmu, input sigma, and p through maxp. # AC1 is the A coefficient array, Eq. (26). Build_Ibar := proc(Ibar1,maxmu,sigma,maxp,alpha) local Aa1,k,mu,p; global AC1; Ibar1 := array(sigma..maxmu+maxp,0..maxp,[]); Build_AC(maxmu+maxp,sigma); Build_Aa(Aa1,maxmu+sigma+maxp,alpha); for mu from sigma to maxmu+maxp do Ibar1[mu,0] := sum('AC1[mu,sigma,mu+sigma-2*k]*Aa1[mu+sigma-2*k]', k=0..(mu+sigma)/2) od; for p from 1 to maxp do Ibar1[sigma,p] := Ibar1[sigma+1,p-1]; for mu from sigma+1 to maxmu+maxp-p do Ibar1[mu,p] := ((mu+sigma+1)*Ibar1[mu+1,p-1]+(mu-sigma)*Ibar1[mu-1,p-1]) /(2*mu+1) od od end: printf("Build_Ibar loaded\n"); #----------------------------------------------------------------------------- # Make functions kbar, Eq. (56), for argument alpha, put in array # Kbar1[mu,p] for mu through maxmu, input sigma, and p through maxp. # AC1 is the A coefficient array, Eq. (26). Build_Kbar := proc(Kbar1,maxmu,sigma,maxp,alpha) local K1,mu,ss,p; Build_AC(maxmu+maxp+sigma,sigma); for mu from 0 to maxmu+maxp+sigma do K1[mu,0] := sum('(mu-2*k)!*AC1[mu,0,mu-2*k]*alpha^(2*k-mu-1)', k=0..mu/2) od; for ss from 1 to sigma do for mu from ss to maxmu+maxp+sigma-ss do K1[mu,ss] := (K1[mu+1,ss-1]-K1[mu-1,ss-1])/(2*mu+1) od od; for mu from sigma to maxmu+maxp do Kbar1[mu,0] := K1[mu,sigma] od; for p from 1 to maxp do Kbar1[sigma,p] := Kbar1[sigma+1,p-1]; for mu from sigma+1 to maxmu+maxp-p do Kbar1[mu,p] := ((mu+sigma+1)*Kbar1[mu+1,p-1]+(mu-sigma)*Kbar1[mu-1,p-1]) /(2*mu+1) od od end: printf("Build_Kbar loaded\n"); #============================================================================= ################# # Test Programs # ################# #----------------------------------------------------------------------------- # Slow, brute force evaluation of L mu sigma(p,alpha). # So slow its only value is verification of the other formulas. LLL := proc(mu,sigma,p,alpha) local T,t; T := QQ(mu,sigma,t)*t^p*exp(-t*alpha); evalf(((mu-sigma)!/(mu+sigma)!)*int(T,t=1..infinity)) end: printf("LLL loaded\n"); #----------------------------------------------------------------------------- Test_L := proc(mumax,sigma,p,alpha,tmax) local LX,LX1,LY1,mu; printf(" Test of L (Digits=%3.0f) for sigma=%1.0f, ",Digits,sigma); printf("p=%2.0f\n alpha=%6.3f, LY cutoff=%3.0f\n", p,alpha,tmax); printf(" mu Exact LX LY\n"); Digits := Digits+60; Build_LX(LX,mumax,sigma,p,alpha); Digits := Digits-60; Build_LX(LX1,mumax,sigma,p,alpha); Build_LY(LY1,mumax,sigma,p,alpha,tmax); for mu from sigma to mumax do printf("%3.0f %20.12e %20.12e %20.12e\n", mu,LX[mu,p],LX1[mu,p],LY1[mu,p]) od end: printf("Test_L loaded\n"); #----------------------------------------------------------------------------- Test_W := proc(mumax,sigma,p1,p2,alpha1,alpha2,tmax,jmax) local mu,WX,WX1,WY1; printf(" Test of W (Digits=%3.0f) for sigma=%1.0f, ",Digits,sigma); printf("p1=%2.0f, p2=%2.0f\n alpha1=%6.3f, alpha2=%6.3f\n", p1,p2,alpha1,alpha2); printf(" LY cutoff=%3.0f, WY cutoff=%3.0f\n",tmax,jmax); printf(" mu Exact WX WY\n"); Digits := Digits+60; Build_WX(WX,mumax,sigma,p1,p2,alpha1,alpha2); Digits := Digits-60; Build_WX(WX1,mumax,sigma,p1,p2,alpha1,alpha2); Build_WY(WY1,mumax,sigma,p1,p2,alpha1,alpha2,tmax,jmax); for mu from sigma to mumax do printf("%3.0f %20.12e %20.12e %20.12e\n", mu,WX[mu,p1,p2],WX1[mu,p1,p2],WY1[mu,p1,p2]) od end: printf("Test_W loaded\n"); #----------------------------------------------------------------------------- Test_WYY:= proc(mumax,alpha1,alpha2,tmax,jmax,pmax) local mu,WX,WYY; printf(" Test of WYY (Digits=%3.0f) for sigma=p1=p2=0\n",Digits); printf(" alpha1=%6.3f, alpha2=%6.3f\n",alpha1,alpha2); printf(" Cutoffs: LY=%3.0f j=%3.0f p=%3.0f\n", tmax,jmax,pmax); printf(" mu Exact WYY\n"); Digits := Digits+60; Build_WX(WX,mumax,0,0,0,alpha1,alpha2); Digits := Digits-60; Build_WYY(WYY,mumax,0,0,0,alpha1,alpha2,tmax,jmax,pmax); for mu from 0 to mumax do printf("%3.0f %20.12e %20.12e\n", mu,WX[mu,0,0],WYY[mu]) od end: printf("Test_WYY loaded\n"); #----------------------------------------------------------------------------- Test_L_recur := proc(mumax,alpha) local Lexact,Lup,Ldown,Ladj,mu; Digits := Digits+60; Build_LX(Lexact,mumax,0,0,alpha); Digits := Digits-60; L_recur_up(Lup,mumax,alpha); L_recur_dn(Ldown,mumax,alpha,0); L_recur_dn(Ladj,mumax,alpha,1); printf(" Test of recursive process for L (Digits=%3.0f) for sigma=p=0", Digits); printf(", alpha=%6.3f\n",alpha); printf(" mu L exact L up"); printf(" L down L adj\n"); for mu from 0 to mumax do printf("%3.0f %18.11e %18.11e %18.11e %18.11e\n",mu,Lexact[mu,0], Lup[mu],Ldown[mu],Ladj[mu]) od end: printf("Test_L_recur loaded\n"); #----------------------------------------------------------------------------- Test_Lp_recur := proc(mumax,p,alpha) local Lexact,Lvert,Lobl,mu; Digits := Digits+80; Build_LX(Lexact,mumax,0,p,alpha); Digits := Digits-80; Lp_recur_vert(Lvert,mumax,p,alpha); Lp_recur_obl(Lobl,mumax,p,alpha); printf("Test of recursive process for L (Digits=%3.0f) for sigma=0, p=%3.0f" ,Digits,p); printf(", alpha=%6.3f\n",alpha); printf(" mu L exact L vertical L oblique\n"); for mu from 0 to mumax do printf("%3.0f %18.11e %18.11e %18.11e\n",mu,Lexact[mu,p], Lvert[mu,p],Lobl[mu,p]) od end: printf("Test_Lp_recur loaded\n"); #----------------------------------------------------------------------------- Test_Ls_recur := proc(mumax,sigma,p,alpha) local LL0,LL1,mu; Digits := Digits+80; Build_LX(LL0,mumax,sigma,p,alpha); Digits := Digits-80; Ls_recur(LL1,mumax,sigma,p,alpha); printf("Test of recursive process for L (Digits=%3.0f) for sigma=%2.0f", Digits,sigma); printf(", p=%2.0f, alpha=%6.3f\n",p,alpha); printf(" mu L exact L recur\n"); for mu from sigma to mumax do printf("%3d%22.14e %22.14e\n",mu,LL0[mu,p],LL1[mu,sigma]) od end: printf("Test_Ls_recur loaded\n"); #----------------------------------------------------------------------------- Test_W_recur := proc(mumax,alpha1,alpha2) local Wex,WK,Wup,Wdn,Wadj,mu; Digits := Digits+60; Build_WX(Wex,mumax,0,0,0,alpha1,alpha2); Digits := Digits-60; Build_W_Kot(WK,mumax,alpha1,alpha2); Recur_W_up(Wup,mumax,alpha1,alpha2); Recur_W_dn(Wdn,mumax,alpha1,alpha2,0); Recur_W_dn(Wadj,mumax,alpha1,alpha2,1); printf(" Test of recursive processes for W (Digits =%3d)",Digits); printf(" for sigma=p1=p2=0\n", Digits); printf(" alpha1 =%6.3f, alpha2 =%6.3f\n", alpha1,alpha2); printf("mu W exact W Kotani W up"); printf(" W down W adj\n"); for mu from 0 to mumax do printf("%2d %16.9e %16.9e %16.9e %16.9e %16.9e\n",mu,Wex[mu,0,0], WK[mu,0,0],Wup[mu,mu],Wdn[mu,mu],Wadj[mu,mu]) od end: printf("Test_W_recur loaded\n"); #----------------------------------------------------------------------------- Test_Recur_Wp := proc(mumax,p1,p2,alpha1,alpha2) local mu,p,WX1,WYP; Digits := Digits+60; Build_WX(WX1,mumax,0,p1,p2,alpha1,alpha2); Digits := Digits-60; Recur_Wp(WYP,mumax,p1,p2,alpha1,alpha2); printf(" Test of recursion in p1 and p2 for W (Digits =%3d)\n",Digits); printf(" sigma=0, p1 =%2d, p2 =%2d",p1,p2); printf(", alpha1 =%6.3f, alpha2 =%6.3f\n",alpha1,alpha2); printf(" mu W exact W recursion\n"); for mu from 0 to mumax do printf("%3d %20.12e %20.12e\n",mu,WX1[mu,p1,p2],WYP[mu,p1,p2]) od end: printf("Test_Recur_Wp loaded\n"); #----------------------------------------------------------------------------- Test_Recur_Ws := proc(mumax,sigma,p1,p2,alpha1,alpha2) local WX1,WY1,mu; Digits := Digits+80; Build_WX(WX1,mumax,sigma,p1,p2,alpha1,alpha2); Digits := Digits-80; Recur_Ws(WY1,mumax,sigma,p1,p2,alpha1,alpha2); printf(" Test of recursion in sigma for W (Digits =%3d)\n",Digits); printf(" sigma=%2d, p1 =%2d, p2 =%2d",sigma,p1,p2); printf(", alpha1 =%6.3f, alpha2 =%6.3f\n",alpha1,alpha2); printf(" mu W exact W recursion\n"); for mu from sigma to mumax do printf("%3d %20.12e %20.12e\n",mu,WX1[mu,p1,p2],WY1[mu,p1,p2]) od end: printf("Test_Recur_Ws loaded\n"); #----------------------------------------------------------------------------- Check_TX := proc(mumax,alpha) local TX1,TX2,mu,nu; Digits := Digits+40; Build_TX(TX1,mumax,alpha); Digits := Digits-40; Build_TX(TX2,mumax,alpha); printf(" Qnu Pmu TX exact TX approx"); printf(" Digits =%3.0f alpha=%6.3f\n",Digits,alpha); for mu from 0 to mumax do for nu from 0 to mumax-mu do printf(" %2.0f %2.0f %20.12e %20.12e\n", nu,mu,TX1[nu,mu],TX2[nu,mu]) od od end: printf("Check_TX loaded\n"); #=============================end of file=====================================