%star_homology
% Calculates scaled (mass / solar mass) powers of (solar scaled)
% luminosity, radius, temperature and density for a quasi-static Main
% Sequence star. Ignore effect of radiation pressue.
%
% LAST UPDATED by Andy French Dec 2025

function star_homology

%Determine density and pressure sacling factors (SI units)
Msun = 1.989e30;
Rsun = 696340e3;
Tsun = 5778;
kB = 1.381e-23;
mp = 1.673e-27;
rho_scale = Msun/( (4/3)*pi*Rsun^3 );
P_scale = rho_scale*kB*Tsun/mp;

%

fusion_type = 'pp chain low temp';

%Energy per unit mass powers of density and temperature
alpha = 1;  beta = 4;  %pp

%Opacity parameterization of density and 1/temperature
y = 1; z = 3.5;  %Kramer's opacity

%Solve for radiative energy transport
[A,B,C,D,dlogL_by_minus_dlogT] = disp_homological_powers( fusion_type, alpha, beta, y, z );

%

fusion_type = 'pp chain low temp but constant opacity';

%Energy per unit mass powers of density and temperature
alpha = 1;  beta = 4;  %pp

%Opacity parameterization of density and 1/temperature
y = 0; z = 0;  %Fixed opacity

%Solve for radiative energy transport
[A,B,C,D,dlogL_by_minus_dlogT] = disp_homological_powers( fusion_type, alpha, beta, y, z );

%

fusion_type = 'CNO chain high temp';

%Energy per unit mass powers of density and temperature
alpha = 1;  beta = 17;  %CNO

%Opacity parameterization of density and 1/temperature
y = 0; z = 0;  %Constant opacity

%Solve for radiative energy transport
[A,B,C,D,dlogL_by_minus_dlogT] = disp_homological_powers( fusion_type, alpha, beta, y, z );

%

fusion_type = 'triple alpha higher temp';

%Energy per unit mass powers of density and temperature
alpha = 2;  beta = 40;  %Triple alpha

%Opacity parameterization of density and 1/temperature
y = 0; z = 0;  %Constant opacity

%Solve for radiative energy transport
[A,B,C,D,dlogL_by_minus_dlogT] = disp_homological_powers( fusion_type, alpha, beta, y, z );

%%

%Prints homological powers in an easy to read format
function [A,B,C,D,dlogL_by_minus_dlogT] = disp_homological_powers( fusion_type, alpha, beta, y, z )

%Solve for radiative energy transport
M_coeffs_ABCD_rad =...
    ( [0,3,0,1 ; 1,0,-beta,-alpha ; 0,4,1,1 ; -1, 4, z+4, -y ]^-1 )*[1;1;2;1];
A = M_coeffs_ABCD_rad(1);
B = M_coeffs_ABCD_rad(2);
C = M_coeffs_ABCD_rad(3);
D = M_coeffs_ABCD_rad(4);

%Predict gradient of Hertzsprung-Russell diagram i.e. d(logL)/d(-logT)
dlogL_by_minus_dlogT = -4/( 1-2*B/A );

%Display results to command window
disp(' ');
disp('*****************************************')
disp(fusion_type);
disp( ['L proportional to M ^ ',num2str(A)] );
disp( ['R proportional to M ^ ',num2str(B)] );
disp( ['T proportional to M ^ ',num2str(C)] );
disp( ['core density proportional to M ^ ',num2str(D)] );
disp( ['HR diagram gradient = ',num2str(dlogL_by_minus_dlogT)] );
disp( ['A - 2B -4C = ',num2str(A-2*B-4*C)] );
disp('*****************************************')

%End of code