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