/**********************************************************************/

/**********************************************************************/

#include<math.h>

#include<cstring.h>

#include<vector.h>

#include<iostream.h>

#include<fstream.h>

#include<iomanip.h>

double new_eqn_of_motion(double);

double R[600], V[600];

vector <double> Period, R_equil;

string output_file = "Cepheids.nja";

float T; //period of oscillation

float dt = 10000; //s

float t_elapsed; //s

float R_eq; //Equilibrium radius. 1.5e10 is a standard value

float Ro; //Initial (perturbed) radius

const float M = 1e31; //typical mass of Cepheid 1e31kg

const float temp = 6500; //assumed Cepheid temp

const float sigma = 5.670e-8; //Stefan-Boltzmann constant

const float Mag_sun = 4.80; //

const float L_sun = 3.83e26; //watts

int main (){

cout << "Cepheids: "

<< "\n\nCalculates oscillation data for a range of stars of "

<< "different radii. \nOutput file is Cepheid.nja\n";

for (int c = 1; c <= 60; c++){

R_eq = c*1e9;//corresponds to 1 M km to 60 M km. Typical value 15Mkm

Ro = 1.1*R_eq; //perturb each star by +10% of equilibrium radius

t_elapsed = 0;

T = 0;

//initialise arrays

R[0] = Ro;

V[0] = 0;

int i = 0;

int j = 1; // j acts as R(i+1)

int expand = 0;//expand=1 if init. expansion,-1 if init. contraction

int state = 0; //state is a 3-valued sentinel to determine period

float t1 = 0;

//main calculation loop

while (t_elapsed <= 5e6 && state!=3){//5e6s=2months

V[j] = V[i] + new_eqn_of_motion(R[i])*dt;

R[j] = R[i] + V[j]*dt;

//does star initially expand or contract?

if (R[1]>R[0])expand = 1; else expand = -1;

//calculate period - if any - of initially expanding star

if (expand ==1){

if (state==0 && R[j]<=R[i]){t1=t_elapsed; state = 1;}

if (state==1 && R[j] >= R[i]){state = 2;}

if (state==2 && R[j]<= R[i]){state = 3; T = (t_elapsed - t1);}

}

//calculate period - if any - of initially contracting star

if (expand ==-1){

if (state==0 && R[j]>=R[i]){t1=t_elapsed; state = 1;}

if (state==1 && R[j]<= R[i]){state = 2;}

if (state==2 && R[j]>= R[i]){state = 3; T =(t_elapsed - t1);}

}

t_elapsed = t_elapsed + dt;

i++; j++;

}

R_equil.push_back(R_eq);

Period.push_back(T);

}

//prepare for file-write and screen output

ofstream write(output_file.c_str());

write.setf(ios::fixed);

write.precision(3);

cout.setf(ios::fixed);

cout.precision(3);

//screen headings

cout << "\nR(eq)/Mkm" << setw(15) << "T/days"<< setw(15) << "log(T/days)"

<< setw(15) << "Mag";

//file headings

write<< "R(eq)/Mkm," << setw(15) << "T/days," << setw(15)<< "log(T/days),"

<< setw(15) << "Mag,";

//write to screen and file

for(c = 0; c < 60; c++){

if(Period[c] > 0){

cout<<"\n"<<R_equil[c]/1e9<<setw(15)<<Period[c]/86400 <<

setw(15) << log10(Period[c]/86400) << setw(15) <<

(-2.5*log10( 4*M_PI*pow(R_equil[c],2)*sigma*pow(temp,4)/L_sun)+Mag_sun);

write << "\n" << R_equil[c]/1e9 << "," << setw(15)

<< Period[c]/86400 << "," << setw(15) << log10(Period[c]/86400) << ","

<< setw(15) <<

(-2.5*log10( 4*M_PI*pow(R_equil[c],2)*sigma*pow(temp,4)/L_sun)+Mag_sun);

}

}

return 0;

}

double new_eqn_of_motion(double R){

return 6.7e-11*M*(R_eq*pow(R,-3)-pow(R,-2));

}