00001
00002
00003
00004 #include "cddefines.h"
00005 #include "phycon.h"
00006 #include "lines_service.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 EmLine * 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->nelem -1][t->IonStg-1];
00087
00088
00089 chIonLbl(chLab, t);
00090 boltz = t->EnergyK/phycon.te;
00091
00092
00093 cs1u = t->cs;
00094
00095
00096 t->cs /= 3.f;
00097
00098
00099 if( AbunxIon <= 1e-20 || boltz > 30. )
00100 {
00101
00102 t->PopLo = AbunxIon;
00103 t->PopOpc = AbunxIon;
00104 t->PopHi = 0.;
00105 t->xIntensity = 0.;
00106 t->cool = 0.;
00107 t->phots = 0.;
00108
00109
00110 t->ColOvTot = 0.;
00111 t->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
00123 DEBUG_EXIT( "AtomSeqBeryllium()" );
00124 return;
00125 }
00126 excit = exp(-boltz);
00127
00128
00129 ASSERT( t->gLo == 1. );
00130 ASSERT( t->gHi == 3. );
00131
00132
00133 ASSERT( cs1u > 0. );
00134
00135
00136 ri02 = t->pump;
00137
00138
00139 ri20 = t->pump*1./3.;
00140
00141
00142 a20 = t->Aul*(t->Pesc + t->Pelec_esc + t->Pdest);
00143 tot20 = a20 + ri20;
00144
00145
00146
00147 c10 = cs1u*dense.cdsqte/9.;
00148 c01 = c10*excit;
00149 r01 = c01;
00150 r10 = c10;
00151
00152
00153 c20 = c10;
00154 c02 = c01*3.;
00155 r02 = c02 + ri02;
00156 r20 = c20 + tot20;
00157
00158 c30 = c10;
00159 c03 = c01*5.;
00160 r30 = c30 + a30;
00161 r03 = c03;
00162
00163 c21 = cs12*dense.cdsqte/3.;
00164 c12 = c21*3.;
00165 r12 = c12;
00166 r21 = c21;
00167
00168 c31 = cs13*dense.cdsqte/5.;
00169 c13 = c31*5.;
00170 r31 = c31;
00171 r13 = c13;
00172
00173 c32 = cs23*dense.cdsqte/5.;
00174 c23 = c32*1.667;
00175 r32 = c32;
00176 r23 = c23;
00177
00178
00179 for( i=0; i <= 3; i++ )
00180 {
00181
00182 zz[i][0] = 1.;
00183 zz[4][i] = 0.;
00184 }
00185
00186
00187 zz[4][0] = 1.;
00188
00189
00190
00191 zz[0][1] = -r01;
00192 zz[1][1] = r10 + r12 + r13;
00193 zz[2][1] = -r21;
00194 zz[3][1] = -r31;
00195
00196
00197 zz[0][2] = -r02;
00198 zz[1][2] = -r12;
00199 zz[2][2] = r20 + r21 + r23;
00200 zz[3][2] = -r32;
00201
00202
00203 zz[0][3] = -r03;
00204 zz[1][3] = -r13;
00205 zz[2][3] = -r23;
00206 zz[3][3] = r30 + r31 + r32;
00207
00208
00209 for( j=0; j <= 3; j++ )
00210 {
00211 for( i=0; i <= 3; i++ )
00212 {
00213 amat[i][j] = zz[i][j];
00214 }
00215 bvec[j] = zz[3+1][j];
00216 }
00217
00218 nerror = 0;
00219
00220
00221 getrf_wrapper(4, 4, (double*)amat, 4, ipiv, &nerror);
00222 getrs_wrapper('N', 4, 1, (double*)amat, 4, ipiv, bvec, 4, &nerror);
00223
00224
00225
00226
00227 if( nerror != 0 )
00228 {
00229 fprintf( ioQQQ, " AtomSeqBeryllium: dgetrs finds singular or ill-conditioned matrix\n" );
00230 puts( "[Stop in AtomSeqBeryllium]" );
00231 cdEXIT(EXIT_FAILURE);
00232 }
00233
00234
00235 for( i=0; i <= 3; i++ )
00236 {
00237 zz[3+1][i] = bvec[i];
00238 }
00239
00240 lgNegPop = false;
00241 for( i=0; i <= 3; i++ )
00242 {
00243 atoms.PopLevels[i] = zz[4][i]*AbunxIon;
00244 if( atoms.PopLevels[i] < 0. )
00245 lgNegPop = true;
00246 }
00247
00248
00249 if( lgNegPop )
00250 {
00251 fprintf( ioQQQ, " AtomSeqBeryllium finds non-positive pop,=" );
00252 for( i=0; i <= 3; i++ )
00253 {
00254 fprintf( ioQQQ, "%g ", atoms.PopLevels[i] );
00255 }
00256 fprintf( ioQQQ, "%s \n", chLab );
00257 fprintf( ioQQQ, " te=%g total abund=%g boltz=%g \n",
00258 phycon.te, AbunxIon, boltz );
00259 puts( "[Stop in AtomSeqBeryllium]" );
00260 cdEXIT(EXIT_FAILURE);
00261 }
00262
00263
00264 atoms.DepLTELevels[0] = 1.;
00265 atoms.DepLTELevels[1] = (atoms.PopLevels[1]/atoms.PopLevels[0])/
00266 excit;
00267 atoms.DepLTELevels[2] = (atoms.PopLevels[2]/atoms.PopLevels[0])/
00268 (excit*3.);
00269 atoms.DepLTELevels[3] = (atoms.PopLevels[3]/atoms.PopLevels[0])/
00270 (excit*5.);
00271
00272
00273 t->ColOvTot = (float)(c02/r02);
00274
00275 t->PopLo = AbunxIon;
00276
00277
00278 t->PopHi = atoms.PopLevels[2];
00279
00280 t->PopOpc = AbunxIon - atoms.PopLevels[2]*1./3.;
00281
00282
00283
00284 t->phots = atoms.PopLevels[2] * t->Aul * ( t->Pesc + t->Pelec_esc );
00285
00286 t->xIntensity = t->phots * t->EnergyErg;
00287
00288
00289
00290
00291
00292
00293
00294 Enr01 = atoms.PopLevels[0]*(c01 + c02 + c03)*t->EnergyErg;
00295 Enr10 = (atoms.PopLevels[1]*c10 + atoms.PopLevels[2]*c20 +
00296 atoms.PopLevels[3]*c30)*t->EnergyErg;
00297
00298
00299 t->cool = Enr01 - Enr10*t->ColOvTot;
00300
00301
00302 t->heat = Enr10*(1. - t->ColOvTot);
00303
00304
00305 coolng = t->cool;
00306 CoolAdd(chLab, t->WLAng ,coolng);
00307
00308
00309 thermal.dCooldT += coolng*(t->EnergyK*thermal.tsq1 - thermal.halfte);
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320 DEBUG_EXIT( "AtomSeqBeryllium()" );
00321 return;
00322 }
00323