function legendre_norm_linear
% 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
m=2;
for nquad=5:2:60
    xn=[xn nquad];
rtbx=rtb; rtbx(nquad:n-1,:)=[];  %delete the bottom nquad to n rows
alfx=alf; alfx(nquad+1:n,:)=[];
% Legendre quadrature points and weights
ptwt=legptwt(nquad); pt=ptwt(:,1); wt=ptwt(:,2); 
rmax=8; r=0.5*(pt+1)*rmax;
wtfcn=(r.*r).*exp(-r.*r); srwtfcn=sqrt(wtfcn);
% Calcn of 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:6
    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);
maxpolyn1=0.5*rmax*wtfcn'.*p4.*p4; maxpolyn2=0.5*rmax*wtfcn'.*p2.*p2;
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','markerfacecolor','k','markersize',12);
axis([10 60 -16 0])
set(gca,'FontSize',32)
set(gca,'Ytick',[-16:4:0],'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)
legd1=legend('$<M_4^2>$','$<M_6^2>$','Location','NorthEast');
set(legd1,'Interpreter','latex','fontsize',32);
set(gcf, 'Units','centimeters','Papersize',[36,36])
set(gcf, 'Units','centimeters','Position',[3 3 24 20])
%linear map
str1={'Linear map'};
str2={'$r_{max}=8$'};
text(12,-13,str1,'Interpreter','latex','fontsize',32)
text(12,-14.5,str2,'Interpreter','latex','fontsize',32)
end
% Legendre quadrature points and weights
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);
wt=2*f(1,:).^2;
ptwt=[pt,wt'];