00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "hydrogenic.h"
00007 #include "taulines.h"
00008 #include "wind.h"
00009 #include "coolheavy.h"
00010 #include "radius.h"
00011 #include "conv.h"
00012 #include "h2.h"
00013 #include "rt.h"
00014 #include "doppvel.h"
00015 #include "opacity.h"
00016 #include "ionbal.h"
00017 #include "dense.h"
00018 #include "trace.h"
00019 #include "dynamics.h"
00020 #include "rfield.h"
00021 #include "grainvar.h"
00022 #include "atmdat.h"
00023 #include "atoms.h"
00024 #include "called.h"
00025 #include "mole.h"
00026 #include "hmi.h"
00027 #include "numderiv.h"
00028 #include "magnetic.h"
00029 #include "phycon.h"
00030 #include "lines_service.h"
00031 #include "hyperfine.h"
00032 #include "iso.h"
00033 #include "thermal.h"
00034 #include "cooling.h"
00035 #include "pressure.h"
00036
00037 STATIC void fndneg(void);
00038
00039 STATIC void fndstr(double tot,
00040 double dc);
00041
00042
00043 static const bool PRT_DERIV = false;
00044
00045 void CoolEvaluate(double *tot)
00046 {
00047 static long int nhit = 0,
00048 nzSave=0;
00049
00050 static double TeEvalCS = 0., TeEvalCS_21cm=0.;
00051 static double TeUsedBrems=-1.f;
00052 static int nzoneUsedBrems=-1;
00053
00054 static double electron_rate_21cm,
00055 atomic_rate_21cm,
00056 proton_rate_21cm;
00057
00058 double
00059 cs ,
00060 deriv,
00061 factor,
00062 qn,
00063 rothi=-SMALLFLOAT,
00064 rotlow=-SMALLFLOAT,
00065 x;
00066
00067 static double oltcool=0.,
00068 oldtemp=0.;
00069
00070 long int coolnum, coolcal;
00071
00072 DEBUG_ENTRY( "CoolEvaluate()" );
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086 if( trace.lgTrace )
00087 fprintf( ioQQQ, " COOLR TE:%.4e zone %li %li Cool:%.4e Heat:%.4e eden:%.4e edenTrue:%.4e\n",
00088 phycon.te,
00089 nzone, conv.nPres2Ioniz ,
00090 thermal.ctot , thermal.htot,dense.eden,dense.EdenTrue );
00091
00092
00093
00094 TempChange(phycon.te , false);
00095
00096
00097 CoolZero();
00098 if( PRT_DERIV )
00099 fprintf(ioQQQ,"DEBUG dCdT 0 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00100 if( gv.lgGrainPhysicsOn )
00101 {
00102
00103
00104
00105 CoolAdd("dust",0,MAX2(0.,gv.GasCoolColl));
00106
00107
00108 thermal.dCooldT += MAX2(0.,gv.GasCoolColl)*3./(2.*phycon.te);
00109
00110
00111
00112 if( gv.lgDustOn() && gv.lgDHetOn )
00113 {
00114
00115 thermal.heating[0][13] = gv.GasHeatPhotoEl;
00116
00117
00118
00119
00120 thermal.heating[0][14] = MAX2(0.,-gv.GasCoolColl);
00121
00122
00123 thermal.heating[0][25] = gv.GasHeatTherm;
00124 }
00125 else
00126 {
00127 thermal.heating[0][13] = 0.;
00128 thermal.heating[0][14] = 0.;
00129 thermal.heating[0][25] = 0.;
00130 }
00131 }
00132 else if( gv.lgBakesPAH_heat )
00133 {
00134
00135
00136 thermal.heating[0][13] = gv.GasHeatPhotoEl;
00137 }
00138
00139 if( PRT_DERIV )
00140 fprintf(ioQQQ,"DEBUG dCdT 1 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00141
00142
00143 if( mole_global.lgNoMole )
00144 {
00145
00146 hmi.hmicol = 0.;
00147 CoolHeavy.brems_cool_hminus = 0.;
00148
00149 CoolHeavy.h2line = 0.;
00150
00151 CoolHeavy.H2PlsCool = 0.;
00152 CoolHeavy.HD = 0.;
00153
00154
00155 thermal.heating[0][8] = 0.;
00156
00157 thermal.heating[0][15] = 0.;
00158
00159 thermal.heating[0][16] = 0.;
00160 hmi.HeatH2Dish_used = 0.;
00161 hmi.HeatH2Dexc_used = 0.;
00162 hmi.deriv_HeatH2Dexc_used = 0.;
00163 }
00164
00165 else
00166 {
00167
00168 thermal.heating[0][15] = hmi.hmihet;
00169 thermal.heating[0][16] = hmi.h2plus_heat;
00170
00171 if( h2.lgEnabled && hmi.lgH2_Thermal_BigH2 )
00172 {
00173 if( h2.lgEvaluated )
00174 {
00175
00176
00177
00178
00179 hmi.HeatH2Dish_used = h2.HeatDiss;
00180 hmi.HeatH2Dexc_used = h2.HeatDexc;
00181
00182
00183
00184
00185
00186 hmi.deriv_HeatH2Dexc_used = -h2.HeatDexc_deriv;
00187 }
00188 else
00189 {
00190 hmi.HeatH2Dish_used = 0;
00191 hmi.HeatH2Dexc_used = 0;
00192 hmi.deriv_HeatH2Dexc_used = 0;
00193 }
00194 }
00195
00196 else if( hmi.chH2_small_model_type == 'T' )
00197 {
00198
00199
00200 hmi.HeatH2Dish_used = hmi.HeatH2Dish_TH85;
00201 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_TH85;
00202 hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_TH85;
00203 }
00204 else if( hmi.chH2_small_model_type == 'H' )
00205 {
00206
00207 hmi.HeatH2Dish_used = hmi.HeatH2Dish_BHT90;
00208 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_BHT90;
00209 hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_BHT90;
00210 }
00211 else if( hmi.chH2_small_model_type == 'B')
00212 {
00213
00214 hmi.HeatH2Dish_used = hmi.HeatH2Dish_BD96;
00215 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_BD96;
00216 hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_BD96;
00217 }
00218 else if(hmi.chH2_small_model_type == 'E')
00219 {
00220
00221 hmi.HeatH2Dish_used = hmi.HeatH2Dish_ELWERT;
00222 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_ELWERT;
00223 hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_ELWERT;
00224 }
00225 else
00226 TotalInsanity();
00227
00228
00229 thermal.heating[0][17] = hmi.HeatH2Dish_used;
00230
00231
00232 thermal.heating[0][28] = 0.;
00233 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00234 {
00235 if( (*diatom)->lgEnabled && mole_global.lgStancil )
00236 thermal.heating[0][28] += (*diatom)->Cont_Diss_Heat_Rate();
00237 }
00238
00239
00240
00241 thermal.heating[0][8] = MAX2(0.,hmi.HeatH2Dexc_used);
00242
00243
00244 CoolAdd("H2cX",0,MAX2(0.,-hmi.HeatH2Dexc_used));
00245
00246
00247
00248
00249
00250
00251
00252 if( hmi.HeatH2Dexc_used < 0. )
00253 thermal.dCooldT -= hmi.deriv_HeatH2Dexc_used;
00254
00255
00256 CoolHeavy.H2PlsCool = (realnum)(MAX2((2.325*phycon.te-1875.)*1e-20,0.)*
00257 dense.xIonDense[ipHYDROGEN][0]*dense.xIonDense[ipHYDROGEN][1]*1.66e-11);
00258
00259 if( h2.lgEnabled )
00260 {
00261
00262
00263 CoolHeavy.h2line = 0.;
00264 }
00265 else
00266 {
00267
00268
00269 x = phycon.alogte - 4.;
00270 if( phycon.te > 1087. )
00271 {
00272 rothi = 3.90e-19*sexp(6118./phycon.te);
00273 }
00274 else
00275 {
00276 rothi = pow(10.,-19.24 + 0.474*x - 1.247*x*x);
00277 }
00278
00279
00280
00281 qn = pow(MAX2(hmi.H2_total,1e-37),0.77) + 1.2*pow(MAX2(dense.xIonDense[ipHYDROGEN][0],1e-37),0.77);
00282
00283 if( phycon.te > 4031. )
00284 {
00285 rotlow = 1.38e-22*sexp(9243./phycon.te)*qn;
00286 }
00287 else
00288 {
00289 rotlow = pow(10.,-22.90 - 0.553*x - 1.148*x*x)*qn;
00290 }
00291
00292 CoolHeavy.h2line = 0.;
00293 if( rotlow > 0. )
00294 CoolHeavy.h2line += hmi.H2_total*rothi/(1. + rothi/rotlow);
00295
00296
00297
00298 }
00299
00300 {
00301 enum {DEBUG_LOC=false};
00302 if( DEBUG_LOC && nzone>187&& iteration > 1)
00303 {
00304 fprintf(ioQQQ,"h2coolbug\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00305 phycon.te,
00306 CoolHeavy.h2line,
00307 hmi.H2_total,
00308 findspecieslocal("H-")->den,
00309 hmi.HMinus_photo_rate,
00310 rothi,
00311 rotlow );
00312 }
00313 }
00314
00315 if( hd.lgEnabled )
00316 {
00317 CoolHeavy.HD = 0.;
00318 }
00319 else
00320 {
00321
00322
00323 factor = sexp(128.6/phycon.te);
00324 CoolHeavy.HD = 2.66e-21 * hydro.D2H_ratio * POW2((double)hmi.H2_total) * phycon.sqrte *
00325 factor/(1416.+phycon.sqrte*hmi.H2_total * (1. + 3.*factor));
00326 }
00327 }
00328
00329 fixit();
00330 #if 0
00331 double chemical_heating = mole.chem_heat();
00332 thermal.heating[0][29] = MAX2(0.,chemical_heating);
00333
00334 CoolAdd("Chem",0,MAX2(0.,-chemical_heating));
00335 #endif
00336
00337
00338 CoolAdd("CT C" , 0. , thermal.char_tran_cool );
00339
00340
00341
00342 CoolAdd("H-fb",0,hmi.hmicol);
00343
00344
00345
00346
00347 thermal.dCooldT += 2.*hmi.hmicol*phycon.teinv;
00348 if( PRT_DERIV )
00349 fprintf(ioQQQ,"DEBUG dCdT 2 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00350
00351 CoolAdd("H2ln",0,CoolHeavy.h2line);
00352
00353
00354
00355
00356
00357
00358
00359
00360 thermal.dCooldT += 3.0*CoolHeavy.h2line*phycon.teinv;
00361
00362 {
00363
00364 enum {DEBUG_LOC=false};
00365 if( DEBUG_LOC )
00366 {
00367 fprintf(ioQQQ,"CoolEvaluate debuggg\t%.2f\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n",
00368 fnzone,
00369 phycon.te,
00370 hmi.H2_total ,
00371 CoolHeavy.h2line,
00372 findspecieslocal("H-")->den ,
00373 dense.eden);
00374 }
00375 }
00376
00377 CoolAdd("HDro",0,CoolHeavy.HD);
00378 thermal.dCooldT += CoolHeavy.HD*phycon.teinv;
00379
00380 CoolAdd("H2+ ",0,CoolHeavy.H2PlsCool);
00381 thermal.dCooldT += CoolHeavy.H2PlsCool*phycon.teinv;
00382
00383
00384 thermal.heating[0][3] = 0.;
00385
00386 thermal.heating[0][23] = 0.;
00387
00388 thermal.heating[0][1] = 0.;
00389
00390
00391 for( long int ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00392 {
00393 for( long int nelem=ipISO; nelem < LIMELM; nelem++ )
00394 {
00395
00396
00397 iso_cool( ipISO , nelem );
00398 }
00399 }
00400
00401
00402
00403
00404
00405
00406
00407
00408 if( fabs(CoolHeavy.brems_cool_net)/SDIV(thermal.ctot) > conv.HeatCoolRelErrorAllowed/5. ||
00409 conv.lgSearch ||
00410 !fp_equal(phycon.te,TeUsedBrems) ||
00411 nzone != nzoneUsedBrems )
00412 {
00413 double BremsThisEnergy;
00414
00415 long int limit;
00416
00417
00418 TeUsedBrems = phycon.te;
00419 nzoneUsedBrems = nzone;
00420
00421 limit = MIN2( rfield.ipMaxBolt , rfield.nflux );
00422
00423 CoolHeavy.brems_cool_hminus = 0.;
00424 CoolHeavy.brems_cool_h = 0.;
00425 CoolHeavy.brems_cool_metals = 0.;
00426 CoolHeavy.brems_cool_he = 0.;
00427 CoolHeavy.brems_heat_total = 0.;
00428
00429 {
00430 double bhfac, bhMinusfac;
00431 realnum sumion[LIMELM+1];
00432 long int ion_lo , ion_hi;
00433
00434 ASSERT(rfield.ipEnergyBremsThin < rfield.nupper);
00435 ASSERT(limit < rfield.nupper);
00436
00437
00438
00439
00440
00441
00442
00443 CoolHeavy.brems_cool_h = 0.;
00444 CoolHeavy.brems_cool_hminus = 0.;
00445
00446 for( long int i=rfield.ipEnergyBremsThin; i < limit; i++ )
00447 {
00448 long int ion = 1;
00449
00450
00451
00452 BremsThisEnergy = rfield.gff[ion][i]*rfield.widflx[i]*rfield.ContBoltz[i];
00453
00454 CoolHeavy.brems_cool_h += BremsThisEnergy;
00455
00456
00457
00458 CoolHeavy.brems_cool_hminus += BremsThisEnergy * opac.OpacStack[i-1+opac.iphmra];
00459 }
00460 bhfac = (dense.xIonDense[ipHYDROGEN][1]+findspecieslocal("H2+")->den+findspecieslocal("H3+")->den)*CoolHeavy.lgFreeOn* dense.eden*1.032e-11/phycon.sqrte*EN1RYD;
00461 bhMinusfac = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop()*CoolHeavy.lgFreeOn* dense.eden*1.032e-11/phycon.sqrte*EN1RYD;
00462 CoolHeavy.brems_cool_h *= bhfac;
00463 CoolHeavy.brems_cool_hminus *= bhMinusfac;
00464
00465
00466 CoolHeavy.brems_cool_he = 0.;
00467 for( long int i=rfield.ipEnergyBremsThin; i < limit; i++ )
00468 {
00469 long int nelem = ipHELIUM;
00470
00471 BremsThisEnergy =
00472 (dense.xIonDense[nelem][1]*rfield.gff[1][i] + 4.*dense.xIonDense[nelem][2]*rfield.gff[2][i])*
00473 rfield.widflx[i]*rfield.ContBoltz[i];
00474 CoolHeavy.brems_cool_he += BremsThisEnergy;
00475 }
00476 CoolHeavy.brems_cool_he *=
00477 CoolHeavy.lgFreeOn * dense.eden*1.032e-11/phycon.sqrte*EN1RYD;
00478
00479
00480
00481
00482 sumion[0] = 0.;
00483 for( long int ion=1; ion<=LIMELM; ++ion )
00484 {
00485 sumion[ion] = 0.;
00486 for( long int nelem=ipLITHIUM; nelem < LIMELM; ++nelem )
00487 {
00488 if( dense.lgElmtOn[nelem] && ion<=nelem+1 )
00489 {
00490 sumion[ion] += dense.xIonDense[nelem][ion];
00491 }
00492 }
00493
00494 sumion[ion] *= POW2((realnum)ion);
00495 }
00496
00497
00498 for( long ipMol = 0; ipMol<mole_global.num_calc; ipMol++ )
00499 {
00500 ASSERT( (mole_global.list[ipMol]->n_nuclei() != 1) ==
00501 (!mole_global.list[ipMol]->isMonatomic()));
00502
00503 if( !mole_global.list[ipMol]->isMonatomic() &&
00504 mole_global.list[ipMol]->parentLabel.empty() &&
00505 mole_global.list[ipMol]->charge > 0 &&
00506 mole_global.list[ipMol]->label != "H2+" &&
00507 mole_global.list[ipMol]->label != "H3+" )
00508 {
00509 ASSERT( mole_global.list[ipMol]->charge < LIMELM+1 );
00510 sumion[mole_global.list[ipMol]->charge] += (realnum)mole.species[ipMol].den * POW2((realnum)mole_global.list[ipMol]->charge);
00511 }
00512 }
00513
00514
00515
00516
00517
00518 ion_lo = 1;
00519 while( sumion[ion_lo]==0 && ion_lo<LIMELM-1 )
00520 ++ion_lo;
00521 ion_hi = LIMELM;
00522 while( sumion[ion_hi]==0 && ion_hi>0 )
00523 --ion_hi;
00524
00525
00526 CoolHeavy.brems_cool_metals = 0.;
00527 CoolHeavy.brems_heat_total = 0.;
00528 for( long int i=rfield.ipEnergyBremsThin; i < limit; i++ )
00529 {
00530 BremsThisEnergy = 0.;
00531 for(long int ion=ion_lo; ion<=ion_hi; ++ion )
00532 BremsThisEnergy += sumion[ion]*rfield.gff[ion][i];
00533
00534 CoolHeavy.brems_cool_metals += BremsThisEnergy*rfield.widflx[i]*rfield.ContBoltz[i];
00535
00536 CoolHeavy.brems_heat_total += opac.FreeFreeOpacity[i]*rfield.flux[0][i]*rfield.anu[i]*EN1RYD*CoolHeavy.lgFreeOn;
00537 }
00538 CoolHeavy.brems_cool_metals *=
00539 CoolHeavy.lgFreeOn * dense.eden*1.032e-11/phycon.sqrte*EN1RYD;
00540
00541 {
00542 enum {DEBUG_LOC=false};
00543 if( DEBUG_LOC && nzone>60 )
00544 {
00545 double sumfield = 0., sumtot=0., sum1=0., sum2=0.;
00546 for( long int i=rfield.ipEnergyBremsThin; i<limit; i++ )
00547 {
00548 sumtot += opac.FreeFreeOpacity[i]*rfield.flux[0][i]*rfield.anu[i];
00549 sumfield += rfield.flux[0][i]*rfield.anu[i];
00550 sum1 += opac.FreeFreeOpacity[i]*rfield.flux[0][i]*rfield.anu[i];
00551 sum2 += opac.FreeFreeOpacity[i]*rfield.flux[0][i];
00552 }
00553 fprintf(ioQQQ,"DEBUG brems heat\t%.2f\t%.3e\t%.3e\t%.3e\t%e\t%.3e\t%.3e\n",
00554 fnzone,
00555 CoolHeavy.brems_heat_total,
00556 sumtot/SDIV(sumfield) ,
00557 sum1/SDIV(sum2),
00558 phycon.te ,
00559 rfield.gff[1][1218],
00560 opac.FreeFreeOpacity[1218]);
00561 }
00562 }
00563 }
00564 }
00565
00566
00567 CoolHeavy.brems_cool_net =
00568 CoolHeavy.brems_cool_h +
00569 CoolHeavy.brems_cool_he +
00570 CoolHeavy.brems_cool_hminus +
00571 CoolHeavy.brems_cool_metals -
00572 CoolHeavy.brems_heat_total;
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584 CoolAdd( "FF c" , 0, MAX2(0.,CoolHeavy.brems_cool_net) );
00585
00586
00587 thermal.heating[0][11] = MAX2(0.,-CoolHeavy.brems_cool_net);
00588
00589
00590
00591 thermal.dCooldT += CoolHeavy.brems_cool_h*thermal.halfte;
00592 if( PRT_DERIV )
00593 fprintf(ioQQQ,"DEBUG dCdT 3 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00594
00595
00596
00597
00598
00599
00600
00601 CoolHeavy.heavfb = 0.;
00602 for( long int nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
00603 {
00604 if( dense.lgElmtOn[nelem] )
00605 {
00606
00607 long limit_lo = MAX2( 1 , dense.IonLow[nelem] );
00608 long limit_hi = MIN2( nelem-NISO+1, dense.IonHigh[nelem] );
00609 for( long int ion=limit_lo; ion<=limit_hi; ++ion )
00610 {
00611
00612
00613
00614
00615
00616
00617 double one = dense.xIonDense[nelem][ion] * ionbal.RR_rate_coef_used[nelem][ion-1]*
00618 dense.eden * phycon.te * BOLTZMANN;
00619
00620
00621 CoolHeavy.heavfb += one;
00622 }
00623 }
00624 }
00625
00626
00627
00628 CoolAdd("hvFB",0,CoolHeavy.heavfb);
00629 thermal.dCooldT += CoolHeavy.heavfb*.113*phycon.teinv;
00630 if( PRT_DERIV )
00631 fprintf(ioQQQ,"DEBUG dCdT 4 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00632
00633
00634
00635
00636 CoolHeavy.eebrm = POW2(dense.eden*phycon.te*1.84e-21);
00637
00638
00639 thermal.dCooldT += CoolHeavy.eebrm*thermal.halfte;
00640 CoolAdd("eeff",0,CoolHeavy.eebrm);
00641
00642
00643
00644 CoolAdd("adve",0,dynamics.Cool() );
00645
00646 thermal.dCooldT += dynamics.dCooldT();
00647
00648 thermal.heating[1][5] = dynamics.Heat();
00649 thermal.dHeatdT += dynamics.dHeatdT;
00650
00651
00652 CoolHeavy.tccool = rfield.cmcool*phycon.te;
00653 CoolAdd("Comp",0,CoolHeavy.tccool);
00654 thermal.dCooldT += rfield.cmcool;
00655
00656
00657
00658 if( thermal.lgCExtraOn )
00659 {
00660 CoolHeavy.cextxx =
00661 (realnum)(thermal.CoolExtra*pow(phycon.te/1e4,(double)thermal.cextpw));
00662 }
00663 else
00664 {
00665 CoolHeavy.cextxx = 0.;
00666 }
00667 CoolAdd("Extr",0,CoolHeavy.cextxx);
00668
00669 realnum dDensityDT;
00670
00671
00672 if( wind.lgBallistic() )
00673 {
00674 dDensityDT = -(realnum)(wind.AccelTotalOutward/wind.windv + 2.*wind.windv/
00675 radius.Radius);
00676 CoolHeavy.expans = -2.5*pressure.PresGasCurr*dDensityDT;
00677 }
00678 else if( dynamics.lgTimeDependentStatic &&
00679 iteration > dynamics.n_initial_relax)
00680 {
00681 realnum dens = scalingDensity();
00682 dDensityDT =
00683 (realnum)((dens-dynamics.Upstream_density)/
00684 (dynamics.timestep*0.5*(dens+dynamics.Upstream_density)));
00685
00686 CoolHeavy.expans = -pressure.PresGasCurr*dDensityDT;
00687 }
00688 else
00689 {
00690 dDensityDT = 0.;
00691 CoolHeavy.expans = 0.;
00692 }
00693 CoolAdd("Expn",0,CoolHeavy.expans);
00694 thermal.dCooldT += CoolHeavy.expans/phycon.te;
00695
00696
00697
00698 CoolHeavy.cyntrn = 4.5433e-25f*magnetic.pressure*PI8*dense.eden*phycon.te;
00699 CoolAdd("Cycl",0,CoolHeavy.cyntrn);
00700 thermal.dCooldT += CoolHeavy.cyntrn/phycon.te;
00701
00702
00703
00704
00705 CoolAdd("Hvin",0,CoolHeavy.colmet);
00706
00707
00708 coolnum = thermal.ncltot;
00709 if( !fp_equal(phycon.te,TeEvalCS_21cm) )
00710 {
00711 {
00712
00713 enum {DEBUG_LOC=false};
00714 if( DEBUG_LOC )
00715 {
00716 # define N21CM_TE 16
00717 int n;
00718 double teval[N21CM_TE]={2.,5.,10.,20.,50.,100.,200.,500.,1000.,
00719 2000.,3000.,5000.,7000.,10000.,15000.,20000.};
00720 for( n = 0; n<N21CM_TE; ++n )
00721 {
00722 fprintf(
00723 ioQQQ,"DEBUG 21 cm deex Te=\t%.2e\tH0=\t%.2e\tp=\t%.2e\te=\t%.2e\n",
00724 teval[n],
00725 H21cm_H_atom( teval[n] ),
00726 H21cm_proton( teval[n] ),
00727 H21cm_electron( teval[n] ) );
00728 }
00729 cdEXIT(EXIT_FAILURE);
00730 # undef N21CM_TE
00731 }
00732 }
00733
00734
00735 atomic_rate_21cm = H21cm_H_atom( phycon.te );
00736 proton_rate_21cm = H21cm_proton( phycon.te );
00737 electron_rate_21cm = H21cm_electron( phycon.te );
00738 TeEvalCS_21cm = phycon.te;
00739 }
00740
00741
00742 cs = (electron_rate_21cm * dense.eden +
00743 atomic_rate_21cm * dense.xIonDense[ipHYDROGEN][0] +
00744 proton_rate_21cm * dense.xIonDense[ipHYDROGEN][1] ) *
00745 3./ dense.cdsqte;
00746 PutCS( cs , HFLines[0] );
00747
00748
00749 if( !fp_equal(phycon.te,TeEvalCS) )
00750 {
00751
00752 for( long int i=1; i < nHFLines; i++ )
00753 {
00754 cs = HyperfineCS( i );
00755
00756 PutCS( cs , HFLines[i] );
00757 }
00758 TeEvalCS = phycon.te;
00759 }
00760
00761
00762 RT_line_one( HFLines[0], true,0.f, GetDopplerWidth(dense.AtomicWeight[(*HFLines[0].Hi()).nelem()-1]) );
00763 H21_cm_pops();
00764 if( PRT_DERIV )
00765 fprintf(ioQQQ,"DEBUG dCdT 5 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00766
00767
00768 hyperfine.cooling_total = HFLines[0].Coll().cool();
00769
00770
00771 for( long int i=1; i < nHFLines; i++ )
00772 {
00773
00774 realnum save = dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1];
00775
00776
00777 if( save<=0. )
00778 continue;
00779
00780 RT_line_one( HFLines[i], true,0.f, GetDopplerWidth(dense.AtomicWeight[(*HFLines[i].Hi()).nelem()-1]) );
00781
00782
00783 dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1] *=
00784 hyperfine.HFLabundance[i];
00785
00786
00787 atom_level2( HFLines[i] );
00788
00789
00790 dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1] = save;
00791
00792
00793 hyperfine.cooling_total += HFLines[i].Coll().cool();
00794 }
00795 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00796 thermal.elementcool[ipHYDROGEN] += thermal.cooling[coolcal];
00797
00798 if( PRT_DERIV )
00799 fprintf(ioQQQ,"DEBUG dCdT 6 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00800
00801 double xIonDenseSave[LIMELM][LIMELM+1];
00802 if( atmdat.lgChiantiOn ||atmdat.lgStoutOn)
00803 {
00804 for( int nelem=0; nelem < LIMELM; nelem++ )
00805 {
00806 for( int ion=0; ion<=nelem+1; ++ion )
00807 {
00808 xIonDenseSave[nelem][ion] = dense.xIonDense[nelem][ion];
00809
00810 if( dense.lgIonChiantiOn[nelem][ion] || dense.lgIonStoutOn[nelem][ion] )
00811 dense.xIonDense[nelem][ion] = 0.;
00812 }
00813 }
00814 }
00815
00816
00817 coolnum = thermal.ncltot;
00818 CoolCarb();
00819 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00820 thermal.elementcool[ipCARBON] += thermal.cooling[coolcal];
00821 if( PRT_DERIV )
00822 fprintf(ioQQQ,"DEBUG dCdT C %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00823
00824
00825 coolnum = thermal.ncltot;
00826 CoolNitr();
00827 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00828 thermal.elementcool[ipNITROGEN] += thermal.cooling[coolcal];
00829 if( PRT_DERIV )
00830 fprintf(ioQQQ,"DEBUG dCdT N %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00831
00832
00833 coolnum = thermal.ncltot;
00834 CoolOxyg();
00835 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00836 thermal.elementcool[ipOXYGEN] += thermal.cooling[coolcal];
00837 if( PRT_DERIV )
00838 fprintf(ioQQQ,"DEBUG dCdT 7 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00839
00840
00841 coolnum = thermal.ncltot;
00842 CoolNeon();
00843 if( PRT_DERIV )
00844 fprintf(ioQQQ,"DEBUG dCdT Ne %.3e dHdT %.3e\n",thermal.dCooldT
00845 , thermal.dHeatdT);
00846 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00847 thermal.elementcool[ipNEON] += thermal.cooling[coolcal];
00848
00849
00850 coolnum = thermal.ncltot;
00851 CoolMagn();
00852 if( PRT_DERIV )
00853 fprintf(ioQQQ,"DEBUG dCdT 8 %.3e dHdT %.3e\n",thermal.dCooldT
00854 , thermal.dHeatdT);
00855 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00856 thermal.elementcool[ipMAGNESIUM] += thermal.cooling[coolcal];
00857
00858
00859 coolnum = thermal.ncltot;
00860 CoolSodi();
00861 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00862 thermal.elementcool[ipSODIUM] += thermal.cooling[coolcal];
00863 if( PRT_DERIV )
00864 fprintf(ioQQQ,"DEBUG dCdT Na %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00865
00866
00867 coolnum = thermal.ncltot;
00868 CoolAlum();
00869 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00870 thermal.elementcool[ipALUMINIUM] += thermal.cooling[coolcal];
00871 if( PRT_DERIV )
00872 fprintf(ioQQQ,"DEBUG dCdT Al %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00873
00874
00875 coolnum = thermal.ncltot;
00876 CoolSili();
00877 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00878 thermal.elementcool[ipSILICON] += thermal.cooling[coolcal];
00879 if( PRT_DERIV )
00880 fprintf(ioQQQ,"DEBUG dCdT 9 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00881
00882
00883 coolnum = thermal.ncltot;
00884 CoolPhos();
00885 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00886 thermal.elementcool[ipPHOSPHORUS] += thermal.cooling[coolcal];
00887
00888
00889 coolnum = thermal.ncltot;
00890 CoolSulf();
00891 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00892 thermal.elementcool[ipSULPHUR] += thermal.cooling[coolcal];
00893
00894
00895 coolnum = thermal.ncltot;
00896 CoolChlo();
00897 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00898 thermal.elementcool[ipCHLORINE] += thermal.cooling[coolcal];
00899
00900
00901 coolnum = thermal.ncltot;
00902 CoolArgo();
00903 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00904 thermal.elementcool[ipARGON] += thermal.cooling[coolcal];
00905 if( PRT_DERIV )
00906 fprintf(ioQQQ,"DEBUG dCdT 10 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00907
00908
00909 coolnum = thermal.ncltot;
00910 CoolPota();
00911 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00912 thermal.elementcool[ipPOTASSIUM] += thermal.cooling[coolcal];
00913
00914
00915 coolnum = thermal.ncltot;
00916 CoolCalc();
00917 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00918 thermal.elementcool[ipCALCIUM] += thermal.cooling[coolcal];
00919
00920
00921 coolnum = thermal.ncltot;
00922 CoolScan();
00923 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00924 thermal.elementcool[ipSCANDIUM] += thermal.cooling[coolcal];
00925
00926
00927 coolnum = thermal.ncltot;
00928 CoolChro();
00929 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00930 thermal.elementcool[ipCHROMIUM] += thermal.cooling[coolcal];
00931
00932
00933
00934 coolnum = thermal.ncltot;
00935 CoolIron();
00936 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00937 thermal.elementcool[ipIRON] += thermal.cooling[coolcal];
00938 if( PRT_DERIV )
00939 fprintf(ioQQQ,"DEBUG dCdT 12 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
00940
00941
00942 coolnum = thermal.ncltot;
00943 CoolCoba();
00944 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00945 thermal.elementcool[ipCOBALT] += thermal.cooling[coolcal];
00946
00947
00948 coolnum = thermal.ncltot;
00949 CoolNick();
00950 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00951 thermal.elementcool[ipNICKEL] += thermal.cooling[coolcal];
00952
00953 coolnum = thermal.ncltot;
00954
00955
00956 if( atmdat.lgChiantiOn || atmdat.lgStoutOn)
00957 {
00958
00959 for( int nelem=0; nelem < LIMELM; nelem++ )
00960 {
00961 for( int ion=0; ion<=nelem+1; ++ion )
00962 {
00963 dense.xIonDense[nelem][ion] = xIonDenseSave[nelem][ion];
00964 }
00965 }
00966 }
00967
00968
00969 CoolDima();
00970
00971 for( int coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00972 thermal.dima += thermal.cooling[coolcal];
00973
00974
00975 dBase_solve();
00976
00977
00978 {
00979 enum {DEBUG_LOC=false};
00980 if( DEBUG_LOC )
00981 {
00982 static bool lgMustPrintHeader=true;
00983 if( lgMustPrintHeader )
00984 {
00985 lgMustPrintHeader = false;
00986 printf("DEBUG Levels\t%s",dBaseSpecies[0].chLabel );
00987 for( long ipSpecies=1; ipSpecies<nSpecies; ipSpecies++ )
00988 {
00989 printf("\t%s",dBaseSpecies[ipSpecies].chLabel );
00990 }
00991 printf("\n" );
00992 printf("DEBUG Max\t%li" ,dBaseSpecies[0].numLevels_max );
00993 for( long ipSpecies=1; ipSpecies<nSpecies; ipSpecies++ )
00994 {
00995 printf( "\t%li" ,dBaseSpecies[ipSpecies].numLevels_max );
00996 }
00997 printf("\n");
00998 }
00999 printf("DEBUG Local\t%li" ,dBaseSpecies[0].numLevels_local );
01000 for( long ipSpecies=1; ipSpecies<nSpecies; ipSpecies++ )
01001 {
01002 printf("\t%li" ,dBaseSpecies[ipSpecies].numLevels_local );
01003 }
01004 printf("\n");
01005 }
01006 }
01007
01008
01009 CoolSum(tot);
01010 if( PRT_DERIV )
01011 fprintf(ioQQQ,"DEBUG dCdT 14 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
01012
01013
01014 if( *tot <= 0. )
01015 {
01016 fprintf( ioQQQ, " COOLR; cooling is <=0, this is impossible.\n" );
01017 ShowMe();
01018 cdEXIT(EXIT_FAILURE);
01019 }
01020
01021
01022 if( thermal.dCooldT == 0. )
01023 {
01024 fprintf( ioQQQ, " COOLR; cooling slope <=0, this is impossible.\n" );
01025 if( *tot > 0. && dense.gas_phase[ipHYDROGEN] < 1e-4 )
01026 {
01027 fprintf( ioQQQ, " Probably due to very low density.\n" );
01028 }
01029 ShowMe();
01030 cdEXIT(EXIT_FAILURE);
01031 }
01032
01033 if( trace.lgTrace )
01034 {
01035 fndstr(*tot,thermal.dCooldT);
01036 }
01037
01038
01039 if( (((((!thermal.lgTemperatureConstant) && *tot < 0.) && called.lgTalk) &&
01040 !conv.lgSearch) && thermal.lgCNegChk) && nzone > 0 )
01041 {
01042 fprintf( ioQQQ,
01043 " NOTE Negative cooling, zone %4ld, =%10.2e coola=%10.2e CHION=%10.2e Te=%10.2e\n",
01044 nzone,
01045 *tot,
01046 iso_sp[ipH_LIKE][ipHYDROGEN].cLya_cool,
01047 iso_sp[ipH_LIKE][ipHYDROGEN].coll_ion,
01048 phycon.te );
01049 fndneg();
01050 }
01051
01052
01053
01054 if( NumDeriv.lgNumDeriv )
01055 {
01056 if( ((nzone > 2 && nzone == nzSave) && ! fp_equal( oldtemp, phycon.te )) && nhit > 4 )
01057 {
01058
01059
01060 deriv = (oltcool - *tot)/(oldtemp - phycon.te);
01061 thermal.dCooldT = deriv;
01062 }
01063 else
01064 {
01065 deriv = thermal.dCooldT;
01066 }
01067 if( nzone != nzSave )
01068 nhit = 0;
01069
01070 nzSave = nzone;
01071 nhit += 1;
01072 oltcool = *tot;
01073 oldtemp = phycon.te;
01074 }
01075 return;
01076 }
01077
01078
01079 #ifdef EPS
01080 # undef EPS
01081 #endif
01082 #define EPS 0.01
01083
01084
01085 STATIC void fndneg(void)
01086 {
01087 long int i;
01088 double trig;
01089
01090 DEBUG_ENTRY( "fndneg()" );
01091
01092 trig = fabs(thermal.htot*EPS);
01093 for( i=0; i < thermal.ncltot; i++ )
01094 {
01095 if( thermal.cooling[i] < 0. && fabs(thermal.cooling[i]) > trig )
01096 {
01097 fprintf( ioQQQ, " negative line=%s %.2f fraction of heating=%.3f\n",
01098 thermal.chClntLab[i], thermal.collam[i], thermal.cooling[i]/
01099 thermal.htot );
01100 }
01101
01102 if( thermal.heatnt[i] > trig )
01103 {
01104 fprintf( ioQQQ, " heating line=%s %.2f fraction of heating=%.3f\n",
01105 thermal.chClntLab[i], thermal.collam[i], thermal.heatnt[i]/
01106 thermal.htot );
01107 }
01108 }
01109 return;
01110 }
01111
01112
01113 STATIC void fndstr(double tot,
01114 double dc)
01115 {
01116 char chStrngLab[NCOLNT_LAB_LEN+1];
01117 long int i;
01118 realnum wl;
01119 double str,
01120 strong;
01121
01122 DEBUG_ENTRY( "fndstr()" );
01123
01124 strong = 0.;
01125 wl = -FLT_MAX;
01126 for( i=0; i < thermal.ncltot; i++ )
01127 {
01128 if( fabs(thermal.cooling[i]) > strong )
01129 {
01130
01131 wl = thermal.collam[i];
01132
01133
01134 ASSERT( strlen( thermal.chClntLab[i] ) <= NCOLNT_LAB_LEN );
01135 strcpy( chStrngLab, thermal.chClntLab[i] );
01136 strong = fabs(thermal.cooling[i]);
01137 }
01138 }
01139
01140 str = strong;
01141
01142 fprintf( ioQQQ,
01143 " fndstr cool: TE=%10.4e Ne=%10.4e C=%10.3e dC/dT=%10.2e ABS(%s %.1f)=%.2e nz=%ld\n",
01144 phycon.te, dense.eden, tot, dc, chStrngLab
01145 , wl, str, nzone );
01146
01147
01148 if( trace.lgCoolTr )
01149 {
01150 realnum ratio;
01151
01152
01153 coolpr(ioQQQ,(char*)thermal.chClntLab[0],1,0.,"ZERO");
01154
01155
01156 for( i=0; i < thermal.ncltot; i++ )
01157 {
01158
01159
01160 ratio = (realnum)(thermal.cooling[i]/thermal.ctot);
01161 if( ratio >= EPS )
01162 {
01163
01164 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i], ratio,"DOIT");
01165 }
01166 }
01167
01168
01169 coolpr(ioQQQ,"DONE",1,0.,"DONE");
01170
01171
01172 if( thermal.heating[0][22]/thermal.ctot > 0.05 )
01173 {
01174 fprintf( ioQQQ,
01175 " All coolant heat greater than%6.2f%% of the total will be printed.\n",
01176 EPS*100. );
01177
01178 coolpr(ioQQQ,"ZERO",1,0.,"ZERO");
01179 for( i=0; i < thermal.ncltot; i++ )
01180 {
01181 ratio = (realnum)(thermal.heatnt[i]/thermal.ctot);
01182 if( fabs(ratio) >=EPS )
01183 {
01184 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i],
01185 ratio,"DOIT");
01186 }
01187 }
01188 coolpr(ioQQQ,"DONE",1,0.,"DONE");
01189 }
01190 }
01191 return;
01192 }