1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
|
double mu;
if (ok_surface){
//Load array of first derivatives to start the general inward integration
Y = Helium(X, Z);
mu = Mean_Molecular_Weight(X, Y, Z);
gamma = Specific_Heat_Ratio();
dlnPdlnT = PTgradient(Pm, P, Tm, T);
dfdr0[0] = dPdr(M_r, rho, r);
dfdr0[1] = dMdr(r, rho);
dfdr0[2] = dLdr(r, rho, epsilon);
dfdr0[3] = dTdr(kappa, rho, T, L_r, r, mu, M_r, gamma, dlnPdlnT, rc_flag);
//Main inward integration loop
bool exit_main = false;
while (!exit_main){
i++;
//Update last zone values
Mm = M_r;
Lm = L_r;
Pm = P;
Tm = T;
Xm = X;
Zm = Z;
rm = r;
taum = tau;
rhom = rho;
kappam = kappa;
epsilonm = epsilon;
double PMLT0[4] = {Pm, Mm, Lm, Tm};
double PMLT[4];
RK_4(n, rm, dr, PMLT0, PMLT, dfdr0, ok_Runge, X, Z, Pm, Tm, step_size_condition, rc_flag);
if (!ok_Runge) exit_main = true;
if (ok_Runge){
//Results from the current step
P = PMLT[0];
M_r = PMLT[1];
L_r = PMLT[2];
T = PMLT[3];
//Compute current step quantities for output and next step
Y = Helium(X, Z);
mu = Mean_Molecular_Weight(X, Y, Z);
rho = Density(T, P, mu, step_size_condition);
kappa = Opacity(T, rho, X, Z);
dlnPdlnT = PTgradient(Pm, P, Tm, T);
epsilon = Nuclear(T, rho, X, Z);
dfdr0[0] = dPdr(M_r, rho, r);
dfdr0[1] = dMdr(r, rho);
dfdr0[2] = dLdr(r, rho, epsilon);
dfdr0[3] = dTdr(kappa, rho, T, L_r, r, mu, M_r, gamma, dlnPdlnT, rc_flag);
tau = taum + Optical_Depth_Change(kappa, kappam, rho, rhom, rm + dr, rm);
outdata << setw(5) << i
<< scientific << showpoint << right << setprecision(3)
<< setw(11) << r
<< setw(11) << tau
<< setw(11) << 1 - M_r/Ms
<< setw(11) << L_r
<< setw(11) << T
<< setw(11) << P
<< setw(11) << rho
<< setw(11) << kappa
<< setw(11) << epsilon
<< " " << rc_flag
<< fixed << showpoint << right << setprecision(1)
<< setw(5) << dlnPdlnT << endl;
if ((M_r/Ms < M_fraction_limit && L_r/Ls < L_fraction_limit && r/Rs < r_fraction_limit)
|| T < 0 || P < 0) exit_main = true;
else if (i > maximum_zones) {
Too_Many_Zones(i, Msolar, Ms, M_r, Lsolar, Ls, L_r, r, Rs, Rsolar, Teff, X, Y, Z, P_0, T_0,
rho_0, kappa_0, epsilon_0, rc_flag);
ok_Runge = false;
exit_main = true;}
if (!exit_main){
//Is it time to change step size?
if (adjust_step_size){
switch (step_size_condition){
case 0:
if (M_r < 0.99*Ms) {
dr = -Rs/100;
step_size_condition = 1;}
break;
case 1:
if (fabs(dr) > 5*r) {
dr /= 10;
step_size_condition = 2;}
break;
}
}
}
}
r += dr;
}
//Core_Extrapolation
if (ok_Runge) {
//Determine core conditions
i++;
Core(M_r, L_r, P, T, X, Z, r, P_0, T_0, rho_0, kappa_0, epsilon_0, rc_flag, dlnPdlnT, ok_core);
if (!ok_core) {
cout << "\nWARNING: There was a problem with the core extrapolation routine" << endl;}
tau += Optical_Depth_Change(kappa_0, kappa, rho_0, rho, 0, r);
outdata << setw(5) << i
<< scientific << showpoint << right << setprecision(3)
<< setw(11) << 0
<< setw(11) << tau
<< setw(11) << 1 - M_r/Ms
<< setw(11) << L_r
<< setw(11) << T_0
<< setw(11) << P_0
<< setw(11) << rho_0
<< setw(11) << kappa_0
<< setw(11) << epsilon_0
<< " " << rc_flag
<< fixed << showpoint << right << setprecision(1)
<< setw(5) << dlnPdlnT << endl;
//Write initial and final conditions to the screen
Final_Results(i, Msolar, Ms, M_r, Lsolar, Ls, L_r, r, Rs, Rsolar, Teff, X, Y, Z,
P, T, rho, kappa, epsilon, P_0, T_0, rho_0, kappa_0, epsilon_0, rc_flag);
}
}
//Does the user want to compute a new model?
all_new = 'Y';
Change_Model(Msolar, Lsolar, Rsolar, Ms, Ls, Rs, Teff, X, Y, Z, all_new);
outdata.close();
if (outdata.fail()){
cout << " Unable to close the output file - the new model calculation is being aborted"
<< "\n\n"
<< "Enter any character to quit: ";
cin >> xpause;
exit(1);
}
}
return 0;
}
|