00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "embesq.h"
00007 #include "phycon.h"
00008 #include "taulines.h"
00009 #include "dense.h"
00010 #include "hmi.h"
00011 #include "h2.h"
00012 #include "co.h"
00013 #include "ligbar.h"
00014 #include "mole.h"
00015 #include "thermal.h"
00016 #include "colden.h"
00017 #include "lines_service.h"
00018 #include "atoms.h"
00019 #include "carb.h"
00020 #include "cooling.h"
00021
00022 void CoolCarb(void)
00023 {
00024 double SaveAbun,
00025 a21,
00026 a31,
00027 a32,
00028 cs,
00029 cs01,
00030 cs02,
00031 cs12,
00032 cs13,
00033 cs23,
00034 cs2s2p,
00035 cs2s3p ,
00036 ortho_frac ,
00037 popup,
00038 popratio;
00039
00040
00041
00042 double cse01,
00043 cse12,
00044 cse02,
00045 csh01,
00046 csh12,
00047 csh02,
00048 csp01,
00049 csp12,
00050 csp02,
00051 csh201,
00052 csh212,
00053 csh202 ,
00054 csh2p01,
00055 csh2p12,
00056 csh2p02,
00057 csh2o01,
00058 csh2o12,
00059 csh2o02,
00060 temp;
00061 double cs_c2_h12=-1.;
00062 realnum pciexc;
00063 int i;
00064
00065 DEBUG_ENTRY( "CoolCarb()" );
00066
00067 (*TauDummy).Zero();
00068 (*(*TauDummy).Hi()).g() = 0.;
00069 (*(*TauDummy).Lo()).g() = 0.;
00070 (*(*TauDummy).Hi()).IonStg() = 0;
00071 (*(*TauDummy).Lo()).IonStg() = 0;
00072 (*(*TauDummy).Hi()).nelem() = 0;
00073 (*(*TauDummy).Lo()).nelem() = 0;
00074 (*TauDummy).Emis().Aul() = 0.;
00075 (*TauDummy).EnergyWN() = 0.;
00076 (*TauDummy).WLAng() = 0.;
00077
00078
00079
00080
00081
00082
00083
00084
00085 PutCS(7.3, TauLines[ipT1656] );
00086 atom_level2(TauLines[ipT1656]);
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096 if( phycon.te<=3.0e3 )
00097 {
00098
00099
00100 cse01 = MAX2(4.80E-06*phycon.te*phycon.te20/phycon.te03,
00101 8.24E-07*phycon.te32/phycon.te01);
00102
00103 cse12 = MAX2(7.67E-05*phycon.te/phycon.te10/phycon.te03,
00104 1.47E-06*phycon.te32*phycon.te10/phycon.te03);
00105
00106 cse02 = MAX2(4.72E-05*phycon.te70*phycon.te03,
00107 3.05E-07*phycon.te32*phycon.te10);
00108 }
00109 else
00110 {
00111
00112
00113 cse01 = MIN2(8.24E-07*phycon.te32/phycon.te01,
00114 0.0035*phycon.sqrte*phycon.te01);
00115
00116 cse12 = MIN2(1.47E-06*phycon.te32*phycon.te10/phycon.te03,
00117 0.0088*phycon.sqrte*phycon.te01*phycon.te005);
00118
00119 cse02 = MIN2(3.05E-07*phycon.te32*phycon.te10,
00120 0.00448*phycon.sqrte/phycon.te10*phycon.te03*phycon.te005);
00121 }
00122
00123
00124
00125
00126
00127
00128
00129
00130 csh01 = MAX2(1.61e-10,5.66e-11*phycon.te20);
00131
00132
00133 csh12 = MAX2(1.93e-10*phycon.te05*phycon.te03,
00134 5.64e-11*phycon.te30*phycon.te02);
00135
00136
00137 csh02 = MAX2(1.08e-10/phycon.te03,
00138 1.67e-11*phycon.te30*phycon.te02*phycon.te02);
00139
00140
00141
00142
00143
00144 if( phycon.te < 25. )
00145 temp = 25.;
00146 else if( phycon.te >20000. )
00147 temp = 20000.;
00148 else
00149 temp = phycon.te;
00150 csp01 = 1e-9*pow2(4.3671821 - 14.39018/log(temp))*(1./3.)*exp(16.4*T1CM/temp);
00151 csp12 = 1e-9*exp(3.2823932 - 60.99754*(log(temp)/temp))*(1./5.)*exp(37.1*T1CM/temp);
00152 csp02 = 1e-9/(0.033932579+ (1503.4042/pow(temp,1.5)))*(3./5.)*exp(43.5*T1CM/temp);
00153
00154
00155
00156
00157
00158
00159
00160
00161 ortho_frac = h2.ortho_density/SDIV(hmi.H2_total);
00162
00163
00164
00165
00166
00167 if( phycon.te<=30. )
00168 {
00169 csh2p01 = MIN2(8.38E-11*phycon.te05*phycon.te01,
00170 2.12e-10/phycon.te20/phycon.te05/phycon.te01);
00171
00172 csh2o01 = MIN2(5.17E-11*phycon.te10*phycon.te05,
00173 1.07e-10/phycon.te10*phycon.te01);
00174 }
00175 else if( phycon.te<=150. )
00176 {
00177 csh2p01 = MAX2(6.60e-11,
00178 2.12e-10/phycon.te20/phycon.te05/phycon.te01);
00179
00180 csh2o01 = MAX2(7.10e-11,
00181 1.07e-10/phycon.te10*phycon.te01);
00182 }
00183 else
00184 {
00185
00186
00187 csh2p01 = MAX2(6.60e-11,3.38e-11*phycon.te10*phycon.te03);
00188 csh2p01 = MIN2(8.10e-11,csh2p01);
00189
00190 csh2o01 = MAX2(7.1e-11,3.37e-11*phycon.te10*phycon.te02*phycon.te02);
00191 csh2o01 = MIN2(8.57e-11,csh2o01);
00192 }
00193
00194
00195 csh201 = ortho_frac*csh2o01 + (1.-ortho_frac)*csh2p01;
00196 if( phycon.te<=30. )
00197 {
00198 csh2p12 = MIN2(1.48E-10*phycon.te05*phycon.te02,
00199 2.25e-10/phycon.te03/phycon.te03);
00200 }
00201 else if( phycon.te <= 100. )
00202 {
00203 csh2p12 = MAX2(1.75e-10,
00204 2.25e-10/phycon.te03/phycon.te03);
00205 }
00206 else
00207 {
00208 csh2p12 = MAX2(1.75e-10,6.23e-11*phycon.te20*phycon.te01);
00209 csh2p12 = MIN2(2.61e-10,csh2p12);
00210 }
00211
00212 csh2o12 = MIN2(2.83e-10,4.46e-11*phycon.te30/phycon.te03);
00213 {
00214
00215 csh212 = ortho_frac*csh2o12 + (1.-ortho_frac)*csh2p12;
00216 }
00217
00218 if( phycon.te<=30 )
00219 {
00220 csh2p02 = MIN2(8.67E-11*phycon.te02*phycon.te02,
00221 1.35e-10/phycon.te10);
00222 }
00223 else if( phycon.te<=150. )
00224 {
00225 csh2p02 = MAX2(8.40e-11,
00226 1.35e-10/phycon.te10);
00227 }
00228 else
00229 {
00230 csh2p02 = MAX2(8.4e-11,4.04e-11*phycon.te10*phycon.te02*phycon.te02);
00231 csh2p02 = MIN2(1.04e-10,csh2p02);
00232 }
00233
00234 csh2o02 = MIN2(1.11e-10,3.16e-11*phycon.te20/phycon.te02);
00235
00236 csh202 = ortho_frac*csh2o02 + (1.-ortho_frac)*csh2p02;
00237
00240
00241
00242
00243
00244 cs01 = cse01 + 3.*(csh01*(dense.xIonDense[ipHYDROGEN][0]+dense.xIonDense[ipHELIUM][0]) + csp01*dense.xIonDense[ipHYDROGEN][1] + csh201*hmi.H2_total)/dense.cdsqte;
00245 cs12 = cse12 + 5.*(csh12*(dense.xIonDense[ipHYDROGEN][0]+dense.xIonDense[ipHELIUM][0]) + csp12*dense.xIonDense[ipHYDROGEN][1] + csh212*hmi.H2_total)/dense.cdsqte;
00246 cs02 = cse02 + 5.*(csh02*(dense.xIonDense[ipHYDROGEN][0]+dense.xIonDense[ipHELIUM][0]) + csp02*dense.xIonDense[ipHYDROGEN][1] + csh202*hmi.H2_total)/dense.cdsqte;
00247
00248 PutCS( cs01 , TauLines[ipT610] );
00249 PutCS( cs12 , TauLines[ipT370] );
00250 PutCS( cs02 , *TauDummy );
00251
00252
00253
00254 atom_level3(
00255 TauLines[ipT610],
00256 TauLines[ipT370],
00257 *TauDummy);
00258
00259
00260 for( i=0; i<3; ++i)
00261 {
00262
00263 colden.C1Pops[i] = (realnum)atoms.PopLevels[i];
00264 }
00265
00266
00267
00268
00269 if( dense.xIonDense[ipCARBON][0] > 0. && phycon.te < 40000. )
00270 {
00271 cs12 = 1.156e-4*phycon.te*(1.09 - 7.5e-6*phycon.te - 2.1e-10*
00272 phycon.te*phycon.te);
00273 cs13 = 2.8e-3*phycon.sqrte;
00274 cs23 = 2.764e-3*phycon.sqrte;
00275
00276 a21 = 3.26e-4*TauLines[ipT9830].Emis().Pesc();
00277 a31 = 2.73e-3;
00278 a32 = 0.528*TauLines[ipT8727].Emis().Pesc();
00280 carb.c8727 = atom_pop3(9.,5.,1.,cs12,cs13,cs23,a21,a31,a32,
00281 1.417e4,1.255e4,&pciexc,dense.xIonDense[ipCARBON][0],0.,0.,0.)*a32*
00282 2.28e-12;
00283 TauLines[ipT9830].Emis().PopOpc() = dense.xIonDense[ipCARBON][0];
00284 (*TauLines[ipT9830].Lo()).Pop() = dense.xIonDense[ipCARBON][0];
00285 (*TauLines[ipT9830].Hi()).Pop() = 0.;
00286 TauLines[ipT9830].Coll().col_str() = (realnum)cs12;
00287 TauLines[ipT8727].Emis().PopOpc() = (carb.c8727/(a32*2.28e-12));
00288 (*TauLines[ipT8727].Lo()).Pop() = (carb.c8727/(a32*2.28e-12));
00289 (*TauLines[ipT8727].Hi()).Pop() = 0.;
00290 TauLines[ipT8727].Coll().col_str() = (realnum)cs23;
00291
00292 carb.c9850 = pciexc*a21*2.02e-12;
00293 thermal.dCooldT += carb.c9850*(1.468e4*thermal.tsq1 + thermal.halfte);
00294 thermal.dCooldT += carb.c8727*(1.255e4*thermal.tsq1 + thermal.halfte);
00295
00296
00297 carb.r9850 = (realnum)(a21/(a21 + cs12/5.*COLL_CONST/phycon.sqrte*dense.eden));
00298 }
00299
00300 else
00301 {
00302 carb.c9850 = 0.;
00303 carb.c8727 = 0.;
00304 carb.r9850 = 0.;
00305 TauLines[ipT9830].Emis().PopOpc() = 0.;
00306 (*TauLines[ipT9830].Lo()).Pop() = 0.;
00307 (*TauLines[ipT9830].Hi()).Pop() = 0.;
00308 TauLines[ipT8727].Emis().PopOpc() = 0.;
00309 (*TauLines[ipT8727].Lo()).Pop() = 0.;
00310 (*TauLines[ipT8727].Hi()).Pop() = 0.;
00311 }
00312 CoolAdd("C 1",8727,carb.c8727);
00313 CoolAdd("C 1",9850,carb.c9850);
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329 cse12 = MIN2(2.20,0.403*phycon.te20/phycon.te02*phycon.te001*phycon.te001);
00330
00331
00332
00333
00334
00335
00336
00337 temp = MIN2(2e3, phycon.te);
00338
00339
00340
00341 cs_c2_h12 = 1e-10*(4.4716028+ 0.69658785*pow(temp, 0.31692387));
00342
00343 if(phycon.te > 2e3)
00344 {
00345
00346 cs_c2_h12 *= pow(phycon.te/2e3, 0.31692387);
00347 }
00348
00349
00350 cs_c2_h12 *= 4.*(dense.xIonDense[ipHYDROGEN][0] + findspecieslocal("H2")->den)/dense.cdsqte;
00351
00352
00353 ASSERT( fabs(dense.eden + dense.xIonDense[ipHYDROGEN][0]*1.7e-4 * dense.HCorrFac -dense.EdenHCorr )/
00354 dense.EdenHCorr < 1e-8 );
00355
00356 PutCS( cse12+cs_c2_h12 ,TauLines[ipT157]);
00357
00358
00359
00360
00361 cs = MIN2(6.73,2.316*phycon.te10);
00362 PutCS(cs,TauLines[ipT1335]);
00363 atom_level2(TauLines[ipT1335]);
00364
00365 static vector< pair<TransitionList::iterator,double> > C2Pump;
00366 C2Pump.reserve(32);
00367
00368
00369 if( C2Pump.empty() )
00370 {
00371
00372 pair<TransitionList::iterator,double> pp( TauLines.begin()+ipT1335, 1./6. );
00373 C2Pump.push_back( pp );
00374
00375 for( i=0; i < nWindLine; ++i )
00376 {
00377
00378 if( (*TauLine2[i].Hi()).nelem() == 6 && (*TauLine2[i].Hi()).IonStg() == 2 )
00379 {
00380 # if 0
00381 DumpLine( &TauLine2[i] );
00382 # endif
00383 double branch_ratio;
00384
00385
00386 if( fp_equal( (*TauLine2[i].Hi()).g(), realnum(2.) ) )
00387 branch_ratio = 2./3.;
00388 else if( fp_equal( (*TauLine2[i].Hi()).g(), realnum(6.) ) )
00389 branch_ratio = 1./2.;
00390 else if( fp_equal( (*TauLine2[i].Hi()).g(), realnum(10.) ) )
00391 branch_ratio = 1./6.;
00392 else
00393 TotalInsanity();
00394 pair<TransitionList::iterator,double> pp2( TauLine2.begin()+i, branch_ratio );
00395 C2Pump.push_back( pp2 );
00396 }
00397 }
00398 }
00399
00400
00401 double pump_rate = 0.;
00402 vector< pair<TransitionList::iterator,double> >::const_iterator c2p;
00403 for( c2p=C2Pump.begin(); c2p != C2Pump.end(); ++c2p )
00404 {
00405 const TransitionProxy::iterator t = c2p->first;
00406 double branch_ratio = c2p->second;
00407 pump_rate += (*t).Emis().pump()*branch_ratio;
00408 # if 0
00409 dprintf( ioQQQ, "C II %.3e %.3e\n",
00410 (*t).WLAng , (*t).Emis().pump()*branch_ratio );
00411 # endif
00412 }
00413
00414
00415
00416
00417
00418
00419 AtomSeqBoron(TauLines[ipT157],
00420 TauLines[ipC2_2325],
00421 TauLines[ipC2_2324],
00422 TauLines[ipC2_2329],
00423 TauLines[ipC2_2328],
00424 TauLines[ipC2_2327],
00425 0.2349 , 0.8237 , 0.8533 , 1.9818 , pump_rate , "C 2");
00426
00427
00428 enum {DEBUG_LOC=false};
00429 if( DEBUG_LOC && nzone > 80 )
00430 {
00431 fprintf(ioQQQ,"DEBUG\t%.2f\t%.3e\t%.3e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00432 fnzone ,
00433 phycon.te,
00434 TauLines[ipT157].Coll().cool() ,
00435 cse12,
00436 csh12,
00437 dense.eden,
00438 (dense.xIonDense[ipHYDROGEN][0] + findspecieslocal("H2")->den)/dense.cdsqte);
00439 }
00440
00441 double sum = 0.;
00442
00443 for( i=0; i < 5; ++i )
00444 {
00445 colden.C2Pops[i] = (realnum)atoms.PopLevels[i];
00446 sum += atoms.PopLevels[i];
00447 }
00448 if( dense.xIonDense[ipCARBON][1] > SMALLFLOAT )
00449 ASSERT( fabs(sum-dense.xIonDense[ipCARBON][1])/dense.xIonDense[ipCARBON][1] < 1e-4 );
00450
00451
00452 PutCS(.1,TauLines[ipT386]);
00453 atom_level2(TauLines[ipT386]);
00454
00455 PutCS(.1,TauLines[ipT310]);
00456 atom_level2(TauLines[ipT310]);
00457
00458 PutCS(.1,TauLines[ipT291]);
00459 atom_level2(TauLines[ipT291]);
00460
00461 PutCS(.1,TauLines[ipT280]);
00462 atom_level2(TauLines[ipT280]);
00463
00464 PutCS(.1,TauLines[ipT274]);
00465 atom_level2(TauLines[ipT274]);
00466
00467 PutCS(.1,TauLines[ipT270]);
00468 atom_level2(TauLines[ipT270]);
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478 cs = MIN2(1.1,2.67/phycon.te10);
00479 a21 = 5.149e-3;
00480 PutCS(cs,TauLines[ipT1909]);
00481
00482
00483
00484 AtomSeqBeryllium(.96,.73,2.8,TauLines[ipT1909],a21);
00485 embesq.em1908 = (realnum)(atoms.PopLevels[3]*a21*1.04e-11);
00486
00487
00488
00489
00490
00491 a21 = 6.87e-4;
00492
00493
00494 cs = 2.8*dense.cdsqte/5.*1.667;
00495 popratio = cs/(cs + a21);
00496 embesq.em13C1910 = (realnum)(a21 * atoms.PopLevels[1]*popratio* 1.04e-11 / co.C12_C13_isotope_ratio);
00497
00498
00499
00500
00501 popup = 0.;
00502 colden.C3Pops[0] = (realnum)atoms.PopLevels[0];
00503 for( i=1; i<4; ++i)
00504 {
00505 popup += atoms.PopLevels[i];
00506 colden.C3Pops[i] = (realnum)atoms.PopLevels[i];
00507 }
00508
00509 SaveAbun = dense.xIonDense[ipCARBON][2];
00510 dense.xIonDense[ipCARBON][2] = (realnum)popup;
00511
00512
00513
00514 cs = MIN2(30.,4.806*phycon.te10*phycon.te05/phycon.te01/phycon.te003);
00515 PutCS(18.45,TauLines[ipc31175]);
00516 atom_level2(TauLines[ipc31175]);
00517 dense.xIonDense[ipCARBON][2] = (realnum)SaveAbun;
00518
00519
00520
00521 cs = MIN2(7.0,1.556*phycon.te10);
00522 PutCS(cs,TauLines[ipT977]);
00523 atom_level2(TauLines[ipT977]);
00524
00525
00526
00527 ligbar(
00528 6,
00529 TauLines[ipT1548],
00530 TauLines[ipT312],
00531 &cs2s2p,&cs2s3p);
00532 PutCS(cs2s2p,TauLines[ipT1548]);
00533 PutCS(cs2s2p*0.5,TauLines[ipT1550]);
00534 PutCS(1.0,*TauDummy);
00535 atom_level3(
00536 TauLines[ipT1550],
00537 *TauDummy,
00538 TauLines[ipT1548]);
00539
00540 PutCS(cs2s3p,TauLines[ipT312]);
00541 atom_level2(TauLines[ipT312]);
00542 return;
00543 }