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