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 = iso_sp[ipH_LIKE][ipHYDROGEN].trans(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 = iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH3p,ipH1s).Emis().Pelec_esc() +
00069 iso_sp[ipH_LIKE][ipHYDROGEN].trans(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*iso_sp[ipH_LIKE][ipHYDROGEN].st[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*iso_sp[ipH_LIKE][ipHYDROGEN].st[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 TauLines[ipT1304].Emis().PopOpc() = atoms.popoi[0] - atoms.popoi[1]*3.0;
00144 TauLines[ipTO1025].Emis().PopOpc() = (atoms.popoi[0] - atoms.popoi[4]*0.6);
00145 TauLines[ipT1039].Emis().PopOpc() = (atoms.popoi[0] - atoms.popoi[3]*3.0);
00146
00147
00149 TauLines[ipT8446].Emis().PopOpc() =
00150 (MAX2(0.,atoms.popoi[1]-atoms.popoi[2]*.33));
00151 TauLines[ipT4368].Emis().PopOpc() = (atoms.popoi[1] - atoms.popoi[5]*.33);
00152 TauLines[ipTOI13].Emis().PopOpc() = (atoms.popoi[2] - atoms.popoi[3]*3.0);
00153 TauLines[ipTOI11].Emis().PopOpc() = (atoms.popoi[2] - atoms.popoi[4]*0.6);
00154 TauLines[ipTOI29].Emis().PopOpc() = (atoms.popoi[3] - atoms.popoi[5]*.33);
00155 TauLines[ipTOI46].Emis().PopOpc() = (atoms.popoi[4] - atoms.popoi[5]*1.67);
00156 return;
00157 }
00158
00159
00160 STATIC void oi_level_pops(double abundoi,
00161 double *coloi)
00162 {
00163 bool lgNegPop;
00164
00165 long int i, j;
00166
00167 int32 ipiv[6], ner;
00168
00169 double a21,
00170 a32,
00171 a41,
00172 a43,
00173 a51,
00174 a53,
00175 a62,
00176 a64,
00177 a65,
00178 c12,
00179 c13,
00180 c14,
00181 c15,
00182 c16,
00183 c21,
00184 c23,
00185 c24,
00186 c25,
00187 c26,
00188 c31,
00189 c32,
00190 c34,
00191 c35,
00192 c36,
00193 c41,
00194 c42,
00195 c43,
00196 c45,
00197 c46,
00198 c51,
00199 c52,
00200 c53,
00201 c54,
00202 c56,
00203 c61,
00204 c62,
00205 c63,
00206 c64,
00207 c65,
00208 cs,
00209 deptoi[6],
00210 e12,
00211 e23,
00212 e34,
00213 e45,
00214 e56,
00215 simple;
00216
00217 double amat[6][6],
00218 bvec[6],
00219 zz[7][7];
00220
00221 static double g[6]={9.,3.,9.,3.,15.,9};
00222
00223 DEBUG_ENTRY( "oilevl()" );
00224
00225
00226
00227
00228
00229
00230
00231 e12 = rfield.ContBoltz[atoms.ipoiex[0]-1];
00232 e23 = rfield.ContBoltz[atoms.ipoiex[1]-1];
00233 e34 = rfield.ContBoltz[atoms.ipoiex[2]-1];
00234 e45 = rfield.ContBoltz[atoms.ipoiex[3]-1];
00235 e56 = rfield.ContBoltz[atoms.ipoiex[4]-1];
00236
00237
00238 a21 = TauLines[ipT1304].Emis().Aul()*(TauLines[ipT1304].Emis().Pdest()+ TauLines[ipT1304].Emis().Pesc() + TauLines[ipT1304].Emis().Pelec_esc());
00239 a41 = TauLines[ipT1039].Emis().Aul()*(TauLines[ipT1039].Emis().Pdest()+ TauLines[ipT1039].Emis().Pesc() + TauLines[ipT1039].Emis().Pelec_esc());
00240 a51 = TauLines[ipTO1025].Emis().Aul()*(TauLines[ipTO1025].Emis().Pdest()+ TauLines[ipTO1025].Emis().Pesc() + TauLines[ipTO1025].Emis().Pelec_esc());
00241 a51 = atoms.pmpo51;
00242 a32 = TauLines[ipT8446].Emis().Aul()*(TauLines[ipT8446].Emis().Pdest()+ TauLines[ipT8446].Emis().Pesc() + TauLines[ipT8446].Emis().Pelec_esc());
00243 a62 = TauLines[ipT4368].Emis().Aul()*(TauLines[ipT4368].Emis().Pdest()+ TauLines[ipT4368].Emis().Pesc() + TauLines[ipT4368].Emis().Pelec_esc());
00244 a43 = TauLines[ipTOI13].Emis().Aul()*(TauLines[ipTOI13].Emis().Pdest()+ TauLines[ipTOI13].Emis().Pesc() + TauLines[ipTOI13].Emis().Pelec_esc());
00245 a53 = TauLines[ipTOI11].Emis().Aul()*(TauLines[ipTOI11].Emis().Pdest()+ TauLines[ipTOI11].Emis().Pesc() + TauLines[ipTOI11].Emis().Pelec_esc());
00246 a64 = TauLines[ipTOI29].Emis().Aul()*(TauLines[ipTOI29].Emis().Pdest()+ TauLines[ipTOI29].Emis().Pesc() + TauLines[ipTOI29].Emis().Pelec_esc());
00247 a65 = TauLines[ipTOI46].Emis().Aul()*(TauLines[ipTOI46].Emis().Pdest()+ TauLines[ipTOI46].Emis().Pesc() + TauLines[ipTOI46].Emis().Pelec_esc());
00248
00249
00250
00251
00252
00253 cs = 2.151e-5*phycon.te/phycon.te03;
00254 PutCS(cs,TauLines[ipT1304]);
00255
00256
00257 cs = 9.25e-7*phycon.te*phycon.te10/phycon.te01/phycon.te01;
00258 PutCS(cs,TauLines[ipTO1025]);
00259 c21 = dense.cdsqte*TauLines[ipT1304].Coll().col_str()/g[1];
00260 c51 = dense.cdsqte*TauLines[ipTO1025].Coll().col_str()/g[4];
00261
00262
00263 c31 = dense.cdsqte*1.0/g[2];
00264 PutCS(0.27,TauLines[ipT1039]);
00265 c41 = dense.cdsqte*TauLines[ipT1039].Coll().col_str()/g[3];
00266 c61 = dense.cdsqte*1./g[5];
00267
00268 c12 = c21*g[1]/g[0]*e12;
00269 c13 = c31*g[2]/g[0]*e12*e23;
00270 c14 = c41*g[3]/g[0]*e12*e23*e34;
00271 c15 = c51*g[4]/g[0]*e12*e23*e34*e45;
00272 c16 = c61*g[5]/g[0]*e12*e23*e34*e45*e56;
00273
00274 c32 = dense.cdsqte*85./g[2];
00275 c42 = dense.cdsqte*85./g[3];
00276 c52 = dense.cdsqte*85./g[4];
00277 c62 = dense.cdsqte*85./g[5];
00278
00279 c23 = c32*g[2]/g[1]*e23;
00280 c24 = c42*g[3]/g[1]*e23*e34;
00281 c25 = c52*g[4]/g[1]*e23*e34*e45;
00282 c26 = c62*g[5]/g[1]*e23*e34*e45*e56;
00283
00284 c43 = dense.cdsqte*70./g[3];
00285 c53 = dense.cdsqte*312./g[4];
00286 c63 = dense.cdsqte*1./g[5];
00287
00288 c34 = c43*g[3]/g[2]*e34;
00289 c35 = c53*g[4]/g[2]*e34*e45;
00290 c36 = c63*g[5]/g[2]*e34*e45*e56;
00291
00292 c54 = dense.cdsqte*50./g[4];
00293 c64 = dense.cdsqte*415./g[5];
00294
00295 c45 = c54*g[4]/g[3]*e45;
00296 c46 = c64*g[5]/g[3]*e45*e56;
00297
00298 c65 = dense.cdsqte*400./g[5];
00299 c56 = c65*g[5]/g[4]*e56;
00300
00303
00304 simple = (c16 + atoms.pmpo15)/(c61 + c62 + c64 + a65 + a64 + a62);
00305 if( simple < 1e-19 )
00306 {
00307 atoms.popoi[0] = abundoi;
00308 for( i=1; i < 6; i++ )
00309 {
00310 atoms.popoi[i] = 0.;
00311 }
00312 *coloi = 0.;
00313 return;
00314 }
00315
00316
00317
00318 for( i=0; i < 6; i++ )
00319 {
00320 zz[i][0] = 1.0;
00321 zz[6][i] = 0.;
00322 }
00323
00324
00325 zz[6][0] = abundoi;
00326
00327
00328 zz[0][1] = -c12;
00329 zz[1][1] = c21 + c23 + c24 + c25 + c26 + a21;
00330 zz[2][1] = -c32 - a32;
00331 zz[3][1] = -c42;
00332 zz[4][1] = -c52;
00333 zz[5][1] = -c62 - a62;
00334
00335
00336 zz[0][2] = -c13;
00337 zz[1][2] = -c23;
00338 zz[2][2] = c31 + c32 + c34 + c35 + c36 + a32;
00339 zz[3][2] = -c43 - a43;
00340 zz[4][2] = -c53 - a53;
00341 zz[5][2] = -c63;
00342
00343
00344 zz[0][3] = -c14;
00345 zz[1][3] = -c24;
00346 zz[2][3] = -c34;
00347 zz[3][3] = c41 + c42 + c43 + c45 + c46 + a41 + a43;
00348 zz[4][3] = -c54;
00349 zz[5][3] = -c64 - a64;
00350
00351
00352 zz[0][4] = -c15 - atoms.pmpo15;
00353 zz[1][4] = -c25;
00354 zz[2][4] = -c35;
00355 zz[3][4] = -c45;
00356 zz[4][4] = c51 + c52 + c53 + c54 + c56 + a51 + a53;
00357 zz[5][4] = -c65 - a65;
00358
00359
00360 zz[0][5] = -c16;
00361 zz[1][5] = -c26;
00362 zz[2][5] = -c36;
00363 zz[3][5] = -c46;
00364 zz[4][5] = -c56;
00365 zz[5][5] = c61 + c62 + c63 + c64 + c65 + a65 + a64 + a62;
00366
00367
00368 for( j=0; j < 6; j++ )
00369 {
00370 for( i=0; i < 6; i++ )
00371 {
00372 amat[i][j] = zz[i][j];
00373 }
00374 bvec[j] = zz[6][j];
00375 }
00376
00377 ner = 0;
00378
00379 getrf_wrapper(6, 6, (double*)amat, 6, ipiv, &ner);
00380 getrs_wrapper('N', 6, 1, (double*)amat, 6, ipiv, bvec, 6, &ner);
00381
00382
00383
00384 if( ner != 0 )
00385 {
00386 fprintf( ioQQQ, " oi_level_pops: dgetrs finds singular or ill-conditioned matrix\n" );
00387 cdEXIT(EXIT_FAILURE);
00388 }
00389
00390
00391
00392 for( i=0; i < 6; i++ )
00393 {
00394 zz[6][i] = bvec[i];
00395 }
00396
00397 lgNegPop = false;
00398 for( i=0; i < 6; i++ )
00399 {
00400 atoms.popoi[i] = zz[6][i];
00401 if( atoms.popoi[i] < 0. )
00402 lgNegPop = true;
00403 }
00404
00405
00406
00407 if( trace.lgTrace && trace.lgTr8446 )
00408 {
00409 deptoi[0] = 1.;
00410 deptoi[1] = atoms.popoi[1]/atoms.popoi[0]/(g[1]/g[0]*
00411 e12);
00412 deptoi[2] = atoms.popoi[2]/atoms.popoi[0]/(g[2]/g[0]*
00413 e12*e23);
00414 deptoi[3] = atoms.popoi[3]/atoms.popoi[0]/(g[3]/g[0]*
00415 e12*e23*e34);
00416 deptoi[4] = atoms.popoi[4]/atoms.popoi[0]/(g[4]/g[0]*
00417 e12*e23*e34*e45);
00418 deptoi[5] = atoms.popoi[5]/atoms.popoi[0]/(g[5]/g[0]*
00419 e12*e23*e34*e45*e56);
00420
00421 fprintf( ioQQQ, " oilevl finds levl pop" );
00422 for(i=0; i < 6; i++)
00423 fprintf( ioQQQ, "%11.3e", atoms.popoi[i] );
00424 fprintf( ioQQQ, "\n" );
00425
00426 fprintf( ioQQQ, " oilevl finds dep coef" );
00427 for(i=0; i < 6; i++)
00428 fprintf( ioQQQ, "%11.3e", deptoi[i] );
00429 fprintf( ioQQQ, "\n" );
00430 }
00431
00432
00433 if( lgNegPop )
00434 {
00435 atoms.nNegOI += 1;
00436
00437 fprintf( ioQQQ, " OILEVL finds negative population" );
00438 for(i=0;i < 6; i++)
00439 fprintf( ioQQQ, "%10.2e", atoms.popoi[i] );
00440 fprintf( ioQQQ, "\n" );
00441
00442 fprintf( ioQQQ, " simple 5 =%10.2e\n", simple );
00443
00444 atoms.popoi[5] = 0.;
00445 atoms.popoi[4] = abundoi*(c15 + atoms.pmpo15)/(a51 + a53 + c51 + c53);
00446 atoms.popoi[3] = 0.;
00447 atoms.popoi[2] = (atoms.popoi[4]*(a53 + c53)) / (a32 + c32);
00448 atoms.popoi[1] = (atoms.popoi[4]*(a53 + c53) + abundoi*c12) / (a21 + c21);
00449 atoms.popoi[0] = abundoi;
00450
00451
00452 }
00453
00454
00455 *coloi =
00456 (atoms.popoi[0]*c12 - atoms.popoi[1]*c21)*1.53e-11 +
00457 (atoms.popoi[0]*c14 - atoms.popoi[3]*c41)*1.92e-11 +
00458 (atoms.popoi[0]*c15 - atoms.popoi[4]*c51)*1.94e-11 +
00459 (atoms.popoi[1]*c23 - atoms.popoi[2]*c32)*2.36e-12 +
00460 (atoms.popoi[1]*c26 - atoms.popoi[5]*c62)*4.55e-12 +
00461 (atoms.popoi[2]*c35 - atoms.popoi[4]*c53)*1.76e-12 +
00462 (atoms.popoi[2]*c34 - atoms.popoi[3]*c43)*1.52e-12 +
00463 (atoms.popoi[3]*c46 - atoms.popoi[5]*c64)*6.86e-13 +
00464 (atoms.popoi[4]*c56 - atoms.popoi[5]*c65)*4.32e-13;
00465
00466 return;
00467 }