Stellarcomputer/functions.h
/*
PC4181 MPhys Project - Stellar Computer
Functions for star.c
Started: 23/12/2005
Last Modified: 11/01/2006
Matthew Dennison (matthew.dennison@student.manchester.ac.uk) & Mike Peel (http://www.mikepeel.net/)
*/
/*
Include the standard functions
*/
#include <math.h>
#include <stdlib.h>
/*
Calculates the radius difference
See section 2.1 of the report for details.
*/
long double calc_radius_dif(long double prev_radius,
long double density)
{
long double output;
output = (1.0 / (4.0 * CONST_PI * pow(prev_radius, 2.0) * density));
return (output);
};
/*
Calculates the pressure difference using the equation of Hydrostatic Equilibrium
See section 2.2 of the report for details.
*/
long double calc_pressure_dif(long double mass,
long double radius)
{
long double output;
output = ((CONST_G * mass) / (4.0 * CONST_PI * pow(radius, 4.0)));
return (output);
};
/*
Calculates the luminosity
See section 2.3 of the report for details.
*/
long double calc_luminosity_dif(long double energy_generation)
{
long double output;
output = energy_generation;
return (output);
};
/*
Calculates the temperature
See section 2.4 of the report for details.
*/
long double calc_temperature_dif(long double temperature,
long double opacity,
long double luminosity,
long double radius,
long double mass,
long double pressure,
long double prev_temperature,
long double prev_pressure,
int i)
{
long double output;
// This calculates the start of the differential equation, i.e. before del.
output = (CONST_G * mass * temperature) / (4.0 * CONST_PI * pow(radius, 4.0) * pressure);
// Calculate the adiabatic gradient, del_ad
del_adiabatic[i] = ((log(temperature/prev_temperature)) / (log(pressure / prev_pressure)));
// Calculate the radiative gradient, del_rad
del_radiative[i] = (radiative_constant * opacity * luminosity * pressure) / (mass * pow(temperature, 4.0));
if (del_radiative[i] > del_adiabatic[i] && del_adiabatic[i] > 0.0)
{
// This is convective.
output = output * del_adiabatic[i];
}
else
{
// This is radiative
output = output * del_radiative[i];
};
return (output);
};
/*
Calculate mu, i.e. the average molecular mass
See section 2.5 of the report for details.
*/
long double calc_mu(long double X,
long double Y,
long double Z,
long double A_Z)
{
long double output;
output = MASS_H / ((X + (Y / 4.0) + (Z / A_Z)));
return (output);
};
/*
Calculates the density
Equation Of State - uses the Ideal Gas law plus the Radiation Pressure.
See section 2.5 of the report for details.
Note that mu contains m_h.
*/
long double calc_density(long double X,
long double temperature,
long double pressure)
{
long double output;
output = (mu / (CONST_K * temperature)) * (pressure + (1.0 / 3.0) * CONST_A * pow(temperature, 4.0));
return (output);
};
/*
Calculates the energy generation per unit volume at the specified temperature
See section 2.6 of the report for details.
*/
long double calc_energy_generation(long double temperature,
long double density,
long double X,
long double Y)
{
long double output, temp, temp_1_3, temp_2_3, alpha, psi, F_PPI, F_PPII, F_PPIII;
/*
For ease of calculation, work out temp = temperature / 10^6
Also, temp_1_3 and temp_2_3 will also be needed several times over...
*/
temp = temperature * 1e-6;
temp_1_3 = pow(temp, (1.0 / 3.0));
temp_2_3 = pow(temp, (2.0 / 3.0));
/*
F_ denotes fraction. The three variables determine the fractions of PPI, PPII and PPIII chains happening at the time. Used to adjust the output energy appropriately via psi.
Alpha is the number of alpha particles available; it impacts psi, F_PPI, F_PPII and F_PPIII.
This only becomes relevant when there's helium present in the star (which is nearly always the case
*/
if (Y > 0)
{
alpha = 1.2e17 * pow((Y / X), 2.0) * exp(-100 / temp_1_3);
F_PPI = (pow((1.0 + (2.0 / alpha)), 0.5) - 1.0) / (pow((1.0 + (2.0 / alpha)), 0.5) + 3.0);
F_PPII = (1 - F_PPI) / (1 + (6.3 / 7.05) * 1e16 * pow(temp, (-1.0 / 6.0)) * exp(-102.65 * pow(temp, -1.0 / 3.0)));
F_PPIII = (1 - F_PPI - F_PPII);
psi = (1.0 - alpha + alpha * pow((1.0 + (2.0 / alpha)), 0.5)) * (0.981 * F_PPI + 0.961 * F_PPII * 0.727 * F_PPIII);
}
else
{
psi = 1.0;
};
/*
This deals with the pp chain
*/
output = 2.36e2 * density * pow(X, 2.0) * pow(temp_2_3, -1) * exp(-33.81 / temp_1_3) * (1 + 0.0123 * temp_1_3 + 0.0109 * temp_2_3 * 0.00095 * temp) * psi;
/*
For debugging...
printf("%e %e %e %e %e %e", alpha, F_PPI, F_PPII, F_PPIII, psi, output);
*/
return (output);
};
/*
Calculates the opacity at the specified temperature
See section 2.7 of the report for details.
*/
long double calc_opacity(long double X,
long double density,
long double temperature)
{
long double output = 0.0;
long double opacity_1, opacity_2, opacity_3;
long double logR = 0.0;
long double logT = 0.0;
int k, k_use, l, l_use;
// Set up a switch between the three opacities.
int opacity_switch = 3;
switch(opacity_switch)
{
case 1:
/*
Use temperature boundaries
*/
if (temperature < 4e4)
{
output = 1.85e-14 * pow(density, 0.5) * pow(temperature, 4.0);
}
else if (temperature < 1.26e6)
{
output = 4e21 * (X + Y) * (1 + X) * density * pow(temperature, -3.5);
}
else
{
output = 0.02 * (1.0 + X);
};
break;
case 2:
/*
Use the middle one of the three opacity equations - this should in general be the appropriate one, and avoids the jump between the third and second equations when going outwards. Unfortunately, it seems to do this by ignoring the third equation completely in the sun...
*/
opacity_1 = 1.85e-14 * pow(density, 0.5) * pow(temperature, 4.0);
opacity_2 = 8e21 * density * pow(temperature, -3.5);
opacity_3 = 0.02 * (1.0 + X);
if ((opacity_1 < opacity_2) && (opacity_1 > opacity_3))
{
output = opacity_1;
}
else if ((opacity_2 < opacity_1) && (opacity_2 > opacity_3))
{
output = opacity_2;
}
else
{
output = opacity_3;
};
break;
case 3:
/*
Use OPAL opacity data from opacity.txt.
Note that we start at k = l = 1, not 0 - the 0th row/column is the value for logT and logR.
*/
logT = log10(temperature);
// 1000 changes to cgs units, which the OPAL data uses
logR = log10((density / 1000.0) * pow((1e-6 * temperature), 3.0));
/*
For debugging
printf("%Lf, %Lf\n", logT, logR);
*/
// Find the closest T
for(k = 1; k < NUM_OPAL_I; k++)
{
// printf("%d\n", k);
if ((opal_data[k][0] - logT) >= 0.0 || k == (NUM_OPAL_I - 1))
{
// Check and see if the previous value was closer or not
k_use = k;
if (k != 1)
{
if ((abs(opal_data[k-1][0] - logT)) < abs(opal_data[k][0] - logT))
{
k_use = k - 1;
};
};
// Find the closest R
for(l = 1; l < NUM_OPAL_J; l++)
{
if ((opal_data[0][l] - logR) >= 0.0 || l == (NUM_OPAL_J - 1))
{
// Check and see if the previous value was closer or not
l_use = l;
if (l != 1)
{
if ((abs(opal_data[0][j-1] - logR)) < abs(opal_data[0][l] - logR))
{
l_use = l - 1;
};
};
output = opal_data[k_use][l_use];
break;
};
};
};
};
// 10 changes back to SI units.
output = pow(10, output) * 10.0;
break;
};
return (output);
};
/*
The Runge-Kutta function. Calculates the next set of differences depending on the previous set.
All functions other than dir, i and j are globally accessible, so aren't explicitely passed.
dT[j], dP[j] etc. are defined previous to this function being called.
dir determines the direction we're going in - it's either -1 or 1.
*/
void rungekutta(int dir, int i, int j)
{
temp_density[j] = calc_density(X, temperature[i - dir] + dT[j], pressure[i - dir] + dP[j]);
temp_radius[j] = dm[i] * calc_radius_dif(radius[i - dir] + dR[j], *ptr_density);
temp_pressure[j] = dm[i] * calc_pressure_dif(mass[i - dir] + dM[j], radius[i - dir] + dR[j]);
temp_opacity[j] = calc_opacity(X, *ptr_density, temperature[i - dir] + dT[j]);
temp_temperature[j] = dm[i] * calc_temperature_dif(temperature[i - dir] + dT[j], *ptr_opacity, luminosity[i - dir] + dL[j], radius[i - dir] + dR[j], mass[i - dir] + dM[j], pressure[i - dir] + dP[j], old_temperature, old_pressure, i);
temp_energy[j] = calc_energy_generation(temperature[i - dir] + dT[j], *ptr_density, X, Y);
temp_luminosity[j] = dm[i] * calc_luminosity_dif(*ptr_energy);
};
