14 int nlev, ind_start,
i,
j,
k, nv, nrt;
15 double d, n_el, n_el_old, sT;
23 int ind_starts[] = {
X_HI,
X_HeI, X_CI, X_NI, X_OI, X_NeI, X_SI, X_SI+
S_IONS};
41 for (nv = 0; nv < 8; nv++) ind_starts[nv] -=
NFLX;
49 for (j = 0; j <
NIONS; j++) v[
NFLX+j] = 0.0;
56 if (T > 5000.) n_el = N*0.1;
68 while ( fabs(n_el - n_el_old) / n_el_old > 1.e-5 && nrt<500 && n_el > 1.0e-14) {
73 for (k = 1; k < 8; k++) {
75 ind_start = ind_starts[
k];
77 if (nlev == 0)
continue;
87 for (j = 0 ; j < nlev-1 ; j++) {
98 for (j = 0; j < nlev; j++) {
111 for (i = 0 ; i < nlev ; i++) {
112 if (rhs[i] >= 0. && rhs[i] <= 1.0) v[
NFLX+i+ind_start] = rhs[
i];
113 else if (rhs[i] < 0.0) v[
NFLX + i + ind_start] = 0.0;
114 else v[
NFLX + i + ind_start] = 1.0;
126 print (
"\n!CompEquil: GR(H) = %12.6e\nDR(H) = %12.6e\n\n",
136 for (nv = 0; nv <
NIONS; nv++) {
141 print1 (
"! CompEquil: error!! n_el NaN\n\n");
void print1(const char *fmt,...)
void CompEquil(double n, double T, double *v)
const double rad_rec_z[31]
void LUBackSubst(double **a, int n, int *indx, double b[])
void print(const char *fmt,...)
#define ARRAY_1D(nx, type)
void Find_Rates(double T, double Ne, double N, double *v)
#define ARRAY_2D(nx, ny, type)
int LUDecompose(double **a, int n, int *indx, double *d)
#define QUIT_PLUTO(e_code)