Stellarcomputer/star.c
/*
PC4181 MPhys Project - Stellar Computer
Calculating the properties of a star given a set of input numbers
Started: 11/10/2005
Last Modified: 09/01/2006
Matthew Dennison (matthew.dennison@student.manchester.ac.uk) & Mike Peel (http://www.mikepeel.net/)
*/
/*
Include the standard functions
*/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
/*
Set up the constants
N is the number of times we're doing the calculation
NUM_OPAL_I and NUM_OPAL_J are the number of rows and columns respectively (inc. logT and logR) of the opacity.txt file.
*/
#define N 100000
#define CONST_PI 3.14159
#define CONST_G 6.674e-11
#define CONST_K 1.381e-23
#define CONST_A 7.5657e-16
#define CONST_H 6.626e-34
#define CONST_C 299792458.0
#define MASS_H 1.673e-27
#define MASS_e 9.109e-31
#define NUM_OPAL_I 71
#define NUM_OPAL_J 20
/*
Set the input/output directory
*/
char directory[500] = "/star/Current/";
/*
Set up the arrays
*/
long double radius[N], err_radius[N];
long double mass[N];
long double temperature[N], err_temperature[N];
long double pressure[N], err_pressure[N];
long double density[N];
long double luminosity[N], err_luminosity[N];
long double opacity[N];
long double energy[N];
long double dm[N], min_dm, new_dm;
long double degenerate_pressure[N], fermi_pressure;
long double del_adiabatic[N], del_radiative[N];
long double mu;
long double X, Y, Z, A_Z;
long double opal_data[71][20];
/*
Set up the variables used throughout the application
The temp_ variables are used in the Runge-Kutta routines.
The set number variables are the coefficients for the various elements within the RUk w/ 5th order error checking
*/
long double temp_density[6],
temp_radius[6],
temp_pressure[6],
temp_opacity[6],
temp_temperature[6],
temp_energy[6],
temp_luminosity[6];
long double *ptr_density, *ptr_opacity, *ptr_energy;
long double a2 = 0.2,
a3 = 3.0/10.0,
a4 = 3.0/5.0,
a5 = 1.0,
a6 = 7.0/8.0;
long double b21 = 0.2,
b31 = 3.0/40.0,
b32 = 9.0/40.0,
b41 = 3.0/10.0,
b42 = -9.0/10.0,
b43 = 6.0/5.0,
b51 = -11.0/54.0,
b52 = 2.5,
b53 = -70.0/27.0,
b54 = 35.0/27.0,
b61 = 1631.0/55296.0,
b62 = 175.0/512.0,
b63 = 575.0/13824.0,
b64 = 44275.0/110592.0,
b65 = 253.0/4096.0;
long double c1 = 37.0/378.0,
c2 = 0.0,
c3 = 250.0/621.0,
c4 = 125.0/594.0,
c5 = 0.0,
c6 = 512.0/1771.0;
long double C1 = 2825.0/27648.0,
C2 = 0.0,
C3 = 18575.0/48384.0,
C4 = 13525.0/55296.0,
C5 = 277.0/14336.0,
C6 = 0.25;
long double dP[6], dT[6], dR[6], dL[6], dM[6];
long double P1, T1, R1, L1;
// error is used to check the numerical error and change dm
long double error = 1e-8;
long double old_temperature, old_pressure;
int direction, loop, i, i_min, i_max, k, j;
// accuracy determines the number of decimal places the program outputs.
int accuracy = 12;
FILE *input_file, *output_file, *opacity_file;
char filename[100];
/*
Work out the radiative constant for later
*/
long double radiative_constant = (3.0 / (16.0 * CONST_PI * CONST_A * CONST_C * CONST_G));
/*
Get the functions
*/
#include "functions.h"
/*
This is the main program
*/
int main()
{
/*
Get the input data from input.txt
*/
sprintf(filename, "%sinput.txt", directory);
input_file = fopen(filename, "r");
if (!input_file)
{
printf("Failed to open %sinput.txt \n", directory);
return 1;
};
fscanf(input_file, "Radius: %Lf\n", &radius[0]);
fscanf(input_file, "Mass: %Lf\n", &mass[0]);
fscanf(input_file, "Surface Temperature: %Lf\n", &temperature[0]);
fscanf(input_file, "Luminosity: %Lf\n", &luminosity[0]);
fscanf(input_file, "Core Temperature: %Lf\n", &temperature[N-1]);
fscanf(input_file, "Core Pressure: %Lf\n", &pressure[N-1]);
fscanf(input_file, "X: %Lf\n", &X);
fscanf(input_file, "Y: %Lf\n", &Y);
fscanf(input_file, "Z: %Lf\n", &Z);
fscanf(input_file, "A_Z: %Lf\n", &A_Z);
fclose(input_file);
/*
Get the OPAL data
*/
sprintf(filename, "%sopacity.txt", directory);
opacity_file = fopen(filename, "r");
if (!opacity_file)
{
printf("Failed to open %sopacity.txt \n", directory);
return 1;
};
// For each column in the csv file...
for (i = 0; i < NUM_OPAL_I; i++)
{
// For each row in the csv file...
for (j = 0; j < NUM_OPAL_J; j++)
{
fscanf(opacity_file, "%Lf", &opal_data[i][j]);
/*
For debugging...
printf("%Lf", opal_data[i][j]);
*/
};
/*
For debugging...
printf("\n");
*/
};
fclose(opacity_file);
/*
Set a minimum mass to be 100x smaller than that defined by N - that's the minimum we should need to use for accuracy of around 1e-8.
*/
min_dm = mass[0] / (100 * N);
// Calculate mu
mu = calc_mu(X, Y, Z, A_Z);
/*
Set the remaining boundary conditions at the surface and the centre.
Pressure uses the surface pressure equation. See section 3.2 of the report for details. 1.85e-14 is kappa_3.
*/
dm[1] = min_dm;
pressure[0] = pow(((CONST_G * mass[0] * pow(CONST_K, 0.5)) / (pow(radius[0], 2.0) * 1.85e-14 * pow(temperature[0], 3.5) * pow(mu, 0.5))), 0.66);
density[0] = calc_density(X, temperature[0], pressure[0]);
opacity[0] = 1.85e-14 * pow(density[0], 0.5) * pow(temperature[0], 4.0);
/*
Get the values for the core.
Temperature and pressure have been defined when input.txt is read.
*/
dm[N-1] = min_dm;
mass[N-1] = 0.0;
luminosity[N-1] = 0.0;
radius[N-1] = 0.0;
density[N-1] = calc_density(X, temperature[N-1], pressure[N-1]);
opacity[N-1] = calc_opacity(X, density[N-1], temperature[N-1]);
energy[N-1] = calc_energy_generation(temperature[N-1], density[N-1], X, Y);
/*
Get the values for the second innermost shell.
See section 3.1 of the report for details.
*/
dm[N-2] = min_dm;
mass[N-2] = dm[N-1];
radius[N-2] = pow(((3.0 * mass[N-2]) / (4.0 * CONST_PI * density[N-1])), 1.0 / 3.0);
luminosity[N-2] = energy[N-1] * mass[N-2];
pressure[N-2] = pressure[N-1] - 0.5 * pow((4.0 * CONST_PI / 3.0), (1.0 / 3.0)) * CONST_G * pow(density[N-1], (4.0 / 3.0)) * pow(mass[N-2], (2.0 / 3.0));
del_radiative[N-2] = radiative_constant * opacity[N-1] * energy[N-1] * pressure[N-1] * pow(temperature[N-1], -4.0);
temperature[N-2] = temperature[N-1] - (3.0 * opacity[N-1] * luminosity[N-2] * mass[N-2])/(256.0 * CONST_PI * CONST_PI * 5.67e-8 * pow(radius[N-2],4.0) * pow(temperature[N-1], 3.0));
density[N-2] = calc_density(X, temperature[N-2], pressure[N-2]);
opacity[N-2] = calc_opacity(X, density[N-2], temperature[N-2]);
energy[N-2] = calc_energy_generation(temperature[N-2], density[N-2], X, Y);
dm[N-3] = min_dm;
for (loop = 0; loop < 2; loop++)
{
if (loop == 1)
{
/*
This is the loop starting at the core, and working outwards.
Direction is negative as i is decreasing - i=N is the center.
*/
direction = -1;
i_min = N - 3;
i_max = 0.5 * N;
}
else
{
/*
This is the loop starting on the surface, and working inwards
Direction is positive as i is increasing - i=1 is the surface.
*/
direction = +1;
i_min = 1;
i_max = 0.5 * N;
};
for (i = i_min; (loop == 0) ? i < i_max : i > i_max; i += direction)
{
/*
Mass
*/
mass[i] = mass[i - direction] - direction * dm[i];
/*
First Runge-Kutta loop
*/
dP[0] = 0.0;
dT[0] = 0.0;
dR[0] = 0.0;
dL[0] = 0.0;
dM[0] = 0.0;
ptr_density = &density[i - direction];
ptr_energy = &energy[i - direction];
ptr_opacity = &opacity[i - direction];
old_temperature = (temperature[i - direction * 2] + temperature[i - direction]) * b21;
old_pressure = (pressure[i - direction * 2] + pressure[i - direction]) * b21;
rungekutta(direction, i, 0);
/*
Second Runge-Kutta loop
*/
dP[1] = temp_pressure[0] * b21;
dT[1] = temp_temperature[0] * b21;
dR[1] = temp_radius[0] * b21;
dL[1] = temp_luminosity[0] * b21;
dM[1] = dm[i] * a2;
ptr_density = &temp_density[0];
ptr_energy = &temp_energy[0];
ptr_opacity = &temp_opacity[0];
old_temperature = ((temperature[i - direction * 2] + temperature[i - direction]) * b21);
old_pressure = ((pressure[i - direction * 2] + pressure[i - direction]) * b21);
rungekutta(direction, i, 1);
/*
Third Runge-Kutta loop
*/
dP[2] = (temp_pressure[0] * b31) + (temp_pressure[1] * b32);
dT[2] = (temp_temperature[0] * b31) + (temp_temperature[1] * b32);
dR[2] = (temp_radius[0] * b31) + (temp_radius[1] * b32);
dL[2] = (temp_luminosity[0] * b31) + (temp_luminosity[1] * b32);
dM[2] = dm[i] * a3;
ptr_density = &temp_density[1];
ptr_energy = &temp_energy[1];
ptr_opacity = &temp_opacity[1];
old_temperature = (temperature[i - direction * 2] + temperature[i - direction]) / 2;
old_pressure = (pressure[i - direction * 2] + pressure[i - direction]) / 2;
rungekutta(direction, i, 2);
/*
Fourth Runge-Kutta loop
*/
dP[3] = (temp_pressure[0] * b41) + (temp_pressure[1] * b42) + (temp_pressure[2] * b43);
dT[3] = (temp_temperature[0] * b41) + (temp_temperature[1] * b42) + (temp_temperature[2] * b43);
dR[3] = (temp_radius[0] * b41) + (temp_radius[1] * b42) + (temp_radius[2] * b43);
dL[3] = (temp_luminosity[0] * b41) + (temp_luminosity[1] * b42) + (temp_luminosity[2] * b43);
dM[3] = dm[i] * a4;
ptr_density = &temp_density[2];
ptr_energy = &temp_energy[2];
ptr_opacity = &temp_opacity[2];
old_temperature = temperature[i - direction];
old_pressure = pressure[i - direction];
rungekutta(direction, i, 3);
/*
Fifth Runge-Kutta loop
*/
dP[4] = (temp_pressure[0] * b51) + (temp_pressure[1] * b52) + (temp_pressure[2] * b53) + (temp_pressure[3] * b54);
dT[4] = (temp_temperature[0] * 51) + (temp_temperature[1] * b52) + (temp_temperature[2] * b53) + (temp_temperature[3] * b54);
dR[4] = (temp_radius[0] * b51) + (temp_radius[1] * b52) + (temp_radius[2] * b53) + (temp_radius[3] * b54);
dL[4] = (temp_luminosity[0] * b51) + (temp_luminosity[1] * b52) + (temp_luminosity[2] * b53) + (temp_luminosity[3] * b54);
dM[4] = dm[i] * a5;
ptr_density = &temp_density[3];
ptr_energy = &temp_energy[3];
ptr_opacity = &temp_opacity[3];
old_temperature = temperature[i - direction];
old_pressure = pressure[i - direction];
rungekutta(direction, i, 4);
/*
Sixth Runge-Kutta loop
*/
dP[5] = (temp_pressure[0] * b61) + (temp_pressure[1] * b62) + (temp_pressure[2] * b63) + (temp_pressure[3] * b64) + (temp_pressure[4] * b65);
dT[5] = (temp_temperature[0] * 61) + (temp_temperature[1] * b62) + (temp_temperature[2] * b63) + (temp_temperature[3] * b64) + (temp_temperature[4] * b65);
dR[5] = (temp_radius[0] * b61) + (temp_radius[1] * b62) + (temp_radius[2] * b63) + (temp_radius[3] * b64) + (temp_radius[4] * b65);
dL[5] = (temp_luminosity[0] * b61) + (temp_luminosity[1] * b62) + (temp_luminosity[2] * b63) + (temp_luminosity[3] * b64) + (temp_luminosity[4] * b65);
dM[5] = dm[i] * a6;
ptr_density = &temp_density[4];
ptr_energy = &temp_energy[4];
ptr_opacity = &temp_opacity[4];
old_temperature = temperature[i - direction];
old_pressure = pressure[i - direction];
rungekutta(direction, i, 5);
/*
Finally, compute the values.
*/
radius[i] = radius[i - direction] - direction * ((temp_radius[0] * c1) + (temp_radius[1] * c2) + (temp_radius[2] * c3) + (temp_radius[3] * c4) + (temp_radius[4] * c5) + (temp_radius[5] * c6));
pressure[i] = pressure[i - direction] + direction * ((temp_pressure[0] * c1) + (temp_pressure[1] * c2) + (temp_pressure[2] * c3) + (temp_pressure[3] * c4) + (temp_pressure[4] * c5) + (temp_pressure[5] * c6));
temperature[i] = temperature[i - direction] + direction * ((temp_temperature[0] * c1) + (temp_temperature[1] * c2) + (temp_temperature[2] * c3) + (temp_temperature[3] * c4) + (temp_temperature[4] * c5) + (temp_temperature[5] * c6));
luminosity[i] = luminosity[i - direction] - direction * ((temp_luminosity[0] * c1) + (temp_luminosity[1] * c2) + (temp_luminosity[2] * c3) + (temp_luminosity[3] * c4) + (temp_luminosity[4] * c5) + (temp_luminosity[5] * c6));
err_radius[i] = sqrt(pow(radius [i] - (radius[i - direction] - direction *((temp_radius[0] * C1) + (temp_radius[1] * C2) + (temp_radius[2] * C3) + (temp_radius[3] * C4) + (temp_radius[4] * C5) + (temp_radius[5] * C6))), 2.0));
err_pressure[i] = sqrt(pow(pressure[i] - (pressure[i - direction] + direction *((temp_pressure[0] * C1) + (temp_pressure[1] * C2) + (temp_pressure[2] * C3) + (temp_pressure[3] * C4) + (temp_pressure[4] * C5) + (temp_pressure[5] * C6))), 2.0));
err_temperature[i] = sqrt(pow(temperature[i] - (temperature[i - direction] + direction *((temp_temperature[0] * C1) + (temp_temperature[1] * C2) + (temp_temperature[2] * C3) + (temp_temperature[3] * C4) + (temp_temperature[4] * C5) + (temp_temperature[5] * C6))), 2.0));
err_luminosity[i] = sqrt(pow(luminosity[i] - (luminosity[i - direction] - direction *((temp_luminosity[0] * C1) + (temp_luminosity[1] * C2) + (temp_luminosity[2] * C3) + (temp_luminosity[3] * C4) + (temp_luminosity[4] * C5) + (temp_luminosity[5] * C6))), 2.0));
L1 = pow(((error * luminosity[i]) / err_luminosity[i]), 0.2);
T1 = pow(((error * temperature[i]) / err_temperature[i]), 0.2);
R1 = pow(((error * radius[i]) / err_radius[i]), 0.2);
P1 = pow(((error * pressure[i]) / err_pressure[i]), 0.2);
if (L1 < T1 && L1 < R1 && L1 < P1)
{
new_dm = dm[i] * L1;
}
else if (T1 > R1 && T1 > P1)
{
new_dm = dm[i] * T1;
}
else if (R1 > P1)
{
new_dm = dm[i] * R1;
}
else
{
new_dm = dm[i] * P1;
};
// Check and see which one is the greater - new_dm or min_dm - and use that one.
if (new_dm < min_dm)
{
dm[i + direction] = min_dm;
}
else
{
dm[i + direction] = new_dm;
};
density[i] = calc_density(X, temperature[i], pressure[i]);
opacity[i] = calc_opacity(X, density[i], temperature[i]);
energy[i] = calc_energy_generation(temperature[i], density[i], X, Y);
/*
Calculate the degeneracy pressure.
*/
fermi_pressure = pow((3.0 * density[i]) / (8.0 * CONST_PI * X * MASS_H), 1.0/3.0) * CONST_H;
if (fermi_pressure < MASS_e * 3.0e8)
{
degenerate_pressure[i] = (pow(CONST_H, 2.0) * pow(3.0 / CONST_PI, 2.0 / 3.0) * pow(density[i] / 2.0, 5.0 / 3.0)) / ( 20.0 * MASS_e * pow(MASS_H, 5.0 / 3.0));
}
else
{
degenerate_pressure[i] = (CONST_H * 3.0e8 * pow(3.0 / CONST_PI, 1.0 / 3.0) * pow(density[i] / 2.0, 4.0 / 3.0)) / (8.0 * pow (MASS_H, 4.0 / 3.0));
};
/*
Are we there yet? ("there" being half the mass of the star)
Alternatively; have things gone to nan? (Could check more than just temperature for nan, but as all the variables are interlinked there's not much point doing so)
*/
if (loop == 1 && (mass[i] > (mass[0] / 2.0) || isnan(temperature[i]) || temperature[i] == 0.0))
{
j = i;
break;
}
else if (loop != 1 && (mass[i] < (mass[0] / 2.0) || isnan(temperature[i]) || temperature[i] == 0.0))
{
k = i + 1;
break;
};
}; // End calculation loop
}; // End big loop
/*
Create the output file
*/
sprintf(filename, "%sOutput.csv", directory);
output_file = fopen(filename, "w");
if (!output_file)
{
printf("Failed to open %soutput.csv \n", directory);
return 1;
};
fprintf(output_file, "Number,Mass,Density,Radius,Radius Error,Pressure,Pressure Error,Degeneracy,Opacity,Temperature,Temperature Error,del_radiative,del_adiabatic,Energy,Luminosity,Luminosity Error\n");
/*
Output the calculated variables
*/
for (loop = 0; loop < 2; loop++)
{
if (loop == 0)
{
/*
This is the loop starting on the surface, and working inwards
Direction is positive as i is increasing - i=1 is the surface.
*/
direction = +1.0;
i_min = 0;
i_max = k;
}
else
{
/*
This is the loop starting at the core, and working outwards.
Direction is negative as i is decreasing - i=N is the center.
*/
direction = +1.0;
i_min = j;
i_max = N;
};
for (i = i_min; i < i_max; i = i + direction)
{
if (mass[i] != 0.0 || i > N-3)
{
fprintf(output_file, "%d,%.*Le,%.*Le,%.*Le,%.*Le,%.*Le,%.*Le,%.*Le,%.*Le,%.*Le,%.*Le,%.*Le,%.*Le,%.*Le,%.*Le,%.*Le\n", i, accuracy, mass[i], accuracy, density[i], accuracy, radius[i], accuracy, err_radius[i], accuracy, pressure[i], accuracy, err_pressure[i], accuracy, degenerate_pressure[i], accuracy, opacity[i], accuracy, temperature[i], accuracy, err_temperature[i], accuracy, del_radiative[i], accuracy, del_adiabatic[i], accuracy, energy[i], accuracy, luminosity[i], accuracy, err_luminosity[i]);
};
};
};
fclose(output_file);
/*
All is done, so end the program.
*/
printf("\n%d positions calculated (%d inwards, %d outwards)\n", (N - j) + k, k, (N - j));
return 0;
}
