function legendre_norm_MK clc npoly=6; % This calculates the norm or orthogonality of two Maxwell polynomials 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 npoly to n rows alfx=alf; alfx(nquad+1:n,:)=[]; ptwt=legptwt(nquad); pt=ptwt(:,1); wt=ptwt(:,2); rmax=8; %Mura-Knowles-map! r=(1+pt).^2./(4-(1+pt).^2); %plot(r,pt,'-k') 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: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); %Mura-Knowles map! drdx=8*(1+pt)./(4-(1+pt).^2).^2; maxpolyn1=drdx'.*wtfcn'.*p4.*p4; maxpolyn2=drdx'.*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','linewidth',1.6,'markerfacecolor','k','markersize',10); hold on axis([10 60 -6 0]) set(gca,'FontSize',32) set(gca,'Ytick',[-6:2: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) set(gcf, 'Units','centimeters','Papersize',[36,36]) set(gcf, 'Units','centimeters','Position',[3 3 24 20]) %linear map str1={'Mura-Knowles map'}; str2={'$x=2\sqrt{r/(r+1)}-1$'}; text(12,-4,str1,'Interpreter','latex','fontsize',32) text(12,-4.8,str2,'Interpreter','latex','fontsize',32) 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'];