% Construct Hertzsprung-Russell (HR) diagram from star catalogue 
% of about 45,000 'nearby' stars.
% https://astrosci.scimuze.com/stellar_data.htm
%
% LAST UPDATED by Dr Andrew French Dec 2025

function make_HR_diagram

%Solar system parameters
Msun = 1.989e30;  %Solar mass /kg
Rsun = 696340;  %Solar radius in km
Tsun = 5772; %Solar effective temperature in K
AU = 1.496e11; %Astronomical unit (approx Earth-Sun distance) in m
Lsun = 3.846e26; %Solar luminosity (W)
Ly = 9.461e15;  %Light-year in m
parsec = 3.262*Ly;  %Parsec in m

%Set fontsize
fsize = 18;

%Set HR diagram limits for x = -log(T/Tsun) and y = log(L/Lsun)
xlimits = [-0.5,0.4]; ylimits = [-6,7];

%Define masses of stars which represent the Main Sequence
N = 20; M = Msun*10.^(linspace( log10(0.1), log10(20), N ));
NN = 2000; MM = Msun*10.^(linspace( log10(0.01), log10(20), NN ));

%

%Determine star luminosity (W) using emprical relationship determined from
%direct measurements from binary star systems
L = Lsun*mass_luminosity(M/Msun); LL = Lsun*mass_luminosity(MM/Msun);

%Perform a y = m*x line of best fit to L/Lsun vs (MM/Msun)^3
[yfit,xfit,r,m,dm,yupper,ylower,s] = bestfit( (MM/Msun).^3,LL/Lsun );

%Obtain star effective temperature (in K) using MS empirical fit
T = Tsun*T_from_L_MS(L/Lsun); TT = Tsun*T_from_L_MS(LL/Lsun);

%Load star catalogue
load star_data
N = length(stars.ID);

%Determine x,y coordinates of HR diagram
x = -log10( stars.T_K/Tsun ); y = log10( stars.L_by_Lsun );
xx = -log10( T/Tsun ); yy = log10( L/Lsun );
xxx = -log10( TT/Tsun ); yyy = log10( LL/Lsun );

%

%Plot Luminosity vs star mass
figure; plot(MM/Msun,LL/Lsun,'r-','linewidth',2); hold on;
plot( MM/Msun, 5.78*(MM/Msun).^3 , 'r--','linewidth',2);
legend({'Empirical','L/L_{sun} = 5.78(M/M_{sun})^3'});
plot(M/Msun,L/Lsun,'k*','linewidth',1);
set(gca,'fontsize',fsize); grid on; box on;
xlabel('M/M_{sun}'); ylabel('L/L_{sun}'); title('Luminosity vs mass for MS stars');
print(gcf,'L vs M for MS.png','-dpng','-r300'); close(gcf);

%Plot log10( Luminosity ) vs star mass
figure; plot(MM/Msun,log10(LL/Lsun),'r-','linewidth',2); hold on;
plot(M/Msun,log10(L/Lsun),'k*','linewidth',1);
set(gca,'fontsize',fsize); grid on; box on;
xlabel('M/M_{sun}'); ylabel('log_{10}( L/L_{sun} )'); title('log_{10}(L/L_{sun}) vs mass for MS stars');
print(gcf,'log10 L vs M for MS.png','-dpng','-r300'); close(gcf);

%Plot effective temperature vs star mass
figure; plot(MM/Msun,TT/Tsun,'r-','linewidth',2); hold on;
plot(M/Msun,T/Tsun,'k*','linewidth',1);
set(gca,'fontsize',fsize); grid on; box on;
xlabel('M/M_{sun}'); ylabel('T / 5772 K'); title('Effective temperature vs mass for MS stars');
print(gcf,'T vs M for MS.png','-dpng','-r300'); close(gcf);

%Construct HR diagram and overlay predicted gradients 
%for Main Sequence stars using homology arguments
p1 = plot( x, y,'k.', 'markersize',2 ); hold on;
p2 = plot( x, -4.1*x, 'r-'); p3 = plot( x, -5.6*x, 'r--');
p4 = plot( x, -8.6*x, 'g-'); p5 = plot( x, -9.5*x, 'b-');
set(gca,'fontsize',fsize); grid on; xlim(xlimits); ylim(ylimits); box on;
legend([p1,p2,p3,p4,p5],{'Star data','pp','pp constant \kappa','CNO','Triple \alpha'})
xlabel( '-log_{10}(T / 5772 K)'); ylabel('log_{10}(L/L_{sun})');
title(['Hertzsprung-Russell diagram for ',num2str(N),' near stars']);
print(gcf,'HR diagram with homology gradients.png','-dpng','-r300');

%Construct HR diagram, colour coded by temperature
clf; scatter( x, y,1,stars.RGB,'filled' ); hold on
plot( xxx, yyy ,'k-','linewidth',3);
plot( xx, yy ,'m*','linewidth',1);
set(gca,'fontsize',fsize); grid on;  box on;
xlim(xlimits); ylim(ylimits);
xlabel( '-log_{10}(T / 5772 K)'); ylabel('log_{10}(L/L_{sun})');
title(['Hertzsprung-Russell diagram for ',num2str(N),' near stars']);
print(gcf,'HR diagram.png','-dpng','-r300'); close(gcf);

%%

%Empirical Luminosity vs effective temperature variation from interpolation of
%Main Sequence in HR diagram. L is L/Lsun and T is T/sun
function T = T_from_L_MS(L)

%Define approximate yy = log10(L/Lsun),vs xx = -log10(T/Tsun) relationship for MS stars
xx = [-0.94,-0.38,-0.35,-0.32,-0.30,-0.27,-0.23,-0.20,-0.15,-0.12,-0.08,-0.05,-0.01,...
    0,0.01,0.04,0.07,0.10,0.13,0.16,0.17,0.18,0.19,0.21,0.22,0.25,0.26,0.3];
yy = [6.5,3.23,2.81,2.31,2.09,1.96,1.77,1.46,1.30,1.09,0.88,0.60,0.19,0,...
    -0.12,-0.37,-0.61,-0.84,-1.13,-1.37,-1.68,-2.13,-2.56,-2.98,-3.46,-4.03,-4.47,-6];

%Interpolate to find T
x = interp1( yy, xx, log10(L) );
T = 10.^(-x);

%%

%Mass-Luminosity ratio
% Empirical mass-luminosity ratio computed from binary systems, where mass
% and luminosity can be calculated directly from observations
% https://en.wikipedia.org/wiki/Mass%E2%80%93luminosity_relation
function L_by_Lsun = mass_luminosity(M)
L_by_Lsun = 0.23*M.^2.3;
L_by_Lsun( (M>0.43) & ( M<2 ) ) = M( (M>0.43) & ( M<2 ) ).^4;
L_by_Lsun( (M>2) & ( M<55 ) ) = 1.4*( M( (M>2) & ( M<55 ) ).^3.5);
L_by_Lsun( M>55 ) = 32000*M( M>55 );

%%

%Line of best fit function yfit = m*x, with product moment correlation
%coefficient r
function [yfit,xfit,r,m,dm,yupper,ylower,s] = bestfit(x,y)

%Find any x or y values that are NaN or Inf
ignore = isnan(abs(x)) | isnan(abs(y)) | isinf(abs(x)) | isinf(abs(y));
x(ignore) = []; y(ignore) = [];
n = length(x);

%Compute line of best fit
xbar = mean(x); ybar = mean(y); xybar = mean(x.*y);
xxbar = mean(x.^2 ); yybar = mean(y.^2 );
Vx = xxbar - xbar^2; Vy = yybar - ybar^2;
COVxy = xybar - xbar*ybar;
m = xybar/xxbar; r = COVxy/sqrt( Vx*Vy );
[x,i] = sort(x); y = y(i);
yfit = m*x; xfit = x;

%Compute error in gradient m
s = sqrt( (1/(n-1))*sum( (y - yfit).^2 ) );
dm = s/sqrt(n*Vx);

%Determine envelope lines
yupper = (m+dm)*x;
ylower = (m-dm)*x;

%End of code