%Chandrasekhar_energy_calculations
% Alternative method to estimate Chandrasekhar mass, and mass vs radiaus
% relationship for a white dwarf of mass << Chandrasekhar mass. Consider
% energy E = (1/3)*degeneracy pressure x volume + GPE. Assumes constant
% density and therefore constant degeneracy pressure throughout star. But
% GPE is using Betti-Ritter so assumes a polytrope!
%
% LAST UPDATED by Andy French Jan 2026

function Chandrasekhar_energy_calculations

%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;
mu = 2*mp;

%Fontsize for graphs
fsize = 18;

%Set mass of white dwarf
M = 2*Msun;

%

%Classical model (polytrope n=3/2)
Ac = (3/2)*(4/3)*pi*(1/20)*( (3/pi)^(2/3) )*( h^2 ) *...
    ( 1/ ( me*( mu^(5/3) ) ) )*( ( M / (4*pi/3) )^(5/3) );
Bc = (6/7)*G*M^2;
Rmin = 2*Ac/Bc;
Emin = Ac/Rmin.^2 - Bc/Rmin;

%Check algebraic expressions for radius R such that dE/dR = 0

%Exact expression
Rmin_check = (7/5)*( ( 81/(8192*(pi^4)) )^(1/3) )*...
    ( (h^2) / ( G*me*(mu^(5/3)) ) )*( M^(-1/3) );
d1 = Rmin_check/Rmin;

%Approximate expression
Rmin_check_check = 4438.5e3 * ( M/(1.44*Msun) )^(-1/3);
d2 = Rmin_check_check/Rmin;


%

%Relativistic parameters (polytrope n=3)
Ar = (3/2)*(4/3)*(1/8)*pi*((3/pi)^(1/3))*h*c*(1/(mu^(4/3)))*( ( M / (4*pi/3) )^(4/3) ) ;
Br = (3/2)*G*M^2;

%Estimate Chandrasekhar mass i.e. mass such that Ar = Br
M_ch_by_Msun =  (2^(-11/2) )*(3/pi)*( ( h*c/(G*(mu^(4/3))) )^(3/2) )/Msun;
MM = M_ch_by_Msun*Msun;
Ar_check = (3/2)*(4/3)*(1/8)*pi*((3/pi)^(1/3))*h*c*(1/(mu^(4/3)))*( ( MM / (4*pi/3) )^(4/3) ) ;
Br_check = (3/2)*G*MM^2;
d3 = Ar_check/Br_check;
disp([' Chandresekhar mass estimated by energy method is ',num2str(M_ch_by_Msun,4), ' solar masses.']);

%

%Define radius
R = linspace(0.2*Rearth,2*Rearth,1000);

%Calculate total energy /J
Ec = Ac./R.^2 - Bc./R;
Er = (Ar-Br)./R;

%Plot classical vs relativistic energy expressions
plot( R/Rearth, Ec, R/Rearth, Er,'linewidth',2 ); box on; grid on; hold on; 
legend({'Classical','Relativistic'},'fontsize',fsize); axis manual
ylimits = get(gca,'ylim'); plot(Rmin*[1,1]/Rearth,ylimits,'k--','linewidth',2 );
xlimits = get(gca,'xlim'); plot(xlimits,Emin*[1,1],'k--','linewidth',2 );
set(gca,'fontsize',fsize);
xlabel('R / R_{earth}'); ylabel('Total energy /J'); 
title(['Energy of white dwarf: M/M_{sun} = ',num2str(M/Msun)]);
print(gcf,['White dwarf energy vs R M=',strrep(num2str(M/Msun,3),'.','p'),'.png'],'-dpng','-r300'); close(gcf);

%End of code
