00001
00002
00003
00004 #include "cddefines.h"
00005 #include "phycon.h"
00006 #include "transition.h"
00007 #include "thirdparty.h"
00008 #include "dense.h"
00009 #include "cooling.h"
00010 #include "thermal.h"
00011 #include "atoms.h"
00012
00013 void AtomSeqBeryllium(double cs12,
00014 double cs13,
00015 double cs23,
00016 const TransitionProxy& t,
00017 double a30)
00018 {
00019
00020 char chLab[5];
00021
00022 bool lgNegPop;
00023
00024 int32 ipiv[4], nerror;
00025 long int i, j;
00026
00027 double AbunxIon,
00028 Enr01,
00029 Enr10,
00030 a20,
00031 boltz,
00032 c01,
00033 c02,
00034 c03,
00035 c10,
00036 c12,
00037 c13,
00038 c20,
00039 c21,
00040 c23,
00041 c30,
00042 c31,
00043 c32,
00044 coolng,
00045 cs1u,
00046 excit,
00047 r01,
00048 r02,
00049 r03,
00050 r10,
00051 r12,
00052 r13,
00053 r20,
00054 r21,
00055 r23,
00056 r30,
00057 r31,
00058 r32,
00059 ri02,
00060 ri20,
00061 tot20;
00062
00063 double amat[4][4],
00064 bvec[4],
00065 zz[5][5];
00066
00067 DEBUG_ENTRY( "AtomSeqBeryllium()" );
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086 AbunxIon = dense.xIonDense[(*t.Hi()).nelem() -1][(*t.Hi()).IonStg()-1];
00087
00088
00089 chIonLbl(chLab, t);
00090 boltz = t.EnergyK()/phycon.te;
00091
00092
00093 cs1u = t.Coll().col_str();
00094
00095
00096 t.Coll().col_str() /= 3.f;
00097
00098
00099 if( AbunxIon <= 1e-20 || boltz > 30. )
00100 {
00101
00102 (*t.Lo()).Pop() = AbunxIon;
00103 t.Emis().PopOpc() = AbunxIon;
00104 (*t.Hi()).Pop() = 0.;
00105 t.Emis().xIntensity() = 0.;
00106 t.Coll().cool() = 0.;
00107 t.Emis().phots() = 0.;
00108
00109
00110 t.Emis().ColOvTot() = 0.;
00111 t.Coll().heat() = 0.;
00112
00113 atoms.PopLevels[0] = AbunxIon;
00114 atoms.PopLevels[1] = 0.;
00115 atoms.PopLevels[2] = 0.;
00116 atoms.PopLevels[3] = 0.;
00117 atoms.DepLTELevels[0] = 1.;
00118 atoms.DepLTELevels[1] = 0.;
00119 atoms.DepLTELevels[2] = 0.;
00120 atoms.DepLTELevels[3] = 0.;
00121 CoolAdd(chLab, t.WLAng() ,0.);
00122 return;
00123 }
00124 excit = exp(-boltz);
00125
00126
00127 ASSERT( (*t.Lo()).g() == 1. );
00128 ASSERT( (*t.Hi()).g() == 3. );
00129
00130
00131 ASSERT( cs1u > 0. );
00132
00133
00134 ri02 = t.Emis().pump();
00135
00136
00137 ri20 = t.Emis().pump()*1./3.;
00138
00139
00140 a20 = t.Emis().Aul()*(t.Emis().Pesc() + t.Emis().Pelec_esc() + t.Emis().Pdest());
00141 tot20 = a20 + ri20;
00142
00143
00144
00145 c10 = cs1u*dense.cdsqte/9.;
00146 c01 = c10*excit;
00147 r01 = c01;
00148 r10 = c10;
00149
00150
00151 c20 = c10;
00152 c02 = c01*3.;
00153 r02 = c02 + ri02;
00154 r20 = c20 + tot20;
00155
00156 c30 = c10;
00157 c03 = c01*5.;
00158 r30 = c30 + a30;
00159 r03 = c03;
00160
00161 c21 = cs12*dense.cdsqte/3.;
00162 c12 = c21*3.;
00163 r12 = c12;
00164 r21 = c21;
00165
00166 c31 = cs13*dense.cdsqte/5.;
00167 c13 = c31*5.;
00168 r31 = c31;
00169 r13 = c13;
00170
00171 c32 = cs23*dense.cdsqte/5.;
00172 c23 = c32*1.667;
00173 r32 = c32;
00174 r23 = c23;
00175
00176
00177 for( i=0; i <= 3; i++ )
00178 {
00179
00180 zz[i][0] = 1.;
00181 zz[4][i] = 0.;
00182 }
00183
00184
00185 zz[4][0] = 1.;
00186
00187
00188
00189 zz[0][1] = -r01;
00190 zz[1][1] = r10 + r12 + r13;
00191 zz[2][1] = -r21;
00192 zz[3][1] = -r31;
00193
00194
00195 zz[0][2] = -r02;
00196 zz[1][2] = -r12;
00197 zz[2][2] = r20 + r21 + r23;
00198 zz[3][2] = -r32;
00199
00200
00201 zz[0][3] = -r03;
00202 zz[1][3] = -r13;
00203 zz[2][3] = -r23;
00204 zz[3][3] = r30 + r31 + r32;
00205
00206
00207 for( j=0; j <= 3; j++ )
00208 {
00209 for( i=0; i <= 3; i++ )
00210 {
00211 amat[i][j] = zz[i][j];
00212 }
00213 bvec[j] = zz[3+1][j];
00214 }
00215
00216 nerror = 0;
00217
00218 getrf_wrapper(4, 4, (double*)amat, 4, ipiv, &nerror);
00219 getrs_wrapper('N', 4, 1, (double*)amat, 4, ipiv, bvec, 4, &nerror);
00220
00221 if( nerror != 0 )
00222 {
00223 fprintf( ioQQQ, " AtomSeqBeryllium: dgetrs finds singular or ill-conditioned matrix\n" );
00224 cdEXIT(EXIT_FAILURE);
00225 }
00226
00227
00228 for( i=0; i <= 3; i++ )
00229 {
00230 zz[3+1][i] = bvec[i];
00231 }
00232
00233 lgNegPop = false;
00234 for( i=0; i <= 3; i++ )
00235 {
00236 atoms.PopLevels[i] = zz[4][i]*AbunxIon;
00237 if( atoms.PopLevels[i] < 0. )
00238 lgNegPop = true;
00239 }
00240
00241
00242 if( lgNegPop )
00243 {
00244 fprintf( ioQQQ, " AtomSeqBeryllium finds non-positive pop,=" );
00245 for( i=0; i <= 3; i++ )
00246 {
00247 fprintf( ioQQQ, "%g ", atoms.PopLevels[i] );
00248 }
00249 fprintf( ioQQQ, "%s \n", chLab );
00250 fprintf( ioQQQ, " te=%g total abund=%g boltz=%g \n",
00251 phycon.te, AbunxIon, boltz );
00252 cdEXIT(EXIT_FAILURE);
00253 }
00254
00255
00256 atoms.DepLTELevels[0] = 1.;
00257 atoms.DepLTELevels[1] = (atoms.PopLevels[1]/atoms.PopLevels[0])/
00258 excit;
00259 atoms.DepLTELevels[2] = (atoms.PopLevels[2]/atoms.PopLevels[0])/
00260 (excit*3.);
00261 atoms.DepLTELevels[3] = (atoms.PopLevels[3]/atoms.PopLevels[0])/
00262 (excit*5.);
00263
00264
00265 t.Emis().ColOvTot() = c02/r02;
00266
00267 (*t.Lo()).Pop() = AbunxIon;
00268
00269
00270 (*t.Hi()).Pop() = atoms.PopLevels[2];
00271
00272 t.Emis().PopOpc() = AbunxIon - atoms.PopLevels[2]*1./3.;
00273
00274
00275
00276 t.Emis().phots() = atoms.PopLevels[2] * t.Emis().Aul() * ( t.Emis().Pesc() + t.Emis().Pelec_esc() );
00277
00278 t.Emis().xIntensity() = t.Emis().phots() * t.EnergyErg();
00279
00280
00281
00282
00283
00284
00285
00286 Enr01 = atoms.PopLevels[0]*(c01 + c02 + c03)*t.EnergyErg();
00287 Enr10 = (atoms.PopLevels[1]*c10 + atoms.PopLevels[2]*c20 +
00288 atoms.PopLevels[3]*c30)*t.EnergyErg();
00289
00290
00291 t.Coll().cool() = Enr01 - Enr10*t.Emis().ColOvTot();
00292
00293
00294 t.Coll().heat() = Enr10*(1. - t.Emis().ColOvTot());
00295
00296
00297 coolng = t.Coll().cool();
00298 CoolAdd(chLab, t.WLAng() ,coolng);
00299
00300
00301 thermal.dCooldT += coolng*(t.EnergyK()*thermal.tsq1 - thermal.halfte);
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311 return;
00312 }