00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "phycon.h"
00007 #include "thermal.h"
00008 #include "dense.h"
00009 #include "thirdparty.h"
00010 #include "atoms.h"
00011
00012
00013 void atom_pop5(
00014
00015 const double g[],
00016
00017
00018
00019 const double EnerWN[],
00020
00021 double cs12,
00022 double cs13,
00023 double cs14,
00024 double cs15,
00025 double cs23,
00026 double cs24,
00027 double cs25,
00028 double cs34,
00029 double cs35,
00030 double cs45,
00031
00032 double a21,
00033 double a31,
00034 double a41,
00035 double a51,
00036 double a32,
00037 double a42,
00038 double a52,
00039 double a43,
00040 double a53,
00041 double a54,
00042
00043 double p[],
00044
00045 realnum abund,
00046 double *Cooling,
00047 double *CoolingDeriv,
00048 double pump01,
00049 double pump02,
00050 double pump03,
00051 double pump04
00052 )
00053 {
00054
00055 DEBUG_ENTRY( "atom_pop5()" );
00056
00057
00058 ASSERT( abund>=0. );
00059 if( abund == 0. )
00060 {
00061 p[0] = 0.;
00062 p[1] = 0.;
00063 p[2] = 0.;
00064 p[3] = 0.;
00065 p[4] = 0.;
00066 *Cooling = 0.;
00067 *CoolingDeriv = 0.;
00068 return;
00069 }
00070
00071
00072 double tf = T1CM/phycon.te;
00073
00074
00075 double BoltzFac[5][5];
00076
00077 BoltzFac[0][1] = sexp(EnerWN[0]*tf);
00078 BoltzFac[1][2] = sexp(EnerWN[1]*tf);
00079 BoltzFac[2][3] = sexp(EnerWN[2]*tf);
00080 BoltzFac[3][4] = sexp(EnerWN[3]*tf);
00081 BoltzFac[0][2] = BoltzFac[0][1]*BoltzFac[1][2];
00082 BoltzFac[0][3] = BoltzFac[0][2]*BoltzFac[2][3];
00083 BoltzFac[0][4] = BoltzFac[0][3]*BoltzFac[3][4];
00084 BoltzFac[1][3] = BoltzFac[1][2]*BoltzFac[2][3];
00085 BoltzFac[1][4] = BoltzFac[1][3]*BoltzFac[3][4];
00086 BoltzFac[2][4] = BoltzFac[2][3]*BoltzFac[3][4];
00087
00088
00089 if( (BoltzFac[0][4]+pump04) == 0. )
00090 {
00091 p[0] = 0.;
00092 p[1] = 0.;
00093 p[2] = 0.;
00094 p[3] = 0.;
00095 p[4] = 0.;
00096 *Cooling = 0.;
00097 *CoolingDeriv = 0.;
00098 return;
00099 }
00100
00101
00102
00103 double col[5][5];
00104 col[1][0] = dense.cdsqte*cs12/g[1];
00105 col[0][1] = col[1][0]*g[1]/g[0]*BoltzFac[0][1];
00106
00107 col[2][0] = dense.cdsqte*cs13/g[2];
00108 col[0][2] = col[2][0]*g[2]/g[0]*BoltzFac[0][2];
00109
00110 col[3][0] = dense.cdsqte*cs14/g[3];
00111 col[0][3] = col[3][0]*g[3]/g[0]*BoltzFac[0][3];
00112
00113 col[4][0] = dense.cdsqte*cs15/g[4];
00114 col[0][4] = col[4][0]*g[4]/g[0]*BoltzFac[0][4];
00115
00116 col[2][1] = dense.cdsqte*cs23/g[2];
00117 col[1][2] = col[2][1]*g[2]/g[1]*BoltzFac[1][2];
00118
00119 col[3][1] = dense.cdsqte*cs24/g[3];
00120 col[1][3] = col[3][1]*g[3]/g[1]*BoltzFac[1][3];
00121
00122 col[4][1] = dense.cdsqte*cs25/g[4];
00123 col[1][4] = col[4][1]*g[4]/g[1]*BoltzFac[1][4];
00124
00125 col[3][2] = dense.cdsqte*cs34/g[3];
00126 col[2][3] = col[3][2]*g[3]/g[2]*BoltzFac[2][3];
00127
00128 col[4][2] = dense.cdsqte*cs35/g[4];
00129 col[2][4] = col[4][2]*g[4]/g[2]*BoltzFac[2][4];
00130
00131 col[4][3] = dense.cdsqte*cs45/g[4];
00132 col[3][4] = col[4][3]*g[4]/g[3]*BoltzFac[3][4];
00133
00134 double amat[5][5], bvec[5];
00135
00136 for( long i=0; i<5; ++i )
00137 {
00138 amat[i][4] = 1.;
00139 bvec[i] = 0.;
00140 }
00141 bvec[4] = abund;
00142
00143
00144 amat[0][0] = col[0][1] + col[0][2] + col[0][3] + col[0][4] +pump01+pump02+pump03+pump04;
00145 amat[1][0] = -a21 - col[1][0];
00146 amat[2][0] = -a31 - col[2][0];
00147 amat[3][0] = -a41 - col[3][0];
00148 amat[4][0] = -a51 - col[4][0];
00149
00150
00151 amat[0][1] = -col[0][1] - pump01;
00152 amat[1][1] = col[1][0] + a21 + col[1][2] + col[1][3] + col[1][4];
00153 amat[2][1] = -col[2][1] - a32;
00154 amat[3][1] = -col[3][1] - a42;
00155 amat[4][1] = -col[4][1] - a52;
00156
00157
00158 amat[0][2] = -col[0][2] - pump02;
00159 amat[1][2] = -col[1][2];
00160 amat[2][2] = a31 + a32 + col[2][0] + col[2][1] + col[2][3] + col[2][4];
00161 amat[3][2] = -col[3][2] - a43;
00162 amat[4][2] = -col[4][2] - a53;
00163
00164
00165 amat[0][3] = -col[0][3] - pump03;
00166 amat[1][3] = -col[1][3];
00167 amat[2][3] = -col[2][3];
00168 amat[3][3] = a41 + col[3][0] + a42 + col[3][1] + a43 + col[3][2] + col[3][4];
00169 amat[4][3] = -col[4][3] - a54;
00170
00171 # if 0
00172
00173
00174 double dmax = -1e0;
00175 for( i=0; i < 6; i++ )
00176 for( j=0; j < 5; j++ )
00177 dmax = MAX2(zz[i][j],dmax);
00178
00179 for( i=0; i < 6; i++ )
00180 for( j=0; j < 5; j++ )
00181 zz[i][j] /= dmax;
00182 # endif
00183
00184 int32 ipiv[5], ner=0;
00185
00186 getrf_wrapper(5,5,(double*)amat,5,ipiv,&ner);
00187 getrs_wrapper('N',5,1,(double*)amat,5,ipiv,bvec,5,&ner);
00188
00189 if( ner != 0 )
00190 {
00191 fprintf( ioQQQ, "DISASTER PROBLEM atom_pop5: dgetrs finds singular or ill-conditioned matrix\n" );
00192 cdEXIT(EXIT_FAILURE);
00193 }
00194
00195
00196
00197 p[1] = MAX2(0.e0,bvec[1]);
00198 p[2] = MAX2(0.e0,bvec[2]);
00199 p[3] = MAX2(0.e0,bvec[3]);
00200 p[4] = MAX2(0.e0,bvec[4]);
00201
00202 p[0] = abund - p[1] - p[2] - p[3] - p[4];
00203
00204
00205 double Erg[5] , EnergyKelvin[5];
00206 Erg[0] = 0.;
00207 EnergyKelvin[0] = 0.;
00208 for( long i=1; i<5; ++i )
00209 {
00210 Erg[i] = Erg[i-1] + EnerWN[i-1]*ERG1CM;
00211 EnergyKelvin[i] = EnergyKelvin[i-1] + EnerWN[i-1]*T1CM;
00212 }
00213
00214 *Cooling = 0.;
00215 *CoolingDeriv = 0.;
00216 for( long ihi=1; ihi<5; ++ihi )
00217 {
00218 for( long ilo=0; ilo<ihi; ++ilo )
00219 {
00220 double CoolOne = (p[ilo]*col[ilo][ihi] - p[ihi]*col[ihi][ilo])*
00221 (Erg[ihi]-Erg[ilo]);
00222
00223 *Cooling += CoolOne;
00224 *CoolingDeriv += CoolOne*( EnergyKelvin[ihi]*thermal.tsq1 - thermal.halfte );
00225 }
00226 }
00227
00228 return;
00229 }