19 const double EnerWN[],
75 double BoltzFac[5][5];
77 BoltzFac[0][1] =
sexp(EnerWN[0]*tf);
78 BoltzFac[1][2] =
sexp(EnerWN[1]*tf);
79 BoltzFac[2][3] =
sexp(EnerWN[2]*tf);
80 BoltzFac[3][4] =
sexp(EnerWN[3]*tf);
81 BoltzFac[0][2] = BoltzFac[0][1]*BoltzFac[1][2];
82 BoltzFac[0][3] = BoltzFac[0][2]*BoltzFac[2][3];
83 BoltzFac[0][4] = BoltzFac[0][3]*BoltzFac[3][4];
84 BoltzFac[1][3] = BoltzFac[1][2]*BoltzFac[2][3];
85 BoltzFac[1][4] = BoltzFac[1][3]*BoltzFac[3][4];
86 BoltzFac[2][4] = BoltzFac[2][3]*BoltzFac[3][4];
89 if( (BoltzFac[0][4]+pump04) == 0. )
104 col[1][0] =
dense.cdsqte*cs12/
g[1];
105 col[0][1] = col[1][0]*
g[1]/
g[0]*BoltzFac[0][1];
107 col[2][0] =
dense.cdsqte*cs13/
g[2];
108 col[0][2] = col[2][0]*
g[2]/
g[0]*BoltzFac[0][2];
110 col[3][0] =
dense.cdsqte*cs14/
g[3];
111 col[0][3] = col[3][0]*
g[3]/
g[0]*BoltzFac[0][3];
113 col[4][0] =
dense.cdsqte*cs15/
g[4];
114 col[0][4] = col[4][0]*
g[4]/
g[0]*BoltzFac[0][4];
116 col[2][1] =
dense.cdsqte*cs23/
g[2];
117 col[1][2] = col[2][1]*
g[2]/
g[1]*BoltzFac[1][2];
119 col[3][1] =
dense.cdsqte*cs24/
g[3];
120 col[1][3] = col[3][1]*
g[3]/
g[1]*BoltzFac[1][3];
122 col[4][1] =
dense.cdsqte*cs25/
g[4];
123 col[1][4] = col[4][1]*
g[4]/
g[1]*BoltzFac[1][4];
125 col[3][2] =
dense.cdsqte*cs34/
g[3];
126 col[2][3] = col[3][2]*
g[3]/
g[2]*BoltzFac[2][3];
128 col[4][2] =
dense.cdsqte*cs35/
g[4];
129 col[2][4] = col[4][2]*
g[4]/
g[2]*BoltzFac[2][4];
131 col[4][3] =
dense.cdsqte*cs45/
g[4];
132 col[3][4] = col[4][3]*
g[4]/
g[3]*BoltzFac[3][4];
134 double amat[5][5], bvec[5];
136 for(
long i=0; i<5; ++i )
144 amat[0][0] = col[0][1] + col[0][2] + col[0][3] + col[0][4] +pump01+pump02+pump03+pump04;
145 amat[1][0] = -a21 - col[1][0];
146 amat[2][0] = -a31 - col[2][0];
147 amat[3][0] = -a41 - col[3][0];
148 amat[4][0] = -a51 - col[4][0];
151 amat[0][1] = -col[0][1] - pump01;
152 amat[1][1] = col[1][0] + a21 + col[1][2] + col[1][3] + col[1][4];
153 amat[2][1] = -col[2][1] - a32;
154 amat[3][1] = -col[3][1] - a42;
155 amat[4][1] = -col[4][1] - a52;
158 amat[0][2] = -col[0][2] - pump02;
159 amat[1][2] = -col[1][2];
160 amat[2][2] = a31 + a32 + col[2][0] + col[2][1] + col[2][3] + col[2][4];
161 amat[3][2] = -col[3][2] - a43;
162 amat[4][2] = -col[4][2] - a53;
165 amat[0][3] = -col[0][3] - pump03;
166 amat[1][3] = -col[1][3];
167 amat[2][3] = -col[2][3];
168 amat[3][3] = a41 + col[3][0] + a42 + col[3][1] + a43 + col[3][2] + col[3][4];
169 amat[4][3] = -col[4][3] - a54;
175 for( i=0; i < 6; i++ )
176 for( j=0; j < 5; j++ )
177 dmax =
MAX2(zz[i][j],dmax);
179 for( i=0; i < 6; i++ )
180 for( j=0; j < 5; j++ )
184 int32 ipiv[5], ner=0;
191 fprintf(
ioQQQ,
"DISASTER PROBLEM atom_pop5: dgetrs finds singular or ill-conditioned matrix\n" );
197 p[1] =
MAX2(0.e0,bvec[1]);
198 p[2] =
MAX2(0.e0,bvec[2]);
199 p[3] =
MAX2(0.e0,bvec[3]);
200 p[4] =
MAX2(0.e0,bvec[4]);
202 p[0] =
abund - p[1] - p[2] - p[3] - p[4];
205 double Erg[5] , EnergyKelvin[5];
207 EnergyKelvin[0] = 0.;
208 for(
long i=1; i<5; ++i )
210 Erg[i] = Erg[i-1] + EnerWN[i-1]*
ERG1CM;
211 EnergyKelvin[i] = EnergyKelvin[i-1] + EnerWN[i-1]*
T1CM;
216 for(
long ihi=1; ihi<5; ++ihi )
218 for(
long ilo=0; ilo<ihi; ++ilo )
220 double CoolOne = (p[ilo]*col[ilo][ihi] - p[ihi]*col[ihi][ilo])*
224 *CoolingDeriv += CoolOne*( EnergyKelvin[ihi]*
thermal.tsq1 -
thermal.halfte );
void atom_pop5(const double g[], const double EnerWN[], double cs12, double cs13, double cs14, double cs15, double cs23, double cs24, double cs25, double cs34, double cs35, double cs45, double a21, double a31, double a41, double a51, double a32, double a42, double a52, double a43, double a53, double a54, double p[], realnum abund, double *Cooling, double *CoolingDeriv, double pump01, double pump02, double pump03, double pump04)
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)