function rate_O_O2(s) % Convergence of the scaled Laguerre quadratures for the O-O2 % dissociative reaction rate coefficient clc format long e % "Exact" value is calculated with 800 quadrature pointn rexact= 1.25756101; temp=6000; kt=0.8617e-04*temp; yzero=14.5/kt; nn=[]; acc2=[]; % Calculate the rate coefficient for N quadrature points for n=2:2:50 nn=[nn n]; % Laguerre quadratures for alfa = 1 [pt,wt] = lagptwt2(n,1); [sig]=sigma_O_O2((s*pt+yzero)*kt); rate1=s^2*sum((wt.*exp(-(s-1)*pt')).*sig); [pt,wt] = lagptwt2(n,0); [sig]=sigma_O_O2((s*pt+yzero)*kt); rate2=s*sum((wt.*exp(-(s-1)*pt')).*sig); rate=1e12*(exp(-yzero)*rate1+yzero*exp(-yzero)*rate2); acc=log10(abs(1-rate/rexact)); acc2=[acc2 acc]; fprintf('%i %20.12f %13.5e\n',n,rate,acc) end % Plot the accuracy versus N plot(nn,acc2,'-ok','linewidth',1.6,'markersize',7,'markerfacecolor','k') hold on % Run again for different scale factors - change plotting symbol % Format the graph axis([0 50 -6 -2]) set(gca,'FontSize',20) set(gca,'Ytick',[-6:1:-2],'linewidth',1.6) set(gca,'Xtick',[0:10:50],'linewidth',1.6) ylabel('$Accuracy$','Interpreter','LaTex','FontSize',24) xlabel('$N$','Interpreter','LaTex','FontSize',24) % Reactive cross section function [sig]=sigma_O_O2(e) n=length(e); for i=1:n sig(i)=4.51*(e(i)-14.5)^1.03/(0.21+e(i)^1.31); end % Laguerre quadratures function [pt,wt] = lagptwt2(n,alf) format long e xn=[0:1:n-1];a=2*xn+alf+1; b=(xn.*(xn+alf)); rtb=sqrt(b); rtb(1)=[]; t=diag(rtb,-1)+diag(a)+diag(rtb,1); [f,lambda]=eig(t); pt=diag(lambda); wt=gamma(alf+1)*f(1,:).^2; ptwt=[pt,wt'];