cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*CoolEvaluate main routine to call others, to evaluate total cooling */ 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 /*fndneg search cooling array to find negative values */ 00032 STATIC void fndneg(void); 00033 /*fndstr search cooling stack to find strongest values */ 00034 STATIC void fndstr(double tot, 00035 double dc); 00036 00037 /* set true to debug derivative of heating and cooling */ 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 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 DEBUG_ENTRY( "CoolEvaluate()" ); 00071 00072 /* returns tot, the total cooling, 00073 * and dc, the derivative of the cooling */ 00074 00075 /* routine atom_level2( t10 ) 00076 * routine atom_level3( abund , t10,t21,t20) 00077 * tsq1 = 1. / (te**2) 00078 * POPEXC( O12,g1,g2,A21,excit,abund); result already*a21 00079 * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) 00080 * AtomSeqBeryllium(cs23,cs24,cs34,tarray,a41) 00081 * FIVEL( G(1-5) , ex(wn,1-5), cs12,cs13,14,15,23,24,25,34,35,45, 00082 * A21,31,41,51,32,42,52,43,53,54, pop(1-5), abund) */ 00083 00084 if( trace.lgTrace ) 00085 fprintf( ioQQQ, " COOLR TE:%.4e zone %li %li Cool:%.4e Heat:%.4e eden:%.4e edenTrue:%.4e\n", 00086 phycon.te, 00087 nzone, conv.nPres2Ioniz , 00088 thermal.ctot , thermal.htot,dense.eden,dense.EdenTrue ); 00089 00090 /* must call TempChange since ionization has changed, there are some 00091 * terms that affect collision rates (H0 term in electron collision) */ 00092 TempChange(phycon.te , false); 00093 00094 /* now zero out the cooling stack */ 00095 CoolZero(); 00096 if( PRT_DERIV ) 00097 fprintf(ioQQQ,"DEBUG dcdt 0 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00098 if( gv.lgGrainPhysicsOn ) 00099 { 00100 /* grain heating and cooling */ 00101 /* grain recombination cooling, evaluated elsewhere 00102 * can either heat or cool the gas, do cooling here */ 00103 CoolAdd("dust",0,MAX2(0.,gv.GasCoolColl)); 00104 00105 /* grain cooling proportional to temperature ^3/2 */ 00106 thermal.dCooldT += MAX2(0.,gv.GasCoolColl)*3./(2.*phycon.te); 00107 00108 /* these are the various heat agents from grains */ 00109 /* options to force gas heating or cooling by grains to zero - for tests only ! */ 00110 if( gv.lgDustOn && gv.lgDHetOn ) 00111 { 00112 /* rate dust heats gas by photoelectric effect */ 00113 thermal.heating[0][13] = gv.GasHeatPhotoEl; 00114 00115 /* if grains hotter than gas then collisions with gas act 00116 * to heat the gas, add this in here 00117 * a symmetric statement appears in COOLR, where cooling is added on */ 00118 thermal.heating[0][14] = MAX2(0.,-gv.GasCoolColl); 00119 00120 /* this is gas heating due to thermionic emissions */ 00121 thermal.heating[0][25] = gv.GasHeatTherm; 00122 } 00123 else 00124 { 00125 thermal.heating[0][13] = 0.; 00126 thermal.heating[0][14] = 0.; 00127 thermal.heating[0][25] = 0.; 00128 } 00129 } 00130 else if( gv.lgBakesPAH_heat ) 00131 { 00132 /* >>chng 06 jul 21, option to include Bakes PAH hack with grain physics off, 00133 * needed to test dynamics models */ 00134 thermal.heating[0][13] = gv.GasHeatPhotoEl; 00135 } 00136 00137 /* find net heating - cooling due to large H2 molecule */ 00138 H2_Cooling("CoolEvaluate"); 00139 if( PRT_DERIV ) 00140 fprintf(ioQQQ,"DEBUG dcdt 1 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00141 00142 /* molecular molecules molecule cooling */ 00143 if( hmi.lgNoH2Mole ) 00144 { 00145 /* this branch - do not include molecules */ 00146 hmi.hmicol = 0.; 00147 CoolHeavy.brems_cool_hminus = 0.; 00148 /* line cooling within simple H2 molecule - zero when big used */ 00149 CoolHeavy.h2line = 0.; 00150 /* H + H+ => H2+ cooling */ 00151 CoolHeavy.H2PlsCool = 0.; 00152 CoolHeavy.HD = 0.; 00153 00154 /* thermal.heating[0][8] is heating due to collisions within X of H2 */ 00155 thermal.heating[0][8] = 0.; 00156 /* thermal.heating[0][15] is H minus heating*/ 00157 thermal.heating[0][15] = 0.; 00158 /* thermal.heating[0][16] is H2+ heating */ 00159 thermal.heating[0][16] = 0.; 00160 hmi.HeatH2Dish_used = 0.; 00161 hmi.HeatH2Dexc_used = 0.; 00162 } 00163 00164 else 00165 { 00166 /* save various molecular heating/cooling agent */ 00167 thermal.heating[0][15] = hmi.hmihet; 00168 thermal.heating[0][16] = hmi.h2plus_heat; 00169 /* now get heating from H2 molecule, either simple or from big one */ 00170 if( h2.lgH2ON && hmi.lgH2_Thermal_BigH2 ) 00171 { 00172 if( hmi.lgBigH2_evaluated ) 00173 { 00174 /* these are explicitly from big H2 molecule, 00175 * first is heating due to radiative pump of excited states, followed by 00176 * radiative decay into continuum of X, followed by dissociation of molecule 00177 * with kinetic energy, typically 0.25 - 0.5 eV per event */ 00178 hmi.HeatH2Dish_used = hmi.HeatH2Dish_BigH2; 00179 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_BigH2; 00180 /*fprintf(ioQQQ,"DEBUG big %.2f\t%.5e\t%.2e\t%.2e\t%.2e\n", 00181 fnzone , phycon.te, hmi.HeatH2Dexc_used, 00182 hmi.H2_total , hmi.H2_star_BigH2);*/ 00183 /* negative sign because right term is really deriv of heating, 00184 * but will be used below as deriv of cooling */ 00185 hmi.deriv_HeatH2Dexc_used = -hmi.deriv_HeatH2Dexc_BigH2; 00186 } 00187 else 00188 { 00189 hmi.HeatH2Dish_used = 0; 00190 hmi.HeatH2Dexc_used = 0; 00191 hmi.deriv_HeatH2Dexc_used = 0; 00192 } 00193 } 00194 00195 else if( hmi.chH2_small_model_type == 'T' ) 00196 { 00197 /* TH85 dissociation heating */ 00198 /* these come from approximations in TH85, see comments above */ 00199 hmi.HeatH2Dish_used = hmi.HeatH2Dish_TH85; 00200 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_TH85; 00201 hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_TH85; 00202 } 00203 else if( hmi.chH2_small_model_type == 'H' ) 00204 { 00205 /* Burton et al. 1990 */ 00206 hmi.HeatH2Dish_used = hmi.HeatH2Dish_BHT90; 00207 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_BHT90; 00208 hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_BHT90; 00209 } 00210 else if( hmi.chH2_small_model_type == 'B') 00211 { 00212 /* Bertoldi & Draine */ 00213 hmi.HeatH2Dish_used = hmi.HeatH2Dish_BD96; 00214 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_BD96; 00215 hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_BD96; 00216 } 00217 else if(hmi.chH2_small_model_type == 'E') 00218 { 00219 /* this is the default when small H2 used */ 00220 hmi.HeatH2Dish_used = hmi.HeatH2Dish_ELWERT; 00221 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_ELWERT; 00222 hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_ELWERT; 00223 } 00224 else 00225 TotalInsanity(); 00226 00227 /* heating due to photodissociation heating */ 00228 thermal.heating[0][17] = hmi.HeatH2Dish_used; 00229 00230 /* heating (usually cooling in big H2) due to collisions within X */ 00231 /* add to heating is net heating is positive */ 00232 thermal.heating[0][8] = MAX2(0.,hmi.HeatH2Dexc_used); 00233 00234 /* add to cooling if net heating is negative */ 00235 CoolAdd("H2cX",0,MAX2(0.,-hmi.HeatH2Dexc_used)); 00236 /*fprintf(ioQQQ,"DEBUG coolh2\t%.2f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n", 00237 fnzone, phycon.te, dense.eden, hmi.H2_total, thermal.ctot, -hmi.HeatH2Dexc_used );*/ 00238 /* add to net derivative */ 00239 /*thermal.dCooldT += MAX2(0.,-hmi.HeatH2Dexc_used)* ( 30172. * thermal.tsq1 - thermal.halfte );*/ 00240 /* >>chng 04 jan 25, check sign to prevent cooling from entering here, 00241 * also enter neg sign since going into cooling stack (bug), in heatsum 00242 * same term adds to deriv of heating */ 00243 if( hmi.HeatH2Dexc_used < 0. ) 00244 thermal.dCooldT -= hmi.deriv_HeatH2Dexc_used; 00245 00246 /* H + H+ => H2+ cooling */ 00247 CoolHeavy.H2PlsCool = (realnum)(MAX2((2.325*phycon.te-1875.)*1e-20,0.)* 00248 dense.xIonDense[ipHYDROGEN][0]*dense.xIonDense[ipHYDROGEN][1]*1.66e-11); 00249 00250 if( h2.lgH2ON ) 00251 { 00252 /* this is simplified approximation to H2 rotation cooling, 00253 * big molecule goes this far better */ 00254 CoolHeavy.h2line = 0.; 00255 } 00256 else 00257 { 00258 /* rate for rotation lines from 00259 * >>refer h2 cool Lepp, S., & Shull, J.M. 1983, ApJ, 270, 578 */ 00260 x = phycon.alogte - 4.; 00261 if( phycon.te > 1087. ) 00262 { 00263 rothi = 3.90e-19*sexp(6118./phycon.te); 00264 } 00265 else 00266 { 00267 rothi = pow(10.,-19.24 + 0.474*x - 1.247*x*x); 00268 } 00269 00270 /* low density rotation cooling */ 00271 /*&qn = pow(MAX2(hmi.Hmolec[ipMH2g],1e-37),0.77) + 1.2*pow(MAX2(dense.xIonDense[ipHYDROGEN][0],1e-37),0.77);*/ 00272 qn = pow(MAX2(hmi.H2_total,1e-37),0.77) + 1.2*pow(MAX2(dense.xIonDense[ipHYDROGEN][0],1e-37),0.77); 00273 /* these are equations 11 from LS83 */ 00274 if( phycon.te > 4031. ) 00275 { 00276 rotlow = 1.38e-22*sexp(9243./phycon.te)*qn; 00277 } 00278 else 00279 { 00280 rotlow = pow(10.,-22.90 - 0.553*x - 1.148*x*x)*qn; 00281 } 00282 00283 if( rotlow > 0. ) 00284 { 00285 /*CoolHeavy.h2line = hmi.Hmolec[ipMH2g]*rothi/(1. + rothi/rotlow);*/ 00286 CoolHeavy.h2line = hmi.H2_total*rothi/(1. + rothi/rotlow); 00287 } 00288 else 00289 { 00290 CoolHeavy.h2line = 0.; 00291 } 00292 } 00293 00294 { 00295 enum {DEBUG_LOC=false}; 00296 if( DEBUG_LOC && nzone>187&& iteration > 1) 00297 { 00298 fprintf(ioQQQ,"h2coolbug\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 00299 phycon.te, 00300 CoolHeavy.h2line, 00301 hmi.H2_total, 00302 hmi.Hmolec[ipMHm], 00303 hmi.HMinus_photo_rate, 00304 rothi, 00305 rotlow ); 00306 } 00307 } 00308 00309 /* >>chng 02 mar 07, add DH cooling using rates (eqn 6) from 00310 * >>refer HD cooling Puy, D., Grenacher, L, & Jetzer, P., 1999, A&A, 345, 723 */ 00311 factor = sexp(128.6/phycon.te); 00312 /*CoolHeavy.HD = 2.66e-21 * hydro.D2H_ratio * POW2(hmi.Hmolec[ipMH2g]) * phycon.sqrte * 00313 factor/(1416.+phycon.sqrte*hmi.Hmolec[ipMH2g] * (1. + 3.*factor));*/ 00314 CoolHeavy.HD = 2.66e-21 * hydro.D2H_ratio * POW2((double)hmi.H2_total) * phycon.sqrte * 00315 /*factor/(1416.+phycon.sqrte*hmi.Hmolec[ipMH2g] * (1. + 3.*factor));*/ 00316 factor/(1416.+phycon.sqrte*hmi.H2_total * (1. + 3.*factor)); 00317 } 00318 00319 /* cooling due to charge transfer ionization / recombination */ 00320 CoolAdd("CT C" , 0. , thermal.char_tran_cool ); 00321 00322 /* H- FB; H + e -> H- + hnu */ 00323 /* H- FF is in with H ff */ 00324 CoolAdd("H-fb",0,hmi.hmicol); 00325 00326 /* >>chng 96 nov 15, fac of 2 in deriv to help convergence in very dense 00327 * models where H- is important, this takes change in eden into 00328 * partial account */ 00329 thermal.dCooldT += 2.*hmi.hmicol*phycon.teinv; 00330 if( PRT_DERIV ) 00331 fprintf(ioQQQ,"DEBUG dcdt 2 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00332 00333 CoolAdd("H2ln",0,CoolHeavy.h2line); 00334 /* >>chng 00 oct 21, added coef of 3.5, sign had been wrong */ 00335 /*thermal.dCooldT += CoolHeavy.h2line*phycon.teinv;*/ 00336 /* >>chng 03 mar 17, change 3.5 to 15 as per behavior in primal.in */ 00337 /*thermal.dCooldT += 3.5*CoolHeavy.h2line*phycon.teinv;*/ 00338 /* >>chng 03 may 18, from 15 to 30 as per behavior in primal.in - overshoots happen */ 00339 /*thermal.dCooldT += 15.*CoolHeavy.h2line*phycon.teinv;*/ 00340 /*>>chng 03 oct 03, from 30 to 3, no more overshoots in primalin */ 00341 /*thermal.dCooldT += 30.*CoolHeavy.h2line*phycon.teinv;*/ 00342 thermal.dCooldT += 3.0*CoolHeavy.h2line*phycon.teinv; 00343 00344 { 00345 /* problems with H2 cooling */ 00346 enum {DEBUG_LOC=false}; 00347 if( DEBUG_LOC /*&& nzone>300 && iteration > 1*/) 00348 { 00349 fprintf(ioQQQ,"CoolEvaluate debuggg\t%.2f\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n", 00350 fnzone, 00351 phycon.te, 00352 hmi.H2_total , 00353 CoolHeavy.h2line, 00354 hmi.Hmolec[ipMHm] , 00355 dense.eden); 00356 } 00357 } 00358 00359 CoolAdd("HDro",0,CoolHeavy.HD); 00360 thermal.dCooldT += CoolHeavy.HD*phycon.teinv; 00361 00362 CoolAdd("H2+ ",0,CoolHeavy.H2PlsCool); 00363 thermal.dCooldT += CoolHeavy.H2PlsCool*phycon.teinv; 00364 00365 /* cooling due to C12O16 */ 00366 CoolAdd("CO C",12,CoolHeavy.C12O16Rot); 00367 thermal.dCooldT += CoolHeavy.dC12O16Rot; 00368 /* >>chng 00 oct 25, add C13O16 cooling */ 00369 /* cooling due to C13O16 */ 00370 CoolAdd("CO C",13,CoolHeavy.C13O16Rot); 00371 thermal.dCooldT += CoolHeavy.dC13O16Rot; 00372 00373 /* heating due to three-body, will be incremented in iso_cool*/ 00374 thermal.heating[0][3] = 0.; 00375 /* heating due to hydrogen lines */ 00376 thermal.heating[0][23] = 0.; 00377 /* heating due to photoionization of all excited states of hydrogen species */ 00378 thermal.heating[0][1] = 0.; 00379 00380 /* isoelectronic species cooling, mainly lines, and level ionization */ 00381 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00382 { 00383 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00384 { 00385 /* must always call iso_cool since we must zero those variables 00386 * that would have been set had the species been present */ 00387 iso_cool( ipISO , nelem ); 00388 } 00389 } 00390 00391 /* >>chng 02 jun 18, don't reevaluate needlessly */ 00392 /* >>chng 03 nov 28, even faster - special logic for when ff is pretty 00393 * negligible - eval of ff is pretty slow */ 00394 /* >>chng 04 feb 19, must not test on temp not changing, since ionization 00395 * can change at constant temperature 00396 * >>chng 04 sep 14, above introduced bug since brems never reevaluated 00397 * now test is zone or temp has changed */ 00398 if( (fabs(CoolHeavy.brems_cool_net)/SDIV(thermal.ctot)>0.001 ) || 00399 conv.lgSearch || 00400 fabs((phycon.te-TeUsedBrems)/phycon.te-1.)>0.01 || 00401 nzone!=nzoneUsedBrems ) 00402 { 00403 double BremsThisEnergy; 00404 /*double OpacityThisIon;*/ 00405 long int limit; 00406 /* free-free free free brems emission for all ions */ 00407 00408 TeUsedBrems = phycon.te; 00409 nzoneUsedBrems = nzone; 00410 /* highest frequency where we have non-zero Boltzmann factors */ 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 realnum sumion[LIMELM+1]; 00422 long int ion_lo , ion_hi; 00423 00424 ASSERT(rfield.ipEnergyBremsThin < rfield.nupper); 00425 ASSERT(limit < rfield.nupper); 00426 00427 /* ipEnergyBremsThin is index to energy where gas becomes optically thin to brems, 00428 * so this loop is over optically thin frequencies 00429 * do not include optically thick part as net emission since self absorbed */ 00430 00431 /* do hydrogen first, before main loop since want to break out as separate 00432 * coolant, and what to add on H- brems */ 00433 nelem = ipHYDROGEN; 00434 ion = 1; 00435 CoolHeavy.brems_cool_h = 0.; 00436 CoolHeavy.brems_cool_hminus = 0.; 00437 /* this is all done in opacity_addtotal - why do here too? */ 00438 for( i=rfield.ipEnergyBremsThin; i < limit; i++ ) 00439 { 00440 00441 /* in all following CoolHeavy.lgFreeOn is flag set with 'no free-free' to 00442 * turn off brems heating and cooling */ 00443 BremsThisEnergy = rfield.gff[ion][i]*rfield.widflx[i]*rfield.ContBoltz[i]; 00444 /*ASSERT( BremsThisEnergy >= 0. );*/ 00445 CoolHeavy.brems_cool_h += BremsThisEnergy; 00446 00447 /* for case of hydrogen, do H- brems - OpacStack contains the ratio 00448 * of the H- to H brems cross section - multiply by this and the fraction in ground */ 00449 CoolHeavy.brems_cool_hminus += BremsThisEnergy * opac.OpacStack[i-1+opac.iphmra]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop; 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 /* now do helium, both He+ and He++ */ 00457 nelem = ipHELIUM; 00458 CoolHeavy.brems_cool_he = 0.; 00459 for( i=rfield.ipEnergyBremsThin; i < limit; i++ ) 00460 { 00461 /* eff. charge is ion, so first rfield.gff argument must be "ion". */ 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 /* >>chng 05 jul 13, rewrite this for speed */ 00471 /* gaunt factors depend only on photon energy and ion charge, so do 00472 * sum of ions here before entering into loop over photon energy */ 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 /* now include the charge, density, and temperature */ 00485 sumion[ion] *= POW2((realnum)ion); 00486 } 00487 /* now find lowest and highest ion we need to consider - following loop is over 00488 * full continuum and eats time 00489 * >>chng 05 oct 19, bounds check had been on ion, rather than ion_lo and ion_hi, so 00490 * array bounds were exceeded */ 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 /* heavy elements */ 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 /* the total heating due to bremsstrahlung */ 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 enum {DEBUG_LOC=false}; 00516 if( DEBUG_LOC && nzone>60 /*&& iteration > 1*/) 00517 { 00518 double sumfield = 0., sumtot=0., sum1=0., sum2=0.; 00519 for( i=rfield.ipEnergyBremsThin; i<limit; i++ ) 00520 { 00521 sumtot += opac.FreeFreeOpacity[i]*rfield.flux[i]*rfield.anu[i]; 00522 sumfield += rfield.flux[i]*rfield.anu[i]; 00523 sum1 += opac.FreeFreeOpacity[i]*rfield.flux[i]*rfield.anu[i]; 00524 sum2 += opac.FreeFreeOpacity[i]*rfield.flux[i]; 00525 } 00526 fprintf(ioQQQ,"DEBUG brems heat\t%.2f\t%.3e\t%.3e\t%.3e\t%e\t%.3e\t%.3e\n", 00527 fnzone, 00528 CoolHeavy.brems_heat_total, 00529 sumtot/SDIV(sumfield) , 00530 sum1/SDIV(sum2), 00531 phycon.te , 00532 rfield.gff[1][1218], 00533 opac.FreeFreeOpacity[1218]); 00534 } 00535 } 00536 } 00537 } 00538 00539 00540 /* these two terms are both large, nearly canceling, near lte */ 00541 CoolHeavy.brems_cool_net = 00542 CoolHeavy.brems_cool_h + 00543 CoolHeavy.brems_cool_he + 00544 CoolHeavy.brems_cool_hminus + 00545 CoolHeavy.brems_cool_metals - 00546 CoolHeavy.brems_heat_total; 00547 /*fprintf(ioQQQ,"DEBUG brems\t%.2f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n", 00548 fnzone, 00549 phycon.te, 00550 CoolHeavy.brems_cool_net, 00551 CoolHeavy.brems_cool_h , 00552 CoolHeavy.brems_cool_he , 00553 CoolHeavy.brems_cool_hminus, 00554 CoolHeavy.brems_cool_metals , 00555 CoolHeavy.brems_heat_total);*/ 00556 00557 /* net free free brems cooling, count as cooling if positive */ 00558 CoolAdd( "FF c" , 0, MAX2(0.,CoolHeavy.brems_cool_net) ); 00559 00560 /* now stuff into heating array if negative */ 00561 thermal.heating[0][11] = MAX2(0.,-CoolHeavy.brems_cool_net); 00562 00563 /* >>chng 96 oct 30, from HFFNet to just FreeFreeCool, 00564 * since HeatSum picks up CoolHeavy.brems_heat_total */ 00565 thermal.dCooldT += CoolHeavy.brems_cool_h*thermal.halfte; 00566 if( PRT_DERIV ) 00567 fprintf(ioQQQ,"DEBUG dcdt 3 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00568 00569 /* >>chng 02 jun 21, net cooling already includes this */ 00570 /* end of brems cooling */ 00571 00572 /* heavy element recombination cooling, do not count hydrogenic since 00573 * already done above, also helium singlets have been done */ 00574 /* >>chng 02 jul 21, put in charge dependent rec term */ 00575 CoolHeavy.heavfb = 0.; 00576 for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ ) 00577 { 00578 if( dense.lgElmtOn[nelem] ) 00579 { 00580 /*>>chng 02 jul 21, upper bound now uses NISO, 00581 * note that detailed iso seq atoms are done in iso_cool */ 00582 /* >>chng 03 sep 03, do not do zeroed stages of ionization */ 00583 long limit_lo = MAX2( 1 , dense.IonLow[nelem] ); 00584 /* there is a -1 on the NISO since loop will be up to <= not < */ 00585 long limit_hi = MIN2( nelem-NISO+1-1 , dense.IonHigh[nelem] ); 00586 for( ion=limit_lo; ion<=limit_hi; ++ion ) 00587 /*for( ion=1; ion < nelem-NISO+1; ion++ )*/ 00588 { 00589 /* factor of 0.9 is roughly correct for nebular conditions, see 00590 * >>refer H rec cooling LaMothe, J., & Ferland, G.J., 2001, PASP, 113, 165 */ 00591 /* note that ionbal.RR_rate_coef_used is the rate coef, cm3 s-1, needs eden */ 00592 /* >>chng 02 nov 07, move rec arrays around, this now has ONLY rad rec, 00593 * previously had included grain rec and three body */ 00594 /* recombination cooling for iso-seq atoms are done in iso_cool */ 00595 double one = dense.xIonDense[nelem][ion] * ionbal.RR_rate_coef_used[nelem][ion-1]* 00596 dense.eden * phycon.te * BOLTZMANN; 00597 /*fprintf(ioQQQ,"debugggfb\t%li\t%li\t%.3e\t%.3e\t%.3e\n", nelem, ion, one, 00598 dense.xIonDense[nelem][ion] , ionbal.RR_rate_coef_used[nelem][ion]);*/ 00599 CoolHeavy.heavfb += one; 00600 } 00601 } 00602 } 00603 00604 /*fprintf(ioQQQ,"debuggg hvFB\t%i\t%.2f\t%.2e\t%.2e\n",iteration, fnzone,CoolHeavy.heavfb, dense.eden);*/ 00605 00606 CoolAdd("hvFB",0,CoolHeavy.heavfb); 00607 thermal.dCooldT += CoolHeavy.heavfb*.113*phycon.teinv; 00608 if( PRT_DERIV ) 00609 fprintf(ioQQQ,"DEBUG dcdt 4 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00610 00611 /* electron-electron brems, approx form from 00612 * >>refer ee brems Stepney and Guilbert, MNRAS 204, 1269 (1983) 00613 * ok for T<10**9 */ 00614 CoolHeavy.eebrm = POW2(dense.eden*phycon.te*1.84e-21); 00615 00616 /* >>chng 97 mar 12, added deriv */ 00617 thermal.dCooldT += CoolHeavy.eebrm*thermal.halfte; 00618 CoolAdd("eeff",0,CoolHeavy.eebrm); 00619 00620 /* add advective heating and cooling */ 00621 /* this is cooling due to loss of matter from this region */ 00622 CoolAdd("adve",0,dynamics.Cool ); 00623 /* >>chng02 dec 04, rm factor of 8 in front of dCooldT */ 00624 thermal.dCooldT += dynamics.dCooldT; 00625 /* local heating due to matter moving into this location */ 00626 thermal.heating[1][5] = dynamics.Heat; 00627 thermal.dHeatdT += dynamics.dHeatdT; 00628 00629 /* total Compton cooling */ 00630 CoolHeavy.tccool = rfield.cmcool*phycon.te; 00631 CoolAdd("Comp",0,CoolHeavy.tccool); 00632 thermal.dCooldT += rfield.cmcool; 00633 00634 /* option for "extra" cooling, expressed as power-law in temperature, these 00635 * are set with the CEXTRA command */ 00636 if( thermal.lgCExtraOn ) 00637 { 00638 CoolHeavy.cextxx = 00639 (realnum)(thermal.CoolExtra*pow(phycon.te/1e4,(double)thermal.cextpw)); 00640 } 00641 else 00642 { 00643 CoolHeavy.cextxx = 0.; 00644 } 00645 CoolAdd("Extr",0,CoolHeavy.cextxx); 00646 00647 /* cooling due to wind expansion, only for winds expansion cooling */ 00648 if( wind.windv > 0. ) 00649 { 00650 dynamics.dDensityDT = (realnum)(wind.AccelTot/wind.windv + 2.*wind.windv/ 00651 radius.Radius); 00652 CoolHeavy.expans = dense.pden*phycon.te*BOLTZMANN*dynamics.dDensityDT; 00653 } 00654 else 00655 { 00656 dynamics.dDensityDT = 0.; 00657 CoolHeavy.expans = 0.; 00658 } 00659 CoolAdd("Expn",0,CoolHeavy.expans); 00660 thermal.dCooldT += CoolHeavy.expans/phycon.te; 00661 00662 /* cyclotron cooling */ 00663 /* coef is 4/3 /8pi /c * vtherm(elec) */ 00664 CoolHeavy.cyntrn = 4.5433e-25f*magnetic.pressure*PI8*dense.eden*phycon.te; 00665 CoolAdd("Cycl",0,CoolHeavy.cyntrn); 00666 thermal.dCooldT += CoolHeavy.cyntrn/phycon.te; 00667 00668 /* heavy element collisional ionization 00669 * derivative should be zero since increased coll ion rate 00670 * decreases neutral fraction by proportional amount */ 00671 CoolAdd("Hvin",0,CoolHeavy.colmet); 00672 00673 /* evaluate H 21 cm spin changing collisions */ 00674 if( fabs(phycon.te-TeEvalCS_21cm)/phycon.te > 0.001 ) 00675 { 00676 { 00677 /* this prints table of rates at points given in original data paper */ 00678 enum {DEBUG_LOC=false}; 00679 if( DEBUG_LOC ) 00680 { 00681 # define N21CM_TE 16 00682 int n; 00683 double teval[N21CM_TE]={2.,5.,10.,20.,50.,100.,200.,500.,1000., 00684 2000.,3000.,5000.,7000.,10000.,15000.,20000.}; 00685 for( n = 0; n<N21CM_TE; ++n ) 00686 { 00687 fprintf( 00688 ioQQQ,"DEBUG 21 cm deex Te=\t%.2e\tH0=\t%.2e\tp=\t%.2e\te=\t%.2e\n", 00689 teval[n], 00690 H21cm_H_atom( teval[n] ), 00691 H21cm_proton( teval[n] ), 00692 H21cm_electron( teval[n] ) ); 00693 } 00694 cdEXIT(EXIT_FAILURE); 00695 # undef N21CM_TE 00696 } 00697 } 00698 /*only evaluate T dependent part when Te changes, but update 00699 * density part below since densities may constantly change */ 00700 atomic_rate_21cm = H21cm_H_atom( phycon.te ); 00701 proton_rate_21cm = H21cm_proton( phycon.te ); 00702 electron_rate_21cm = H21cm_electron( phycon.te ); 00703 TeEvalCS_21cm = phycon.te; 00704 } 00705 /* H 21 cm emission/population, 00706 * cs will be sum of e cs and H cs converted from rate */ 00707 cs = (electron_rate_21cm * dense.eden + 00708 atomic_rate_21cm * dense.xIonDense[ipHYDROGEN][0] + 00709 proton_rate_21cm * dense.xIonDense[ipHYDROGEN][1] ) * 00710 3./ dense.cdsqte; 00711 PutCS( cs , &HFLines[0] ); 00712 00713 /* fine structure lines */ 00714 if( fabs(phycon.te-TeEvalCS)/phycon.te > 0.05 ) 00715 { 00716 /* H 21 cm done above, now loop over remaining lines to get collision strengths */ 00717 for( i=1; i < nHFLines; i++ ) 00718 { 00719 cs = HyperfineCS( i ); 00720 /* now generate the collision strength and put into the line array */ 00721 PutCS( cs , &HFLines[i] ); 00722 } 00723 TeEvalCS = phycon.te; 00724 } 00725 00726 /* do level pops for H 21 cm which is a special case since Lya pumping in included */ 00727 H21_cm_pops(); 00728 if( PRT_DERIV ) 00729 fprintf(ioQQQ,"DEBUG dcdt 5 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00730 00731 /* find total cooling due to hyperfine structure lines */ 00732 hyperfine.cooling_total = HFLines[0].Coll.cool; 00733 00734 /* now do level pops for all except 21 cm */ 00735 for( i=1; i < nHFLines; i++ ) 00736 { 00737 /* remember current gas-phase abundances */ 00738 realnum save = dense.xIonDense[HFLines[i].Hi->nelem-1][HFLines[i].Hi->IonStg-1]; 00739 00740 /* bail if no abundance */ 00741 if( save<=0. ) continue; 00742 00743 /* set gas-phase abundance to total times isotope ratio */ 00744 dense.xIonDense[HFLines[i].Hi->nelem-1][HFLines[i].Hi->IonStg-1] *= 00745 hyperfine.HFLabundance[i]; 00746 00747 /* use the collision strength generated above and find pops and cooling */ 00748 atom_level2( &HFLines[i] ); 00749 00750 /* put the correct gas-phase abundance back in the array */ 00751 dense.xIonDense[HFLines[i].Hi->nelem-1][HFLines[i].Hi->IonStg-1] = save; 00752 00753 /* find total cooling due to hyperfine structure lines */ 00754 hyperfine.cooling_total += HFLines[i].Coll.cool; 00755 } 00756 if( PRT_DERIV ) 00757 fprintf(ioQQQ,"DEBUG dcdt 6 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00758 00759 /* Carbon cooling */ 00760 CoolCarb(); 00761 if( PRT_DERIV ) 00762 fprintf(ioQQQ,"DEBUG dcdt C %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00763 00764 /* Nitrogen cooling */ 00765 CoolNitr(); 00766 if( PRT_DERIV ) 00767 fprintf(ioQQQ,"DEBUG dcdt N %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00768 00769 /* Oxygen cooling */ 00770 CoolOxyg(); 00771 if( PRT_DERIV ) 00772 fprintf(ioQQQ,"DEBUG dcdt 7 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00773 00774 /* Fluorine cooling */ 00775 CoolFluo(); 00776 00777 /* Neon cooling */ 00778 CoolNeon(); 00779 if( PRT_DERIV ) 00780 fprintf(ioQQQ,"DEBUG dcdt Ne %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00781 00782 /* Magnesium cooling */ 00783 CoolMagn(); 00784 if( PRT_DERIV ) 00785 fprintf(ioQQQ,"DEBUG dcdt 8 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00786 00787 /* Sodium cooling */ 00788 CoolSodi(); 00789 if( PRT_DERIV ) 00790 fprintf(ioQQQ,"DEBUG dcdt Na %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00791 00792 /* Aluminum cooling */ 00793 CoolAlum(); 00794 if( PRT_DERIV ) 00795 fprintf(ioQQQ,"DEBUG dcdt Al %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00796 00797 /* Silicon cooling */ 00798 CoolSili(); 00799 if( PRT_DERIV ) 00800 fprintf(ioQQQ,"DEBUG dcdt 9 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00801 00802 /* Phosphorus */ 00803 CoolPhos(); 00804 00805 /* Sulphur cooling */ 00806 CoolSulf(); 00807 00808 /* Chlorine cooling */ 00809 CoolChlo(); 00810 00811 /* Argon cooling */ 00812 CoolArgo(); 00813 if( PRT_DERIV ) 00814 fprintf(ioQQQ,"DEBUG dcdt 10 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00815 00816 /* Potasium cooling */ 00817 CoolPota(); 00818 00819 /* Calcium cooling */ 00820 CoolCalc(); 00821 00822 /* Scandium cooling */ 00823 CoolScan(); 00824 00825 /* Titanium cooling */ 00826 CoolTita(); 00827 if( PRT_DERIV ) 00828 fprintf(ioQQQ,"DEBUG dcdt 11 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00829 00830 /* Vanadium cooling */ 00831 CoolVana(); 00832 00833 /* Chromium cooling */ 00834 CoolChro(); 00835 00836 /* Iron cooling */ 00837 CoolIron(); 00838 if( PRT_DERIV ) 00839 fprintf(ioQQQ,"DEBUG dcdt 12 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00840 00841 /* Manganese cooling */ 00842 CoolMang(); 00843 00844 /* Cobalt cooling */ 00845 CoolCoba(); 00846 00847 /* Nickel cooling */ 00848 CoolNick(); 00849 00850 /* Zinc cooling */ 00851 CoolZinc(); 00852 if( PRT_DERIV ) 00853 fprintf(ioQQQ,"DEBUG dcdt 13 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00854 00855 /* do all the thousands of lines Dima added with g-bar approximation */ 00856 CoolDima(); 00857 00858 /* do Chianti & Leiden database atoms and molecules here */ 00859 atmol_popsolve(); 00860 00861 /* now add up all the coolants */ 00862 CoolSum(tot); 00863 if( PRT_DERIV ) 00864 fprintf(ioQQQ,"DEBUG dcdt 14 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT); 00865 00866 /* negative cooling */ 00867 if( *tot <= 0. ) 00868 { 00869 fprintf( ioQQQ, " COOLR; cooling is <=0, this is impossible.\n" ); 00870 ShowMe(); 00871 cdEXIT(EXIT_FAILURE); 00872 } 00873 00874 /* bad derivative */ 00875 if( thermal.dCooldT == 0. ) 00876 { 00877 fprintf( ioQQQ, " COOLR; cooling slope <=0, this is impossible.\n" ); 00878 if( *tot > 0. && dense.gas_phase[ipHYDROGEN] < 1e-4 ) 00879 { 00880 fprintf( ioQQQ, " Probably due to very low density.\n" ); 00881 } 00882 ShowMe(); 00883 cdEXIT(EXIT_FAILURE); 00884 } 00885 00886 if( trace.lgTrace ) 00887 { 00888 fndstr(*tot,thermal.dCooldT); 00889 } 00890 00891 /* lgTSetOn true for constant temperature model */ 00892 if( (((((!thermal.lgTemperatureConstant) && *tot < 0.) && called.lgTalk) && 00893 !conv.lgSearch) && thermal.lgCNegChk) && nzone > 0 ) 00894 { 00895 fprintf( ioQQQ, 00896 " NOTE Negative cooling, zone %4ld, =%10.2e coola=%10.2e CHION=%10.2e Te=%10.2e\n", 00897 nzone, 00898 *tot, 00899 iso.cLya_cool[ipH_LIKE][ipHYDROGEN], 00900 iso.coll_ion[ipH_LIKE][ipHYDROGEN], 00901 phycon.te ); 00902 fndneg(); 00903 } 00904 00905 /* possibility of getting empirical cooling derivative 00906 * normally false, set true with 'set numerical derivatives' command */ 00907 if( NumDeriv.lgNumDeriv ) 00908 { 00909 if( ((nzone > 2 && nzone == nzSave) && ! fp_equal( oldtemp, phycon.te )) && nhit > 4 ) 00910 { 00911 /* hnit is number of tries on this zone - use to stop numerical problems 00912 * do not evaluate numerical deriv until well into solution */ 00913 deriv = (oltcool - *tot)/(oldtemp - phycon.te); 00914 thermal.dCooldT = deriv; 00915 } 00916 else 00917 { 00918 deriv = thermal.dCooldT; 00919 } 00920 if( nzone != nzSave ) 00921 nhit = 0; 00922 00923 nzSave = nzone; 00924 nhit += 1; 00925 oltcool = *tot; 00926 oldtemp = phycon.te; 00927 } 00928 return; 00929 } 00930 00931 /* */ 00932 #ifdef EPS 00933 # undef EPS 00934 #endif 00935 #define EPS 0.01 00936 00937 /*fndneg search cooling array to find negative values */ 00938 STATIC void fndneg(void) 00939 { 00940 long int i; 00941 double trig; 00942 00943 DEBUG_ENTRY( "fndneg()" ); 00944 00945 trig = fabs(thermal.htot*EPS); 00946 for( i=0; i < thermal.ncltot; i++ ) 00947 { 00948 if( thermal.cooling[i] < 0. && fabs(thermal.cooling[i]) > trig ) 00949 { 00950 fprintf( ioQQQ, " negative line=%s %.2f fraction of heating=%.3f\n", 00951 thermal.chClntLab[i], thermal.collam[i], thermal.cooling[i]/ 00952 thermal.htot ); 00953 } 00954 00955 if( thermal.heatnt[i] > trig ) 00956 { 00957 fprintf( ioQQQ, " heating line=%s %.2f fraction of heating=%.3f\n", 00958 thermal.chClntLab[i], thermal.collam[i], thermal.heatnt[i]/ 00959 thermal.htot ); 00960 } 00961 } 00962 return; 00963 } 00964 00965 /*fndstr search cooling stack to find strongest values */ 00966 STATIC void fndstr(double tot, 00967 double dc) 00968 { 00969 char chStrngLab[NCOLNT_LAB_LEN+1]; 00970 long int i; 00971 realnum wl; 00972 double str, 00973 strong; 00974 00975 DEBUG_ENTRY( "fndstr()" ); 00976 00977 strong = 0.; 00978 wl = -FLT_MAX; 00979 for( i=0; i < thermal.ncltot; i++ ) 00980 { 00981 if( fabs(thermal.cooling[i]) > strong ) 00982 { 00983 /* this is the wavelength of the coolant, 0 for a continuum*/ 00984 wl = thermal.collam[i]; 00985 /* make sure labels are all valid*/ 00986 /*>>chng 06 jun 06, bug fix, assert length was ==4, should be <=NCOLNT_LAB_LEN */ 00987 ASSERT( strlen( thermal.chClntLab[i] ) <= NCOLNT_LAB_LEN ); 00988 strcpy( chStrngLab, thermal.chClntLab[i] ); 00989 strong = fabs(thermal.cooling[i]); 00990 } 00991 } 00992 00993 str = strong; 00994 00995 fprintf( ioQQQ, 00996 " fndstr cool: TE=%10.4e Ne=%10.4e C=%10.3e dC/dT=%10.2e ABS(%s %.1f)=%.2e nz=%ld\n", 00997 phycon.te, dense.eden, tot, dc, chStrngLab 00998 , wl, str, nzone ); 00999 01000 /* option for extensive printout of lines */ 01001 if( trace.lgCoolTr ) 01002 { 01003 realnum ratio; 01004 01005 /* flag all significant coolants, first zero out the array */ 01006 coolpr(ioQQQ,(char*)thermal.chClntLab[0],1,0.,"ZERO"); 01007 01008 /* push all coolants onto the stack */ 01009 for( i=0; i < thermal.ncltot; i++ ) 01010 { 01011 /* usually positive, although can be neg for coolants that heats, 01012 * only do positive here */ 01013 ratio = (realnum)(thermal.cooling[i]/thermal.ctot); 01014 if( ratio >= EPS ) 01015 { 01016 /*>>chng 99 jan 27, only cal when ratio is significant */ 01017 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i], ratio,"DOIT"); 01018 } 01019 } 01020 01021 /* complete the printout for positive coolants */ 01022 coolpr(ioQQQ,"DONE",1,0.,"DONE"); 01023 01024 /* now do cooling that was counted as a heat source if significant */ 01025 if( thermal.heating[0][22]/thermal.ctot > 0.05 ) 01026 { 01027 fprintf( ioQQQ, 01028 " All coolant heat greater than%6.2f%% of the total will be printed.\n", 01029 EPS*100. ); 01030 01031 coolpr(ioQQQ,"ZERO",1,0.,"ZERO"); 01032 for( i=0; i < thermal.ncltot; i++ ) 01033 { 01034 ratio = (realnum)(thermal.heatnt[i]/thermal.ctot); 01035 if( fabs(ratio) >=EPS ) 01036 { 01037 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i], 01038 ratio,"DOIT"); 01039 } 01040 } 01041 coolpr(ioQQQ,"DONE",1,0.,"DONE"); 01042 } 01043 } 01044 return; 01045 }