function rt_chandra(n) %Legendre quadrature poits and weights [ptwt]=legptwt(n); pt=ptwt(:,1); wt=ptwt(:,2); % Use the positive points pt(1:n/2)=[]; wt(1:n/2)=[]; dx=0.1; jjmax=n/2-1; % Calculate the jj-th eigenvalue of the Radiative Transfer equation for jj=1:jjmax % Grid x=[1/pt(jj+1)+.001:.001:1/pt(jj)-.001]; % Search for the root of the function, f(x), below f=0; for jk=1:n/2 f=f+wt(jk)./(1-x.^2*pt(jk)^2); end f=f-1; % Function is plotted to see that the root has been bracketed plot(x,f,'-k') % Root bracketed between xl (x-"lower") and xu (x-"upper") nx=length(x); nf=length(f); t1=f(1:nf-1); t2=f(2:nf); tt=t1.*t2; indx=find(tt<0)'; xl=x(indx); xu=x(indx+1); fl=f(indx); fu=f(indx+1); root=(fl*xu+fu*xl )/(fl+fu); % Iteration to calculate the root for kk=1:30 xr=(xl+xu)/2; % Print successive roots in the iteration % fprintf('%i %20.12f\n',kk,xu) fu=0; for jk=1:n/2 fu=fu+wt(jk)./(1-xu.^2*pt(jk)^2); end fu=fu-1; fl=0; for jk=1:n/2 fl=fl+wt(jk)./(1-xl.^2*pt(jk)^2); end fl=fl-1; fr=0; for jk=1:n/2 fr=fr+wt(jk)./(1-xr.^2*pt(jk)^2); end fr=fr-1; if(fu*fr) < 0 xl=xr; else xu=xr; end if(fl*fr) < 0 xu=xr; else xl=xr; end end % Store the eigenvalues eig(jjmax-jj+1)=xu; end display('Eigenvalues') for i=1:n/2-1 fprintf('%2i %16.8f\n',i,eig(i)) end % Calculate the extrapolation length one=linspace(1,1,n/2)'; mat=[]; for k=1:(n-2)/2 col=(1./(1-eig(k)*pt)); mat=[mat col]; end mat=[mat one]; sol=mat\pt; display('Extrapolation length') nlast=(n-2)/2; fprintf('%10.5f\n', sol(n/2)) % Legendre quadrature points and weights function [ptwt] = legptwt(n) xn=[1:1:n]; xnsq=xn.*xn; b=xnsq./(4*xnsq-1); rtb=sqrt(b); rtb(n)=[]; t=diag(rtb,-1)+diag(rtb,1); [f,lambda]=eig(t); pt=diag(lambda); wt=2*f(1,:).^2; ptwt=[pt,wt'];