00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "taulines.h"
00007 #include "doppvel.h"
00008 #include "iso.h"
00009 #include "trace.h"
00010 #include "dense.h"
00011 #include "rt.h"
00012 #include "rfield.h"
00013 #include "phycon.h"
00014 #include "lines_service.h"
00015 #include "thirdparty.h"
00016 #include "atoms.h"
00017
00018
00019 STATIC void oi_level_pops(double abundoi,
00020 double *coloi);
00021
00022
00023 void atom_oi_calc(double *coloi)
00024 {
00025 double esab,
00026 eslb,
00027 esoi,
00028 flb,
00029 foi,
00030 opaclb,
00031 opacoi,
00032 xlb,
00033 xoi;
00034 double aoi = TauLines[ipTO1025].Emis->Aul;
00035 double alb = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul;
00036
00037 DEBUG_ENTRY( "atom_oi_calc()" );
00038
00039 fixit();
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 if( dense.xIonDense[ipOXYGEN][0] <= 0. )
00051 {
00052 for( int i=0; i < 6; i++ )
00053 {
00054 atoms.popoi[i] = 0.;
00055 }
00056 *coloi = 0.;
00057 atoms.pmpo15 = 0.;
00058 atoms.pmpo51 = 0.;
00059 return;
00060 }
00061
00062
00063 esab = TauLines[ipTO1025].Emis->Pelec_esc + TauLines[ipTO1025].Emis->Pesc;
00064
00065
00066
00067 esoi = TauLines[ipTO1025].Emis->Pelec_esc + TauLines[ipTO1025].Emis->Pesc;
00068 eslb = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc +
00069 Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc;
00070
00071
00072 if( trace.lgTr8446 && trace.lgTrace )
00073 {
00074 fprintf( ioQQQ,
00075 " P8446 finds Lbeta, OI widths=%10.2e%10.2e and esc prob=%10.2e%10.2e esAB=%10.2e\n",
00076 GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]), GetDopplerWidth(dense.AtomicWeight[ipOXYGEN]), eslb, esoi, esab );
00077 }
00078
00079
00080 opacoi = 2.92e-9*dense.xIonDense[ipOXYGEN][0]*0.5556/GetDopplerWidth(dense.AtomicWeight[ipOXYGEN]);
00081 opaclb = 1.22e-8*StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop/GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]);
00082
00083
00084 xoi = opacoi/(opacoi + opaclb);
00085 xlb = opaclb/(opacoi + opaclb);
00086
00087
00088 foi = MIN2(GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]),GetDopplerWidth(dense.AtomicWeight[ipOXYGEN]))/GetDopplerWidth(dense.AtomicWeight[ipOXYGEN]);
00089 flb = MIN2(GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]),GetDopplerWidth(dense.AtomicWeight[ipOXYGEN]))/GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN])*
00090 MAX2(0.,1.- TauLines[ipTO1025].Emis->Pelec_esc - TauLines[ipTO1025].Emis->Pesc);
00091
00092 if( trace.lgTr8446 && trace.lgTrace )
00093 {
00094 fprintf( ioQQQ,
00095 " P8446 opac Lb, OI=%10.2e%10.2e X Lb, OI=%10.2e%10.2e FLb, OI=%10.2e%10.2e\n",
00096 opaclb, opacoi, xlb, xoi, flb, foi );
00097 }
00098
00099
00100
00101
00102 if( rfield.lgInducProcess )
00103 {
00104 atoms.pmpo15 = (realnum)((flb*alb*StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH3p].Pop*
00105 xoi*(1. - esab)/dense.xIonDense[ipOXYGEN][0]));
00106
00107 atoms.pmpo51 = (realnum)(aoi*(1. - (1. - foi)*(1. - esoi) - xoi*(1. - esab)*foi));
00108 }
00109 else
00110 {
00111 atoms.pmpo15 = 0.;
00112 atoms.pmpo51 = 0.;
00113 }
00114
00115
00116 oi_level_pops(dense.xIonDense[ipOXYGEN][0],coloi);
00117
00118
00119
00120
00121
00122 TauLines[ipT1304].Lo->Pop = atoms.popoi[0];
00123 TauLines[ipTO1025].Lo->Pop = atoms.popoi[0];
00124 TauLines[ipT1039].Lo->Pop = atoms.popoi[0];
00125 TauLines[ipT8446].Lo->Pop = atoms.popoi[1];
00126 TauLines[ipT4368].Lo->Pop = atoms.popoi[1];
00127 TauLines[ipTOI13].Lo->Pop = atoms.popoi[2];
00128 TauLines[ipTOI11].Lo->Pop = atoms.popoi[2];
00129 TauLines[ipTOI29].Lo->Pop = atoms.popoi[3];
00130 TauLines[ipTOI46].Lo->Pop = atoms.popoi[4];
00131
00132
00133 TauLines[ipT1304].Hi->Pop = atoms.popoi[1];
00134 TauLines[ipTO1025].Hi->Pop = atoms.popoi[4];
00135 TauLines[ipT1039].Hi->Pop = atoms.popoi[3];
00136 TauLines[ipT8446].Hi->Pop = atoms.popoi[2];
00137 TauLines[ipT4368].Hi->Pop = atoms.popoi[5];
00138 TauLines[ipTOI13].Hi->Pop = atoms.popoi[3];
00139 TauLines[ipTOI11].Hi->Pop = atoms.popoi[4];
00140 TauLines[ipTOI29].Hi->Pop = atoms.popoi[5];
00141 TauLines[ipTOI46].Hi->Pop = atoms.popoi[5];
00142
00143
00144
00146 TauLines[ipT1304].Emis->PopOpc = atoms.popoi[0];
00147 TauLines[ipTO1025].Emis->PopOpc = (atoms.popoi[0] - atoms.popoi[4]*0.6);
00148 TauLines[ipT1039].Emis->PopOpc = (atoms.popoi[0] - atoms.popoi[3]*3.0);
00149
00150
00152 TauLines[ipT8446].Emis->PopOpc =
00153 (MAX2(0.,atoms.popoi[1]-atoms.popoi[2]*.33));
00154 TauLines[ipT4368].Emis->PopOpc = (atoms.popoi[1] - atoms.popoi[5]*.33);
00155 TauLines[ipTOI13].Emis->PopOpc = (atoms.popoi[2] - atoms.popoi[3]*3.0);
00156 TauLines[ipTOI11].Emis->PopOpc = (atoms.popoi[2] - atoms.popoi[4]*0.6);
00157 TauLines[ipTOI29].Emis->PopOpc = (atoms.popoi[3] - atoms.popoi[5]*.33);
00158 TauLines[ipTOI46].Emis->PopOpc = (atoms.popoi[4] - atoms.popoi[5]*1.67);
00159 return;
00160 }
00161
00162
00163 STATIC void oi_level_pops(double abundoi,
00164 double *coloi)
00165 {
00166 bool lgNegPop;
00167
00168 long int i, j;
00169
00170 int32 ipiv[6], ner;
00171
00172 double a21,
00173 a32,
00174 a41,
00175 a43,
00176 a51,
00177 a53,
00178 a62,
00179 a64,
00180 a65,
00181 c12,
00182 c13,
00183 c14,
00184 c15,
00185 c16,
00186 c21,
00187 c23,
00188 c24,
00189 c25,
00190 c26,
00191 c31,
00192 c32,
00193 c34,
00194 c35,
00195 c36,
00196 c41,
00197 c42,
00198 c43,
00199 c45,
00200 c46,
00201 c51,
00202 c52,
00203 c53,
00204 c54,
00205 c56,
00206 c61,
00207 c62,
00208 c63,
00209 c64,
00210 c65,
00211 cs,
00212 deptoi[6],
00213 e12,
00214 e23,
00215 e34,
00216 e45,
00217 e56,
00218 simple;
00219
00220 double amat[6][6],
00221 bvec[6],
00222 zz[7][7];
00223
00224 static double g[6]={9.,3.,9.,3.,15.,9};
00225
00226 DEBUG_ENTRY( "oilevl()" );
00227
00228
00229
00230
00231
00232
00233
00234 e12 = rfield.ContBoltz[atoms.ipoiex[0]-1];
00235 e23 = rfield.ContBoltz[atoms.ipoiex[1]-1];
00236 e34 = rfield.ContBoltz[atoms.ipoiex[2]-1];
00237 e45 = rfield.ContBoltz[atoms.ipoiex[3]-1];
00238 e56 = rfield.ContBoltz[atoms.ipoiex[4]-1];
00239
00240
00241 a21 = TauLines[ipT1304].Emis->Aul*(TauLines[ipT1304].Emis->Pdest+ TauLines[ipT1304].Emis->Pesc + TauLines[ipT1304].Emis->Pelec_esc);
00242 a41 = TauLines[ipT1039].Emis->Aul*(TauLines[ipT1039].Emis->Pdest+ TauLines[ipT1039].Emis->Pesc + TauLines[ipT1039].Emis->Pelec_esc);
00243 a51 = TauLines[ipTO1025].Emis->Aul*(TauLines[ipTO1025].Emis->Pdest+ TauLines[ipTO1025].Emis->Pesc + TauLines[ipTO1025].Emis->Pelec_esc);
00244 a51 = atoms.pmpo51;
00245 a32 = TauLines[ipT8446].Emis->Aul*(TauLines[ipT8446].Emis->Pdest+ TauLines[ipT8446].Emis->Pesc + TauLines[ipT8446].Emis->Pelec_esc);
00246 a62 = TauLines[ipT4368].Emis->Aul*(TauLines[ipT4368].Emis->Pdest+ TauLines[ipT4368].Emis->Pesc + TauLines[ipT4368].Emis->Pelec_esc);
00247 a43 = TauLines[ipTOI13].Emis->Aul*(TauLines[ipTOI13].Emis->Pdest+ TauLines[ipTOI13].Emis->Pesc + TauLines[ipTOI13].Emis->Pelec_esc);
00248 a53 = TauLines[ipTOI11].Emis->Aul*(TauLines[ipTOI11].Emis->Pdest+ TauLines[ipTOI11].Emis->Pesc + TauLines[ipTOI11].Emis->Pelec_esc);
00249 a64 = TauLines[ipTOI29].Emis->Aul*(TauLines[ipTOI29].Emis->Pdest+ TauLines[ipTOI29].Emis->Pesc + TauLines[ipTOI29].Emis->Pelec_esc);
00250 a65 = TauLines[ipTOI46].Emis->Aul*(TauLines[ipTOI46].Emis->Pdest+ TauLines[ipTOI46].Emis->Pesc + TauLines[ipTOI46].Emis->Pelec_esc);
00251
00252
00253
00254
00255
00256 cs = 2.151e-5*phycon.te/phycon.te03;
00257 PutCS(cs,&TauLines[ipT1304]);
00258
00259
00260 cs = 9.25e-7*phycon.te*phycon.te10/phycon.te01/phycon.te01;
00261 PutCS(cs,&TauLines[ipTO1025]);
00262 c21 = dense.cdsqte*TauLines[ipT1304].Coll.col_str/g[1];
00263 c51 = dense.cdsqte*TauLines[ipTO1025].Coll.col_str/g[4];
00264
00265
00266 c31 = dense.cdsqte*1.0/g[2];
00267 PutCS(0.27,&TauLines[ipT1039]);
00268 c41 = dense.cdsqte*TauLines[ipT1039].Coll.col_str/g[3];
00269 c61 = dense.cdsqte*1./g[5];
00270
00271 c12 = c21*g[1]/g[0]*e12;
00272 c13 = c31*g[2]/g[0]*e12*e23;
00273 c14 = c41*g[3]/g[0]*e12*e23*e34;
00274 c15 = c51*g[4]/g[0]*e12*e23*e34*e45;
00275 c16 = c61*g[5]/g[0]*e12*e23*e34*e45*e56;
00276
00277 c32 = dense.cdsqte*85./g[2];
00278 c42 = dense.cdsqte*85./g[3];
00279 c52 = dense.cdsqte*85./g[4];
00280 c62 = dense.cdsqte*85./g[5];
00281
00282 c23 = c32*g[2]/g[1]*e23;
00283 c24 = c42*g[3]/g[1]*e23*e34;
00284 c25 = c52*g[4]/g[1]*e23*e34*e45;
00285 c26 = c62*g[5]/g[1]*e23*e34*e45*e56;
00286
00287 c43 = dense.cdsqte*70./g[3];
00288 c53 = dense.cdsqte*312./g[4];
00289 c63 = dense.cdsqte*1./g[5];
00290
00291 c34 = c43*g[3]/g[2]*e34;
00292 c35 = c53*g[4]/g[2]*e34*e45;
00293 c36 = c63*g[5]/g[2]*e34*e45*e56;
00294
00295 c54 = dense.cdsqte*50./g[4];
00296 c64 = dense.cdsqte*415./g[5];
00297
00298 c45 = c54*g[4]/g[3]*e45;
00299 c46 = c64*g[5]/g[3]*e45*e56;
00300
00301 c65 = dense.cdsqte*400./g[5];
00302 c56 = c65*g[5]/g[4]*e56;
00303
00306
00307 simple = (c16 + atoms.pmpo15)/(c61 + c62 + c64 + a65 + a64 + a62);
00308 if( simple < 1e-19 )
00309 {
00310 atoms.popoi[0] = abundoi;
00311 for( i=1; i < 6; i++ )
00312 {
00313 atoms.popoi[i] = 0.;
00314 }
00315 *coloi = 0.;
00316 return;
00317 }
00318
00319
00320
00321 for( i=0; i < 6; i++ )
00322 {
00323 zz[i][0] = 1.0;
00324 zz[6][i] = 0.;
00325 }
00326
00327
00328 zz[6][0] = abundoi;
00329
00330
00331 zz[0][1] = -c12;
00332 zz[1][1] = c21 + c23 + c24 + c25 + c26 + a21;
00333 zz[2][1] = -c32 - a32;
00334 zz[3][1] = -c42;
00335 zz[4][1] = -c52;
00336 zz[5][1] = -c62 - a62;
00337
00338
00339 zz[0][2] = -c13;
00340 zz[1][2] = -c23;
00341 zz[2][2] = c31 + c32 + c34 + c35 + c36 + a32;
00342 zz[3][2] = -c43 - a43;
00343 zz[4][2] = -c53 - a53;
00344 zz[5][2] = -c63;
00345
00346
00347 zz[0][3] = -c14;
00348 zz[1][3] = -c24;
00349 zz[2][3] = -c34;
00350 zz[3][3] = c41 + c42 + c43 + c45 + c46 + a41 + a43;
00351 zz[4][3] = -c54;
00352 zz[5][3] = -c64 - a64;
00353
00354
00355 zz[0][4] = -c15 - atoms.pmpo15;
00356 zz[1][4] = -c25;
00357 zz[2][4] = -c35;
00358 zz[3][4] = -c45;
00359 zz[4][4] = c51 + c52 + c53 + c54 + c56 + a51 + a53;
00360 zz[5][4] = -c65 - a65;
00361
00362
00363 zz[0][5] = -c16;
00364 zz[1][5] = -c26;
00365 zz[2][5] = -c36;
00366 zz[3][5] = -c46;
00367 zz[4][5] = -c56;
00368 zz[5][5] = c61 + c62 + c63 + c64 + c65 + a65 + a64 + a62;
00369
00370
00371 for( j=0; j < 6; j++ )
00372 {
00373 for( i=0; i < 6; i++ )
00374 {
00375 amat[i][j] = zz[i][j];
00376 }
00377 bvec[j] = zz[6][j];
00378 }
00379
00380 ner = 0;
00381
00382 getrf_wrapper(6, 6, (double*)amat, 6, ipiv, &ner);
00383 getrs_wrapper('N', 6, 1, (double*)amat, 6, ipiv, bvec, 6, &ner);
00384
00385
00386
00387 if( ner != 0 )
00388 {
00389 fprintf( ioQQQ, " oi_level_pops: dgetrs finds singular or ill-conditioned matrix\n" );
00390 cdEXIT(EXIT_FAILURE);
00391 }
00392
00393
00394
00395 for( i=0; i < 6; i++ )
00396 {
00397 zz[6][i] = bvec[i];
00398 }
00399
00400 lgNegPop = false;
00401 for( i=0; i < 6; i++ )
00402 {
00403 atoms.popoi[i] = zz[6][i];
00404 if( atoms.popoi[i] < 0. )
00405 lgNegPop = true;
00406 }
00407
00408
00409
00410 if( trace.lgTrace && trace.lgTr8446 )
00411 {
00412 deptoi[0] = 1.;
00413 deptoi[1] = atoms.popoi[1]/atoms.popoi[0]/(g[1]/g[0]*
00414 e12);
00415 deptoi[2] = atoms.popoi[2]/atoms.popoi[0]/(g[2]/g[0]*
00416 e12*e23);
00417 deptoi[3] = atoms.popoi[3]/atoms.popoi[0]/(g[3]/g[0]*
00418 e12*e23*e34);
00419 deptoi[4] = atoms.popoi[4]/atoms.popoi[0]/(g[4]/g[0]*
00420 e12*e23*e34*e45);
00421 deptoi[5] = atoms.popoi[5]/atoms.popoi[0]/(g[5]/g[0]*
00422 e12*e23*e34*e45*e56);
00423
00424 fprintf( ioQQQ, " oilevl finds levl pop" );
00425 for(i=0; i < 6; i++)
00426 fprintf( ioQQQ, "%11.3e", atoms.popoi[i] );
00427 fprintf( ioQQQ, "\n" );
00428
00429 fprintf( ioQQQ, " oilevl finds dep coef" );
00430 for(i=0; i < 6; i++)
00431 fprintf( ioQQQ, "%11.3e", deptoi[i] );
00432 fprintf( ioQQQ, "\n" );
00433 }
00434
00435
00436 if( lgNegPop )
00437 {
00438 atoms.nNegOI += 1;
00439
00440 fprintf( ioQQQ, " OILEVL finds negative population" );
00441 for(i=0;i < 6; i++)
00442 fprintf( ioQQQ, "%10.2e", atoms.popoi[i] );
00443 fprintf( ioQQQ, "\n" );
00444
00445 fprintf( ioQQQ, " simple 5 =%10.2e\n", simple );
00446
00447 atoms.popoi[5] = 0.;
00448 atoms.popoi[4] = abundoi*(c15 + atoms.pmpo15)/(a51 + a53 + c51 + c53);
00449 atoms.popoi[3] = 0.;
00450 atoms.popoi[2] = (atoms.popoi[4]*(a53 + c53)) / (a32 + c32);
00451 atoms.popoi[1] = (atoms.popoi[4]*(a53 + c53) + abundoi*c12) / (a21 + c21);
00452 atoms.popoi[0] = abundoi;
00453
00454
00455 }
00456
00457
00458 *coloi =
00459 (atoms.popoi[0]*c12 - atoms.popoi[1]*c21)*1.53e-11 +
00460 (atoms.popoi[0]*c14 - atoms.popoi[3]*c41)*1.92e-11 +
00461 (atoms.popoi[0]*c15 - atoms.popoi[4]*c51)*1.94e-11 +
00462 (atoms.popoi[1]*c23 - atoms.popoi[2]*c32)*2.36e-12 +
00463 (atoms.popoi[1]*c26 - atoms.popoi[5]*c62)*4.55e-12 +
00464 (atoms.popoi[2]*c35 - atoms.popoi[4]*c53)*1.76e-12 +
00465 (atoms.popoi[2]*c34 - atoms.popoi[3]*c43)*1.52e-12 +
00466 (atoms.popoi[3]*c46 - atoms.popoi[5]*c64)*6.86e-13 +
00467 (atoms.popoi[4]*c56 - atoms.popoi[5]*c65)*4.32e-13;
00468
00469 return;
00470 }