function legendre_norm_Becke npoly=6; % This calculates the norm or orthogonality of two Maxwell polynomials clc format long e hold off load abspd.dat; n=90; alf=abspd(:,1); b=abspd(:,2); rtb=sqrt(b); xn=[]; acc1=[]; acc2=[]; % Start NQUAD LOOP for nquad=5:2:60 xn=[xn nquad]; rtbx=rtb; rtbx(nquad:n-1,:)=[]; %delete the bottom npoly to n rows alfx=alf; alfx(nquad+1:n,:)=[]; ptwt=legptwt(nquad); pt=ptwt(:,1); wt=ptwt(:,2); rmax=1; %Becke Mapping scale=rmax*(1-pt(nquad))/(1+pt(nquad)); scale=1; r=scale*(1+pt)./(1-pt); wtfcn=(r.*r).*exp(-r.*r); srwtfcn=sqrt(wtfcn); % Calcn of the Maxwell Polynomials gam=sqrt(pi/4); srgam=sqrt(gam); rtpi=sqrt(pi); b1=sqrt(4/(rtpi*(3*pi/8-1))); a1=rtpi*b1/2; p1=(2/sqrt(rtpi))*ones(1,nquad); p2=a1*r'-b1*ones(1,nquad); for n=2:npoly p3=(r'.*p2-alf(n)*p2)/rtb(n) - p1*rtb(n-1)/rtb(n); p1=p2; p2=p3; if n==4 p4=p2; else end end % End Calcn of Maxwell Polynomials npt=length(pt); zero=0*ones(1,nquad); %Becke Mapping maxpolyn1=2.d0*scale*wtfcn'.*p4.*p4./(1-pt').^2; maxpolyn2=2.d0*scale*wtfcn'.*p2.*p2./(1-pt').^2; integ1=sum(wt'.*maxpolyn1); integ2=sum(wt'.*maxpolyn2); err1=log10(abs(1-integ1)); err2=log10(abs(1-integ2)); acc1=[acc1 err1]; acc2=[acc2 err2]; plot(xn,acc1,'-^k',xn,acc2,'-ok','linewidth',1.6,'markerfacecolor','k','markersize',10); hold on axis([10 60 -9 1]) set(gca,'FontSize',32) set(gca,'Ytick',[-9:2:1],'linewidth',1.6) set(gca,'Xtick',[10:10:60],'linewidth',1.6) xlabel('$N$','Interpreter','LaTex','FontSize',32) ylabel('$\log_{10}|1-{\rm I}_{approx}|)$','Interpreter','LaTex','FontSize',32 ) str1={'Boyd-Becke map'}; str2={'${x=\frac{r+1}{r-1}}$'}; text(12,-6.8,str1,'Interpreter','latex','fontsize',32) text(12,-8,str2,'Interpreter','latex','fontsize',32) set(gcf, 'Units','centimeters','Papersize',[36,36]) set(gcf, 'Units','centimeters','Position',[3 3 24 20]) end function ptwt = legptwt(n) xn=[0:1:n-1]; xnsq=xn.*xn; b=xnsq./(4*xnsq-1); rtb=sqrt(b); rtb(1)=[]; t=diag(rtb,-1)+diag(rtb,1); [f,lambda]=eig(t); pt=diag(lambda); for i=1:n wt(i)=2*f(1,i)^2; end ptwt=[pt,wt'];