%lane_emden
% Numeric solutions to the Lane-Emden differential equation, which occurs
% in models of quasi-static stellar structure for a polytropic star of
% pressure vs density profile P = k * rho^(1+1/n) where n is the polytropic
% index.
%
% LAST UPDATED by Andy French Dec 2025

function lane_emden

%Define polytropic index n
n = [1,1.5,2,3,4,5];

%Define Y array increments and maximum
dY = 0.0001; Ymax = 20;

%Graph fontsize
fsize = 18;

%Legend string
for k=1:length(n)
    lgnd_str{k} = ['n = ',num2str(n(k))];
end

%%

%Initialize solution arrays (one row for each n value)
Y = 0:dY:Ymax;
YY = repmat( Y.',[1,length(n)] );
QQ = ones(size(YY));
dQQ_by_dYY = ones(size(YY));

%Solve using Euler method
N = length(n);
for k=1:N
    [Q,dQ_by_dY,Y0(k),K(k)] = lane_emden_solver( n(k),dY,Ymax);
    QQ(:,k) = Q(1,:);
    dQQ_by_dYY(:,k) = dQ_by_dY(1,:);
end

%Print table of outputs
output = [n; Y0; K; (Y0.^3)./(3*K) ]

%Plot outputs
RGB = line_plot_colours(N); hold on;
for n=1:N
    plot(YY(:,n),QQ(:,n),'linewidth',2,'color',RGB(n,:))
end; ylim([-0.5,1]); xlim([0,Ymax]); grid on; box on;
set(gca,'fontsize',fsize);
xlabel('\psi'); ylabel('\theta'); title('Solution to Lane-Emden equation');
legend(lgnd_str,'fontsize',fsize);
print(gcf,'lane-emden theta vs psi.png','-dpng','-r300'); close(gcf);
hold on;
for n=1:N
    plot(YY(:,n),dQQ_by_dYY(:,n),'linewidth',2,'color',RGB(n,:))
end; ylim([-0.5,1]); xlim([0,Ymax]); grid on;xlim([0,Ymax]); ylim([-1,0.5]);
grid on; box on;
set(gca,'fontsize',fsize);
xlabel('\psi'); ylabel('d\theta/d\psi'); title('Solution to Lane-Emden equation');
legend(lgnd_str,'fontsize',fsize,'location','best');
print(gcf,'lane-emden dthetabydpsi vs psi.png','-dpng','-r300'); close(gcf);

%%

%Euler method solver for Lane-Emden equation
%k is -Y^2 * dQ_by_dY evaluated ay Y0, i.e. when Q=0
function [Q,dQ_by_dY,Y0,k] = lane_emden_solver(n,dY,Ymax)

%Determine Y array
Y = 0:dY:Ymax;

%Analytic cases (n=2 has a series expansion, and n=5 is solvable when Y>0,
%but I have ignored these!
%https://en.wikipedia.org/wiki/Lane%E2%80%93Emden_equation)
if n==0
    Q = 1 - (1/6)*Y.^2;
    dQ_by_dY = -(1/3)*Y;
    Y0 = sqrt(6);
    k = 2*sqrt(6);
elseif n==1
    Q = sin(Y)./Y;
    dQ_by_dY = ( Y.*cos(Y) - sin(Y) )./(Y.^2);
    Y0 = pi;
    k = pi;
else
    %Use default for n=5. This will be reset if i0 is not NaN i.e. Q<0
    E = zeros(size(Y)); Q = ones(size(Y)); Y0calc=1; i0 = NaN;
    Y0 = inf; k = 1.7;
    
    %Apply Euler method to solve
    for i = 1:(length(Y)-1)
        E(i+1) = E(i) + ( Y(i)^2 ) * ( Q(i)^n ) * dY;
        if Y(i)==0;
            Q(i+1) = 1;
        else
            Q(i+1) = Q(i) - E(i)*dY/( Y(i)^2 );
        end
        if ( Q(i+1)<0 ) && (Y0calc==1)
            i0 = i; Y0calc = 0;  %This corresponds to the star radius (if it has one!)
        end
    end
    dQ_by_dY = -E./( Y.^2 );
    if ~isnan(i0)
        Y0 = Y(i0); k = - Y0^2 * dQ_by_dY(i0);
    end
end

%%

%N line plot colours from jet colormap
function RGB = line_plot_colours(N)
%Determine plot line colours by interpolating 'jet' colormap
RGB = zeros(N,3); cmap = colormap('jet'); close(gcf);
red = cmap(:,1).'; green = cmap(:,2).'; blue = cmap(:,3).';
dim = size(cmap); nn = linspace(1,N,dim(1));
for n=1:N
    RGB(n,1) = interp1( nn, red, n );
    RGB(n,2) = interp1( nn, green, n );
    RGB(n,3) = interp1( nn, blue, n );
end

%End of code