%MS_stars
% Generates a set of graphs which models the physical parameters of Main
% Sequence stars, using a polytropic model of gas pressure with density.
% Stellar interiors are modelled as ideal gases, and radiation pressure
% and degeneracy effects are ignored.
%
% LAST UPDATED by Andy French Jan 2026

function MS_stars

%% INPUTS %%

%Known solar parameters
Msun = 1.989e30;  %Solar mass /kg
Rsun = 696340e3;  %Solar radius /m
Lsun = 3.846e26;  %Solar luminosity /W
Tsun = 5772;      %Solar surface temperature /K
rho_sun = Msun/( (4/3)*pi*Rsun^3 );  %Solar average density (kg/m^3)
X_H_sun = 0.747; %Hydrogen mass fraction
X_He_sun = 0.236; %Helium mass fraction


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Define a range of masses  which span the Main Sequence. Note log scale.
% MS stars are between 0.08 and about 120 solar masses.
% Based upon the HR diagram, the hottest 'common' MS stars are about 14,500K
% whereas the coolest are about 3250K. But some blue stars can have
% surface temperatures possibly up to 50,000K. But they probably don't last
% very long, as the fusion power per kg rises very steeply with temperature
% for fusion of nuclei which are heavier than hydrogen.
N = 10; M = 10.^(linspace( log10(0.1), log10(90), N ));

%Set hydrogen and helium mass fractions to be the same as the Sun.
X_H = X_H_sun*ones(1,N); X_He = X_He_sun*ones(1,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

%Set fontsize for graphs
fsize = 18;
lgnd_fsize = 0.8*fsize;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Sun input
% N=1; M = 1; X_H = X_H_sun;  X_He = X_He_sun; RGB = [1,0,0];
% lgnd_fsize = fsize;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Make folder for graphs
dirname = ['stars N=',num2str(N)];
mkdir(dirname);

%

%% COMPUTE MODEL %%

%Compute luminosity (W) and surface temperature (K) using empirical relationships
L_by_Lsun = mass_luminosity(M);
Ts = Tsun*T_from_L_MS(L_by_Lsun);

%Run Main Sequence star model for each of the M,Ts inputs
r = []; r_by_R = []; R = []; m = []; M_by_Msun = []; rho = []; rho0 = []; rho_by_rho0 = [];
l = []; k = []; k_calc = [];  P = []; Prad = [];
T = []; T0 = []; dT_by_dr = [];  dT_by_dr_conv = []; dT_by_dr_rad = []; conv_instability = [];
lgnd = {}; fusion_model = {};
for n=1:N
    
    %Run model. Note M in solar masses
    s = MS_star_model( M(n),Ts(n),X_H(n),X_He(n) );
    
    %Add to arrays to enable efficient plotting
    M_by_Msun = [M_by_Msun,M(n)*ones(length(s.r),1)];
    r = [r,s.r.']; R = [R,s.R];
    r_by_R = [r_by_R,s.r'/s.R];
    m = [m,s.m.'];
    rho = [rho,s.rho.'];
    rho_by_rho0 = [rho_by_rho0,s.rho.'/s.rho0];
    l = [l,s.l.'];
    k = [k,s.k.'];
    k_calc = [k_calc,s.k_calc.'];
    T = [T,s.T.'];
    dT_by_dr = [dT_by_dr ,s.dT_by_dr .'];
    dT_by_dr_conv = [dT_by_dr_conv ,s.dT_by_dr_conv .'];
    dT_by_dr_rad = [dT_by_dr_rad ,s.dT_by_dr_rad .'];
    conv_instability = [ conv_instability, ( abs(s.dT_by_dr_rad) > abs(s.dT_by_dr_conv) ).' ];
    P = [P,s.P.'];
    Prad = [Prad,s.Prad.'];
    T0 = [T0,s.T0];
    rho0 = [rho0,s.rho0];
    fusion_model{n} = s.fusion_model;
    
    %Build up legend for graphs
    lgnd{n} = ['M/Msun = ',num2str(M(n),2)] ;
end

%

%% PLOT GRAPHS %%

%Plot surface of convective instability (1 means yes) against m/M vs M/Msun
pcolor( M_by_Msun, m./(Msun*M_by_Msun), conv_instability ); icmap; shading interp; colorbar;
view(2); axis tight
set(gca,'fontsize',fsize); xlabel('M/M_{sun}'); ylabel('m/M'); grid on; box on;
title('Convective (1) vs Radiative (0) stellar interiors');
print(gcf,[dirname,'\convective instability.png'],'-r300','-dpng'); close(gcf);

%Generate mass vs range plot
hold on; for n=1:N; plot(r_by_R(:,n),m(:,n)/Msun,'linewidth',2,'color',RGB(n,:)); end
set(gca,'fontsize',fsize);
xlabel('r/R'); ylabel('mass enclosed / M_{sun}');
if lgnd_fsize>0; legend( lgnd,'fontsize',lgnd_fsize ); end
grid on; title('Enclosed mass vs radius fraction'); box on;
print(gcf,[dirname,'\m vs r_by_R.png'],'-r300','-dpng'); clf;

%Generate fractional density vs fractional range plot
hold on; for n=1:N; plot(r_by_R(:,n),rho_by_rho0(:,n),'linewidth',2,'color',RGB(n,:)); end;
set(gca,'fontsize',fsize);
xlabel('r/R'); ylabel('rho / \rho_{0}');
if lgnd_fsize>0; legend( lgnd,'fontsize',lgnd_fsize ); end
grid on; title('\rho / \rho_{0} vs fractional radius'); box on;
print(gcf,[dirname,'\rho vs r_by_R.png'],'-r300','-dpng'); close(gcf);

%Plot core density vs M
plot(M,rho0,'r.-','linewidth',2); hold on;
set(gca,'fontsize',fsize);
xlabel('M/M_{sun}'); ylabel('Core density \rho_{0} (kgm^{-3})');
grid on; title('Core density \rho_{0} vs M/M_{sun}'); box on;
print(gcf,[dirname,'\rho0 vs M.png'],'-r300','-dpng'); close(gcf);

%Plot star radius vs M
plot(M,R/Rsun,'r.-','linewidth',2); hold on;
set(gca,'fontsize',fsize);
xlabel('M/M_{sun}'); ylabel('R/R_{sun}');
grid on; title('Star radius R/R_{sun} vs star mass M/M_{sun}'); box on;
print(gcf,[dirname,'\R vs M.png'],'-r300','-dpng'); close(gcf);

%Generate log10(luminosity) vs fractional range plot
hold on; for n=1:N; plot(r_by_R(:,n),log10(l(:,n)/Lsun),'linewidth',2,'color',RGB(n,:)); end;
set(gca,'fontsize',fsize);  ylim([-5,5]);
xlabel('r/R'); ylabel('log_{10}( L/L_{sun} )');
if lgnd_fsize>0; legend( lgnd,'fontsize',lgnd_fsize,'location','southeast' ); end
grid on; title('Luminosity vs r/R'); box on;
print(gcf,[dirname,'\l vs r_by_R.png'],'-r300','-dpng'); close(gcf);

%Generate calculated opacity vs fractional range plot
hold on; for n=1:N; plot(r_by_R(:,n),k_calc(:,n),'linewidth',2,'color',RGB(n,:)); end;
set(gca,'fontsize',fsize); 
if lgnd_fsize>0; legend( lgnd,'fontsize',lgnd_fsize,'location','northwest' ); end
ylimits = get(gca,'ylim');
xlabel('r/R'); ylabel('Opacity (m^2kg^{-1})');
grid on; title('Calculated opacity vs fractional radius'); box on;
print(gcf,[dirname,'\opacity vs r_by_R.png'],'-r300','-dpng'); close(gcf);

%Generate Kramer's opacity vs fractional range plot
hold on; for n=1:N; plot(r_by_R(:,n),k(:,n),'linewidth',2,'color',RGB(n,:)); end;
set(gca,'fontsize',fsize); 
if lgnd_fsize>0; legend( lgnd,'fontsize',lgnd_fsize,'location','northwest' ); end
set(gca,'ylim',ylimits);
xlabel('r/R'); ylabel('Opacity (m^2kg^{-1})');
grid on; title('Kramers opacity vs fractional radius'); box on;
print(gcf,[dirname,'\kramers opacity vs r_by_R.png'],'-r300','-dpng'); close(gcf);

%Generate temperature vs range plot
hold on; for n=1:N; plot(r(:,n)/Rsun,T(:,n)/Tsun,'linewidth',2,'color',RGB(n,:)); end;
set(gca,'fontsize',fsize);
xlabel('r/R_{sun}'); ylabel('T / 5772 K');
if lgnd_fsize>0; legend( lgnd,'fontsize',lgnd_fsize ); end
grid on; title('Temperature vs radius'); box on;
print(gcf,[dirname,'\T vs r.png'],'-r300','-dpng'); close(gcf);

%Generate temperature vs fractional range plot
hold on; for n=1:N; plot(r_by_R(:,n),T(:,n)/Tsun,'linewidth',2,'color',RGB(n,:)); end;
set(gca,'fontsize',fsize);
xlabel('r/R'); ylabel('T / 5772 K');
if lgnd_fsize>0; legend( lgnd,'fontsize',lgnd_fsize ); end
grid on; title('Temperature vs fractional radius'); box on;
print(gcf,[dirname,'\T vs r_by_R.png'],'-r300','-dpng'); close(gcf);

%Plot core temperature vs M
plot(M,T0/Tsun,'r.-','linewidth',2); hold on;
xlimits = get(gca,'xlim');
plot(xlimits,(1.5e7/Tsun)*[1,1],'b-','linewidth',2);
plot(xlimits,(1e8/Tsun)*[1,1],'m-','linewidth',2);
legend({'T_{0}/5772 K','CNO','Triple \alpha'},'fontsize',fsize );
set(gca,'fontsize',fsize);
xlabel('M/M_{sun}'); ylabel('T / 5772 K');
grid on; title('Core temperature T_0 vs M/M_{sun}'); box on;
print(gcf,[dirname,'\T0 vs M.png'],'-r300','-dpng'); close(gcf);

%Plot temperature gradient vs fractional range
hold on; for n=1:N; plot(r_by_R(:,n),log10(dT_by_dr(:,n)),'linewidth',2,'color',RGB(n,:)); end;
set(gca,'fontsize',fsize); ylim([-5,0]);
xlabel('r/R'); ylabel('log_{10}( dT/dr (K/m) )');
if lgnd_fsize>0; legend( lgnd,'fontsize',lgnd_fsize ); end
grid on; title('log_{10}( dT/dr ) vs r/R'); box on;
print(gcf,[dirname,'\dT_by_dr vs r_by_R.png'],'-r300','-dpng'); close(gcf);

%Plot |temperature gradient|  - |convective dT/dr| vs fractional range. If positive this
%implies convective instablity.
hold on;
for n=1:N
    plot(r_by_R(:,n), abs( dT_by_dr_rad(:,n) ) - abs( dT_by_dr_conv(:,n) ),...
        'linewidth',2,'color',RGB(n,:));
end
set(gca,'fontsize',fsize);
xlabel('r/R'); ylabel('dT/dr difference (K/m)');
if lgnd_fsize>0; legend( lgnd,'fontsize',lgnd_fsize ); end
grid on; title('Temperature gradient difference from adiabatic'); box on;
print(gcf,[dirname,'\dT_by_dr minus dT_by_dr_conv vs r_by_R.png'],'-r300','-dpng'); close(gcf);

%Generate ratio of radiation pressure to ideal gas pressure vs fractional range plot
hold on; 
for n=1:N; 
    plot(r_by_R(:,n),log10(Prad(:,n)./P(:,n)),'linewidth',2,'color',RGB(n,:)); 
end;
set(gca,'fontsize',fsize); xlim([0,1]); axes_label_N_sf(gca,3);
if lgnd_fsize>0; legend( lgnd,'fontsize',lgnd_fsize ); end
xlabel('r/R'); ylabel('log_{10}( P_{rad}/P_{ideal gas} )');
grid on; title('Radiation pressure / ideal gas pressure vs radius'); box on;
print(gcf,[dirname,'\Prad vs P vs r_by_R.png'],'-r300','-dpng'); close(gcf);

%%

%MS_star_model
% Computes enclosed mass, density, temperature, temperature gradient,
% pressure, pressure gradient, and luminosity vs radius of Main Sequence
% stars. (i.e. model isn't appropriate for Red Giants or White Dwarfs).
% Assumes ideal gases and ignores radiation (and degeneracy)
% pressure. Uses numeric solutions to the Lane-Emden differential equation,
% which occurs in models of quasi-static stellar structure for a polytropic
% star of pressure vs density profile P = k * rho^(1+1/n) where n is the
% polytropic index. For MS stars, n=3.
%
% Inputs are:
% star mass M (in solar masses)
% star surface temperature Ts (in K)
% hydrogen (X_H) and helium (X_He) mass fractions. ('Metalicity' Z = 1 - X_H - X_He).
function s = MS_star_model(M,Ts,X_H,X_He)

%Universal physical constants
kB = 1.38e-23;   %Boltzmann constant (J/K)
G = 6.67e-11;  %Universal Gravitational constant (m^3 kg^-1 s^-2)
mp = 1.67e-27; %Proton mass (kg)
sigma = 5.67e-8;  %Setfan-Boltzmann constant (Wm^-2 K^-4)
c = 2.998e8; %Speed of light (m/s)
gmma = 5/3; %Ratio of specific heats for fully ionized gas
y = 1; z = 3.5;  %Density and temperature powers for Kramer's opacity model

%Set polytropic index
n_poly = 3;

%Known solar parameters
Msun = 1.989e30;  %Solar mass /kg
Rsun = 696340e3;  %Solar radius /m
Lsun = 3.846e26;  %Solar luminosity /W
Tsun = 5772;      %Solar surface temperature /K
rho_sun = Msun/( (4/3)*pi*Rsun^3 );  %Solar average density (kg/m^3)
X_H_sun = 0.747; %Hydrogen mass fraction
Y_H_sun = 0.236; %Helium mass fraction

%

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

%Determine star luminosity (W) and mass M in kg
s.L = Lsun*mass_luminosity(M); s.M = M*Msun;

%Determine star radius (m)
s.R = sqrt( s.L/(4*pi*sigma*Ts^4) );

%Determine average density (kg/m^3)
s.rho_mean = s.M/( (4/3)*pi*s.R^3 );

%Determine surface outputs of Lane-Emden equation. Works up to n=4.
dY = 0.0001; Ymax = 20; [Q,dQ_by_dY,Y0,k,Y] = lane_emden_solver(n_poly,dY,Ymax);

%Determine star core density using n=3 polytrope solution to Lane-Emden
%equation
s.condensation_factor = ( (Y0^3) / (3*k) ); 
s.rho0 = s.condensation_factor*s.rho_mean;

%Remove Lane-Emden outputs where Q<0  (i.e. surface of star has been reached).
i = find(Q<0); Q(i) = []; dQ_by_dY(i) = []; Y(i) = [];

%Assuming ideal gas model, determine core temperature (K)
s.T0 = (4*pi*G*s.mu*s.rho0*s.R^2)/( kB*(n_poly+1)*Y0^2 );

%Determine fusion model powers, given core temperature T0
[AX,AY,AZ,B,C,s.fusion_model] = fusion_model_powers(s.T0);

%Determine Lane Emden radial parameter (in m)
s.alpha = sqrt( kB*s.T0*(n_poly+1)/(4*pi*G*s.mu*s.rho0 ) );

%Star radial extent (m) and star radius. s.R_check should equal s.R.
s.r = s.alpha*Y; s.R_check = s.alpha*Y0;

%Star enclosed mass (kg), and total star mass (kg)
s.m = 4*pi*(s.alpha^3)*s.rho0*( - (Y.^2) .* dQ_by_dY );
s.M = 4*pi*(s.alpha^3)*s.rho0*k;
s.M_check = s.M/(M*Msun); 

%Star density (kg/m^3) vs r and average density
s.rho = s.rho0 * Q.^n_poly ;  s.rho_bar = s.rho0*3*k/(Y0^3);

%Determine pressure (Pa) from density using polytropic model and ideal gas
%equation
s.P0 = s.rho0*kB*s.T0/s.mu; s.P = s.P0*(s.rho/s.rho0).^(1+1/n_poly);

%Assuming ideal gas equation, determine temperature (K)
s.T = s.P*s.mu./(s.rho*kB);

%Check surface temperature (K) with input
s.Ts = s.T(end);

%Determine temperature gradient dT/dr
s.dT_by_dr = diff(s.T)./diff(s.r); s.dT_by_dr = [s.dT_by_dr(1),s.dT_by_dr];

%Calculate radiation pressure (Pa)
s.Prad = (4*sigma/(3*c))*s.T.^4;

%Determine pressure gradient dP/dr
s.dP_by_dr = diff(s.P)./diff(s.r); s.dP_by_dr = [s.dP_by_dr(1),s.dP_by_dr];

%Calculate temperature gradient assuming convection.
%Schwarzschild criterion for convective instability is |dT/dr| is greater
%than |dT_by_dr_conv|
s.dT_by_dr_conv = (1- 1/gmma)*s.T.*s.dP_by_dr./s.P;

%Integrate luminosity generated by nuclear fusion within layers of mass dm
%assuming density and temperature parametrized power law for energy
%generated per unit mass
N = length(s.r); l = zeros(1,N);
for n=1:(N-1)
    
    %Shell width (in m)
    dr = s.r(n+1)-s.r(n);
    
    %Mass of shell of width dr (in kg)
    dm = 4*pi*( s.r(n)^2 )*s.rho(n)*dr;
    
    %Contribution to luminosity (in W) using nuclear fusion power model
    yes = fusion_yes_or_no( s.T(n), s.fusion_model);
    if yes==1;
        %Temperature is hot enough for nuclear fusion
        dl = dm * ( X_H^AX )*( X_He^AY )*( (1-X_H-X_He)^AZ )*( s.rho(n)^B )*( s.T(n)^C );
    else
        %Assume temperature is not hot enough for fusion, so therefore no 
        %further contributions of stellar luminosity
        dl = 0;
    end
    
    %Cumulatively sum luminosity
    l(n+1) = l(n) + dl;
end

%Scale such that luminosity at surface is L
s.l = s.L*l/l(N);

%Calculate opacity. Assume matching of constant opacity (low energy
%Compton scattering) at core
s.k0 = 0.02*(1+X_H)/( (s.rho0^y)/(s.T0^z) );
s.k = s.k0*( s.rho.^y ).*(s.T.^-z);
s.k_check = 0.02*(1+X_H)*( s.rho/s.rho0 ).^( 1 - z/n_poly );

%Determine dT/dr assuming radiative transfer only
s.dT_by_dr_rad = -(3/4)*s.k.*s.rho.*s.l./( 16*pi*(s.r.^2) * sigma .*(s.T.^ 3) );

%Calculate opacity assuming radiative transfer and previously computed
%luminosity
s.k_calc = -64*pi*(s.r.^2).*sigma.*(s.T.^3).*s.dT_by_dr./( 3*s.rho.*s.l );

%%

%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 );

%%

%Fusion model. Rate of energy generation /kg is proportional to
% X^AX, Y^AY, Z^AZ, rho^B, T^C. Assumes a single dominant process.
function [AX,AY,AZ,B,C,fusion_model] = fusion_model_powers(T0)
%Default is pp chain
fusion_model = 'pp';
AX = 2; AY = 0; AZ=0; B = 1; C = 4;
if ( T0>1.5e7 ) & (T0<1e8)
    %CNO cycle
    fusion_model = 'CNO';
    AX = 1; AY = 0; AZ=1; B = 1; C = 17;
elseif T0>1e8
    %Triple-alpha
    fusion_model = 'Triple-alpha';
    AX = 0; AY = 3; AZ=0; B = 2; C = 40;
end

%Ignores carbon, oxygen, silicon etc fusion. One assumes this is
%not possible for main sequence 'quasi-static' stars, since timescale is 300 years or
%less.

%%

%Fusion yes on no? Depending on fusion model, fusion allowed if T is higher
%than a particular threshold. T in K. yes = 1 means fusion occurs.
% Assume a single fusion model for a star. In practice, one might imagine
% CNO or triple alpha near the core, and pp at larger radii?
function yes = fusion_yes_or_no(T, fusion_model)
if strcmp( fusion_model,'pp')
    if T>=4e6
        yes=1;
    else
        yes=0;
    end
elseif strcmp( fusion_model,'CNO')
    if T>=1.5e7
        yes=1;
    else
        yes=0;
    end
elseif strcmp( fusion_model,'Triple-alpha')
    if T>=1e8
        yes=1;
    else
        yes=0;
    end
else
    yes = 0;
end

%%

%Euler method solver for Lane-Emden equation
%k is -Y^2 * dQ_by_dY evaluated ay Y0, i.e. when Q=0
function [Q,dQ_by_dY,Y0,k,Y] = lane_emden_solver(n,dY,Ymax)

%Determine Y array
Y = 0:dY:Ymax;

%Analytic cases (n=2 has a series expansion, and n=5 is solvable when Y>0,
%but I have ignored these!
%https://en.wikipedia.org/wiki/Lane%E2%80%93Emden_equation)
if n==0
    Q = 1 - (1/6)*Y.^2; dQ_by_dY = -(1/3)*Y; Y0 = sqrt(6); k = 2*sqrt(6);
elseif n==1
    Q = sin(Y)./Y; dQ_by_dY = ( Y.*cos(Y) - sin(Y) )./(Y.^2); Y0 = pi; k = pi;
else
    %Use default for n=5. This will be reset if i0 is not NaN i.e. Q<0
    E = zeros(size(Y)); Q = ones(size(Y)); Y0calc=1; i0 = NaN;
    Y0 = inf; k = 1.7;
    
    %Apply Euler method to solve
    for i = 1:(length(Y)-1)
        E(i+1) = E(i) + ( Y(i)^2 ) * ( Q(i)^n ) * dY;
        if Y(i)==0;
            Q(i+1) = 1;
        else
            Q(i+1) = Q(i) - E(i)*dY/( Y(i)^2 );
        end
        if ( Q(i+1)<0 ) && (Y0calc==1)
            i0 = i; Y0calc = 0;  %This corresponds to the star radius (if it has one!)
        end
    end
    dQ_by_dY = -E./( Y.^2 );
    if ~isnan(i0)
        Y0 = Y(i0); k = - Y0^2 * dQ_by_dY(i0);
    end
end

%%

%Function to force axes labels to be to N significant figures
function axes_label_N_sf(ax,N)
%Get x tick labels
a = get(ax);
xticklabels = a.XTickLabel;
for n=1:length(xticklabels);
    xticklabels{n} = num2str( str2num( xticklabels{n} ), N );
end
yticklabels = a.YTickLabel;
for n=1:length(yticklabels);
    yticklabels{n} = num2str( str2num( yticklabels{n} ), N );
end
set(ax,'XTickLabel',xticklabels,'YTickLabel',yticklabels);

%%

%Interpolates current colormap to yield better graduated shading.
function icmap
N = 1000; %Legth of new colormap
map = colormap; new_map = ones(N,3); %Get current colormap

%Get size of current colormap and initalise R,G,B vectors
dim = size(map);
R = ones(1,dim(1)); G = ones(1,dim(1)); B = ones(1,dim(1));
RR = ones(1,N); GG = ones(1,N); BB = ones(1,N);

%Populate R,G,B with current colormap
R(:) = map(:,1); G(:) = map(:,2); B(:) = map(:,3);

%Interpolate R,G,B to yield new colour map
x = linspace( 1, dim(1), N );
RR = interp1( 1:dim(1), R, x ); GG = interp1( 1:dim(1), G, x ); BB = interp1( 1:dim(1), B, x );
new_map(:,1) = RR(:); new_map(:,2) = GG(:); new_map(:,3) = BB(:);

%Set colormap to be new map
colormap( new_map );

%End of code