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