%chandrasekhar
% Performs calculations for the limiting mass of a white dwarf stellar
% remnant, and the radius vs mass relationship for white dwarfs of less
% than this critical mass (when electrons are not relativistic).
%
% LAST UPDATED by Andy French Jan 2026

function chandrasekhar

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

%

%% Calculation for Sun using Lane-Emden polytropic model %%
disp(' ')
disp('*******************************************************************')
disp('Lane Emden & Ideal Gas calculation for the Sun')

rho0 = 53.95*Msun/((4/3)*pi*Rsun^3);
mu=0.6*mp;
T0 = ( pi*G*mu*rho0/kB )*( Rsun/6.89 )^2;
T0_check = (3*53.95/(4*6.89^2))*G*Msun*mu/(kB*Rsun);
disp([' Solar core temperature = ',num2str(T0,3), ' K.'] );

%Check results in correct masses and radii
alpha = sqrt( kB*T0*(3+1)/(4*pi*G*mu*rho0) );
M_by_Msun = 4*pi* (alpha^3) *rho0 * 2.02 /Msun;
R_by_Rsun = alpha*6.89/Rsun;
disp('*******************************************************************')


%% Calculate average spacing of carbon nuclei in white dwarf at 1.44Msun %%
disp(' ')
disp('*******************************************************************')
disp('Lane Emden & Ideal Gas calculation for a white dwarf at the Chandrasekhar limit')

%Mass per ionized electron in white dwarf (assume metalicity of unity, i.e.
%mostly C and O atoms).
mu = 2*mp; 
M = 1.44*Msun;
T0 = 9*6e8;
rho0 = ( 16*(pi^2) * (2.02^2) / (M^2) )*( (kB*T0/(pi*G*mu ) )^3 );
x = ( (9*mp*M^2) *( ( pi*G*mu/(kB*T0) )^3 )/ (16*(pi^3)*2.02^2) )^(1/3);
x_check = ( 12*mp/( (4/3)*pi*rho0 ) )^(1/3);
alpha = sqrt( kB*T0*(3+1)/(4*pi*G*mu*rho0) );
M_by_Msun = 4*pi* (alpha^3) *rho0 * 2.02 /Msun;
R_by_Rsun = alpha*6.89/Rsun;
disp([' Assume core temperature is T0 = ',num2str(T0,3), ' K']);
disp([' White dwarf radius/solar radius = ',num2str(R_by_Rsun,3) ]);
disp([' White dwarf core density = ',num2str(rho0,3), ' kg/m^3'])
disp([' Interatomic spacing x = ',num2str(x,3), 'm']);
disp('*******************************************************************')

%

%% Calculate average spacing of carbon nuclei in white dwarf at 1.44Msun %%
disp(' ')
disp('*******************************************************************')
disp('Lane Emden & Ideal Gas calculation for a white dwarf at the Chandrasekhar limit #2')

%Mass per ionized electron in white dwarf (assume metalicity of unity, i.e.
%mostly C and O atoms).
mu = 2*mp; 
M = 1.44*Msun;
rho0 = 54.16*1.44*Msun/( (4/3)*pi*( (7762e3)^3 ) );
x = ( 12*mp/( (4/3)*pi*rho0 ) )^(1/3);
alpha = 7762e3/6.89;
T0 = 4*pi*G*mu*rho0*(alpha^2)/( kB*(3+1) );
disp([' White dwarf core density = ',num2str(rho0,3), ' kg/m^3'])
disp([' White dwarf core temperature is T0 = ',num2str(T0,3), ' K']);
disp([' Interatomic spacing x = ',num2str(x,3), 'm']);
disp('*******************************************************************')

%

%% Relativistic limit i.e. Chandrasekhar mass %%

%Lane-Emden parameters
n = 3; Y0 = 6.90;  A = 2.02; rho_c_by_rho_av_rel = 54.16;

%Polytropic constant a from degeneracy pressure P = a * rho^(1+1/n)
a = (1/8)*( (3/pi)^(1/3) )*h*c/ ( mu^(4/3) );
k = a*(n+1)/(4*pi*G);
k = k*( (4*pi/3)^(1-(1/n)) );
k = k*( Y0^2 );
k = k*( ( (Y0^3)/(3*A) )^( 1/n - 1 ) );
M_ch = k^(n/(n-1));
M_ch_check =  ( A*sqrt(3)/( pi * sqrt(32) ) )*( ( h*c/( G*( mu^(4/3) ) ) )^(3/2) );
M_ch_by_Msun_check = M_ch_check/Msun;
M_ch_by_Msun = M_ch/Msun;
disp(' ')
disp('************************************************************************************')
disp([' Chandrasekhar mass = ',num2str(M_ch_by_Msun),' solar masses.'])

%

%% Non relativistic limit %%

%Lane-Emden parameters
n = 3/2; Y0 = 3.65; A = 2.71; rho_c_by_rho_av_classical = 5.99;

%Polytropic constant a from degeneracy pressure P = a * rho^(1+1/n)
a = (h^2) /( me * ( mu^(5/3) ) );
a = (1/20)*( (3/pi)^(2/3) )*a;

%Calculate radius of white dwarf 
%R (in km) = ( R_ch in km )* (M/M_ch)^(-1/3)
k = a*(n+1)/(4*pi*G);
k = k*( ( (Y0^3)/(3*A) )^( (1/n) -1) );
k = k*( ( 4*pi/3)^(1 - 1/n ) );
k = k*(Y0^2);
k = k^(n/(3-n));
M_power = (1-n)/(3-n);
k = k*( (M_ch)^M_power );
R_ch_in_km = k/1000;
disp([' White dwarf radius R (in km) = ( ',num2str(R_ch_in_km),' km )* (M/M_ch)^(-1/3)']);

f = 9/( (pi^4)*(2^13) ); f = f^(1/3); 
f = f*( A^(1/3) )*Y0;
R_ch_check = (1e-3)*f*( (h^2)/( G*me*( mu^(5/3) )) )*( M_ch^(-1/3) );

%Calculate average and core density of (non-relativistic) white dwarf at
%M_ch
rho_av = M_ch/ ( (4/3)*pi*( ( R_ch_in_km*1000 )^3 ) );
rho_c_classical = rho_av*rho_c_by_rho_av_classical;
rho_c_rel= rho_av*rho_c_by_rho_av_rel;
disp([' White dwarf density at core assuming classical n=3/2 model: ',num2str(rho_c_classical,3), ' kg/m^3' ] )
disp([' White dwarf density at core assuming relativistic n=3 model: ',num2str(rho_c_rel,3), ' kg/m^3' ] )

%(Mass/solar mass) x (volume/earth volume) product for white dwarf
M_times_V = M_ch_by_Msun*( R_ch_in_km / ( Rearth/1000) )^3;
disp([' White dwarf mass in solar masses x white dwarf volume in Earth volumes = ',...
    num2str(M_times_V) ]);
disp('************************************************************************************')
disp(' ')

%

%Plot Radius vs white dwarf mass (classical limit)
M_by_Msun = linspace(0,1.44,1000);
R_in_km = 7762*( M_by_Msun/1.44 ).^(-1/3);
plot( M_by_Msun,R_in_km/1000,'r-','linewidth',2 );
hold on; xlimits = get(gca,'xlim'); ylimits = get(gca,'ylim');
plot( xlimits,7.762*[1,1],'k--','linewidth',2 );
plot( 1.44*[1,1],ylimits,'k--','linewidth',2 );
plot(1.44,7.762,'k*','markersize',10);
plot( M_by_Msun,R_in_km/1000,'r-','linewidth',2 );
set(gca,'fontsize', 18); xlabel('M/M_{sun}'); ylabel('R /1000km');
title('Classical limit of white dwarf radius vs mass');
box on; grid on;
print(gcf,'white dwarf R vs M.png','-dpng','-r300'); close(gcf);

%End of code
