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