%Eddington model of a star. (Arthur Eddington 1882-1944).
% Assumes a constant ratio of gas to radiation pressure within a star.
% Beta is the fraction of total pressure (gas + radiation) that is gas
% pressure.
%
% LAST UPDATED by Andy French Dec 2025

function edddington

%Fontsize for graphs
fsize = 18;

X_H = 0.747; %Hydrogen mass fraction (for Sun)
X_He = 0.236; %Helium mass fraction (for Sun)

%Determine average molecular mass (in kg) per electron ionized
mu_by_mp = 4/( 6*X_H + X_He + 2);

%Check beta using solar parameters
beta_sun = eddington_quartic( 1, mu_by_mp );
disp(' ');
disp([' Solution to Eddingtons quartic for M = Msun, mu/mp = ',...
    num2str(mu_by_mp,3),': beta = ',num2str(beta_sun)]);
disp(' ');

%Check beta for M = 150, mu_by_mp=2
beta_sun = eddington_quartic( 150, 2 );
disp(' ');
disp([' Solution to Eddingtons quartic for M = 150Msun, mu/mp = 2: beta = ',num2str(beta_sun)]);
disp(' ');

%Star mass in solar masses
M_by_Msun = linspace(0.5,200,1000);

%Molecular mass in proton masses
N = 10;
mu_by_mp = linspace(0.5,2,N);

%

%Calculate beta using Eddington model
for n=1:N
    for m=1:length(M_by_Msun)
        [beta(m,n),psi] = eddington_quartic( M_by_Msun(m), mu_by_mp(n) );
    end
    %Build up legend for graphs
    lgnd{n} = ['\mu/m_{p} = ',num2str(mu_by_mp(n),2)] ;
end

%Plot beta vs M/Msun for various mu_by_mp values
RGB = line_plot_colours(N);
M_by_Msun = repmat( M_by_Msun.',[1,N]);
hold on; warning off
for n=1:N
    plot( M_by_Msun(:,n), beta(:,n),'-','linewidth',2,'color',RGB(n,:));
end
set(gca, 'fontsize',fsize); grid on; box on; legend(lgnd,'fontsize',8);
xlabel(' M / M_{sun}'); ylabel('\beta');
title('Eddington fraction \beta of total pressure that is gas pressure');
print(gcf,'eddington beta vs M by Msun.png','-r300','-dpng'); close(gcf);

%

%Plot Eddington's quartic for a variety of M and mu values
N = 10;
M_by_Msun = 10.^linspace(-1,2,N);  %Log-spaced masses
mu_by_mp = linspace(0.5,2,N); %Full range of mu
RGB = line_plot_colours(N);
hold on; plot([0,1],[0,0],'k-','linewidth',2); %Plot y = 0;
warning off
for n=1:N
    %Solve Eddington's quartic
    [bb,psi] = eddington_quartic( M_by_Msun(n), mu_by_mp(n) );
    b = linspace(0,1,1000);  %Beta values
    y = ( M_by_Msun(n)^2 )*b.^4 + psi*b*( mu_by_mp(n)^(-2) ) -...
        psi*( mu_by_mp(n)^(-2) ); %Eddington's Quartic
    
    %Plot Eddinton's quartic & roots
    p(n) = plot(b,y,'-','linewidth',2,'color',RGB(n,:));
    plot(bb,0,'*','color',RGB(n,:));
    
    %Construct legend entry
    lgnd{n} = ['M/M_{sun}=',num2str(M_by_Msun(n),3),...
        ', \mu/m_{p} = ',num2str(mu_by_mp(n),2)] ;
end
set(gca, 'fontsize',fsize); grid on; box on; ylim([-20,20]); legend(p,lgnd,'fontsize',8);
ylabel('Eddington quartic'); xlabel('\beta'); title('Eddington quartic');
print(gcf,'eddington quartic.png','-r300','-dpng'); close(gcf);

%%

%Solve Eddington's quartic numerically for the fraction beta of total pressure that is gas
%pressure (rather than radiation pressure).
function [beta,psi] = eddington_quartic( M_by_Msun, mu_by_mp )

%Solar parameters and universal constants in SI units
Msun = 1.988416e30;
Rsun = 696340e3;
Rearth = 6371e3;
kB = 1.380649e-23;
G = 6.67430e-11;
sigma = 5.670374419e-8;
mp = 1.67262192e-27;
me = 9.1093837e-31;
h = 6.62607015e-34;
c = 299792458;

%Mass in kg of gas particles per ionized electron in star
mu = mu_by_mp*mp;

%Lane-Emden parameters for n=3 polytrope
A = 2.02;  Y0 = 6.90;

%Solve Edington Quartic for fraction beta of total pressure that is gas
%pressure (rather than radiation pressure).
phi = (12/pi)*( A^2 )*( kB^4 )*c/( sigma *(G^3) * (mp^4) *( Msun^2 ) );
psi = sqrt(phi);
phi = psi/(mu_by_mp^2);
betas = roots( [M_by_Msun^2, 0,0,phi,-phi] );

%Find root which is +ve and real
betas( abs( imag(betas) )> eps ) = [];
betas( betas<0 ) = [];
beta = betas;

%Check psi calculation from polytrope n=3 equation for M
a = ( 3 *(kB^4)*c/( 4*sigma * mp^4 ) )^(1/3);
psi_check = ( ( a/(pi*G) )*(( 1/(3*A) )^(-2/3))*( (4*pi/3)^(2/3) ) )^(3/2);
psi_check = psi_check/Msun;
psi_check2 = (4*A/sqrt(pi))*( ( a/G )^(3/2) )/Msun;

%%

%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