cloudy: /home66/gary/public_html/cloudy/c13_branch/source/cool

00001 /* This file is part of Cloudy and is copyright (C)1978-2013 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 "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 /*fndneg search cooling array to find negative values */
00037 STATIC void fndneg(void);
00038 /*fndstr search cooling stack to find strongest values */
00039 STATIC void fndstr(double tot, 
00040   double dc);
00041 
00042 /* set true to debug derivative of heating and cooling */
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         /* returns tot, the total cooling,
00075          * and dc, the derivative of the cooling */
00076 
00077         /* routine atom_level2( t10 )
00078          * routine atom_level3( abund , t10,t21,t20)
00079          * tsq1 = 1. / (te**2)
00080          * POPEXC( O12,g1,g2,A21,excit,abund); result already*a21
00081          * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2)
00082          * AtomSeqBeryllium(cs23,cs24,cs34,tarray,a41)
00083          * FIVEL( G(1-5) , ex(wn,1-5), cs12,cs13,14,15,23,24,25,34,35,45,
00084          *  A21,31,41,51,32,42,52,43,53,54, pop(1-5), abund) */
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         /* must call TempChange since ionization has changed, there are some
00093          * terms that affect collision rates (H0 term in electron collision) */
00094         TempChange(phycon.te , false);
00095 
00096         /* now zero out the cooling stack */
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                 /* grain heating and cooling */
00103                 /* grain recombination cooling, evaluated elsewhere
00104                 * can either heat or cool the gas, do cooling here */
00105                 CoolAdd("dust",0,MAX2(0.,gv.GasCoolColl));
00106 
00107                 /* grain cooling proportional to temperature ^3/2 */
00108                 thermal.dCooldT += MAX2(0.,gv.GasCoolColl)*3./(2.*phycon.te);
00109 
00110                 /* these are the various heat agents from grains */
00111                 /* options to force gas heating or cooling by grains to zero - for tests only ! */
00112                 if( gv.lgDustOn() && gv.lgDHetOn )
00113                 {
00114                         /* rate dust heats gas by photoelectric effect */
00115                         thermal.heating[0][13] = gv.GasHeatPhotoEl;
00116 
00117                         /* if grains hotter than gas then collisions with gas act
00118                         * to heat the gas, add this in here
00119                         * a symmetric statement appears in COOLR, where cooling is added on */
00120                         thermal.heating[0][14] = MAX2(0.,-gv.GasCoolColl);
00121 
00122                         /* this is gas heating due to thermionic emissions */
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                 /* >>chng 06 jul 21, option to include Bakes PAH hack with grain physics off,
00135                  * needed to test dynamics models */
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         /* molecular molecules molecule cooling */
00143         if( mole_global.lgNoMole )
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                 hmi.deriv_HeatH2Dexc_used = 0.;
00163         }
00164 
00165         else
00166         {
00167                 /* save various molecular heating/cooling agent */
00168                 thermal.heating[0][15] = hmi.hmihet;
00169                 thermal.heating[0][16] = hmi.h2plus_heat;
00170                 /* now get heating from H2 molecule, either simple or from big one */
00171                 if( h2.lgEnabled  && hmi.lgH2_Thermal_BigH2 )
00172                 {
00173                         if( h2.lgEvaluated )
00174                         {
00175                                 /* these are explicitly from big H2 molecule,
00176                                  * first is heating due to radiative pump of excited states, followed by
00177                                  * radiative decay into continuum of X, followed by dissociation of molecule
00178                                  * with kinetic energy, typically 0.25 - 0.5 eV per event */
00179                                 hmi.HeatH2Dish_used = h2.HeatDiss;
00180                                 hmi.HeatH2Dexc_used = h2.HeatDexc;
00181                                 /*fprintf(ioQQQ,"DEBUG big %.2f\t%.5e\t%.2e\t%.2e\t%.2e\n", 
00182                                         fnzone , phycon.te, hmi.HeatH2Dexc_used,
00183                                         hmi.H2_total , hmi.H2_star_BigH2);*/
00184                                 /* negative sign because right term is really deriv of heating,
00185                                  * but will be used below as deriv of cooling */
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                         /* TH85 dissociation heating */
00199                         /* these come from approximations in TH85, see comments above */
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                         /* Burton et al. 1990 */
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                         /* Bertoldi & Draine */
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                         /* this is the default when small H2 used */
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                 /* heating due to photodissociation heating */
00229                 thermal.heating[0][17] = hmi.HeatH2Dish_used;
00230 
00231                 /* heating due to continuum photodissociation */
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                 /* heating (usually cooling in big H2) due to collisions within X */
00240                 /* add to heating is net heating is positive */
00241                 thermal.heating[0][8] = MAX2(0.,hmi.HeatH2Dexc_used);
00242 
00243                 /* add to cooling if net heating is negative */
00244                 CoolAdd("H2cX",0,MAX2(0.,-hmi.HeatH2Dexc_used));
00245                 /*fprintf(ioQQQ,"DEBUG coolh2\t%.2f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
00246                         fnzone, phycon.te, dense.eden, hmi.H2_total, thermal.ctot, -hmi.HeatH2Dexc_used );*/
00247                 /* add to net derivative */
00248                 /*thermal.dCooldT += MAX2(0.,-hmi.HeatH2Dexc_used)* ( 30172. * thermal.tsq1 - thermal.halfte );*/
00249                 /* >>chng 04 jan 25, check sign to prevent cooling from entering here,
00250                  * also enter neg sign since going into cooling stack (bug), in heatsum
00251                  * same term adds to deriv of heating */
00252                 if( hmi.HeatH2Dexc_used < 0. )
00253                         thermal.dCooldT -= hmi.deriv_HeatH2Dexc_used;
00254 
00255                 /*  H + H+ => H2+ cooling */
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                         /* this is simplified approximation to H2 rotation cooling,
00262                          * big molecule goes this far better */
00263                         CoolHeavy.h2line = 0.;
00264                 }
00265                 else
00266                 {
00267                         /*  rate for rotation lines from 
00268                         *  >>refer      h2      cool    Lepp, S., & Shull, J.M. 1983, ApJ, 270, 578 */
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                         /*  low density rotation cooling */
00280                         /*&qn = pow(MAX2(findspecieslocal("H2")->den,1e-37),0.77) + 1.2*pow(MAX2(dense.xIonDense[ipHYDROGEN][0],1e-37),0.77);*/
00281                         qn = pow(MAX2(hmi.H2_total,1e-37),0.77) + 1.2*pow(MAX2(dense.xIonDense[ipHYDROGEN][0],1e-37),0.77);
00282                         /* these are equations 11 from LS83 */
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                         /* \todo 1 add this from LS83 or (better yet) update to another form.  See Galli & Palla 1998, A5-7. */
00296                         //if( viblow > 0. )
00297                         //      CoolHeavy.h2line += hmi.H2_total*vibhi/(1. + vibhi/viblow);
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                         /* >>chng 02 mar 07, add DH cooling using rates (eqn 6) from
00322                          * >>refer      HD      cooling Puy, D., Grenacher, L, & Jetzer, P., 1999, A&A, 345, 723 */
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(); // test and enable this by default
00330 #if 0
00331         double chemical_heating = mole.chem_heat();     
00332         thermal.heating[0][29] = MAX2(0.,chemical_heating);
00333         /* add to cooling if net heating is negative */
00334         CoolAdd("Chem",0,MAX2(0.,-chemical_heating));
00335 #endif
00336 
00337         /* cooling due to charge transfer ionization / recombination */
00338         CoolAdd("CT C" , 0. , thermal.char_tran_cool );
00339 
00340         /*  H- FB; H + e -> H- + hnu */
00341         /*  H- FF is in with H ff */
00342         CoolAdd("H-fb",0,hmi.hmicol);
00343 
00344         /* >>chng 96 nov 15, fac of 2 in deriv to help convergence in very dense
00345          * models where H- is important, this takes change in eden into
00346          * partial account */
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         /* >>chng 00 oct 21, added coef of 3.5, sign had been wrong */
00353         /*thermal.dCooldT += CoolHeavy.h2line*phycon.teinv;*/
00354         /* >>chng 03 mar 17, change 3.5 to 15 as per behavior in primal.in */
00355         /*thermal.dCooldT += 3.5*CoolHeavy.h2line*phycon.teinv;*/
00356         /* >>chng 03 may 18, from 15 to 30 as per behavior in primal.in - overshoots happen */
00357         /*thermal.dCooldT += 15.*CoolHeavy.h2line*phycon.teinv;*/
00358         /*>>chng 03 oct 03, from 30 to 3, no more overshoots in primalin */
00359         /*thermal.dCooldT += 30.*CoolHeavy.h2line*phycon.teinv;*/
00360         thermal.dCooldT += 3.0*CoolHeavy.h2line*phycon.teinv;
00361 
00362         {
00363                 /* problems with H2 cooling */
00364                 enum {DEBUG_LOC=false};
00365                 if( DEBUG_LOC /*&& nzone>300 && iteration > 1*/)
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         /* heating due to three-body, will be incremented in iso_cool*/
00384         thermal.heating[0][3] = 0.;
00385         /* heating due to hydrogen lines */
00386         thermal.heating[0][23] = 0.;
00387         /* heating due to photoionization of all excited states of hydrogen species */
00388         thermal.heating[0][1] = 0.;
00389 
00390         /* isoelectronic species cooling, mainly lines, and level ionization */
00391         for( long int ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00392         {
00393                 for( long int nelem=ipISO; nelem < LIMELM; nelem++ )
00394                 {
00395                         /* must always call iso_cool since we must zero those variables
00396                         * that would have been set had the species been present */
00397                         iso_cool( ipISO , nelem );
00398                 }
00399         }
00400 
00401         /* >>chng 02 jun 18, don't reevaluate needlessly */
00402         /* >>chng 03 nov 28, even faster - special logic for when ff is pretty
00403          * negligible - eval of ff is pretty slow */
00404         /* >>chng 04 feb 19, must not test on temp not changing, since ionization
00405          * can change at constant temperature 
00406          * >>chng 04 sep 14, above introduced bug since brems never reevaluated
00407          * now test is zone or temp has changed */
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                 /*double OpacityThisIon;*/
00415                 long int limit;
00416                 /* free-free free free brems emission for all ions */
00417 
00418                 TeUsedBrems = phycon.te;
00419                 nzoneUsedBrems = nzone;
00420                 /* highest frequency where we have non-zero Boltzmann factors */
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                         /* ipEnergyBremsThin is index to energy where gas becomes optically thin to brems,
00438                          * so this loop is over optically thin frequencies 
00439                          * do not include optically thick part as net emission since self absorbed */
00440 
00441                         /* do hydrogen first, before main loop since want to break out as separate
00442                          * coolant, and what to add on H- brems */
00443                         CoolHeavy.brems_cool_h = 0.;
00444                         CoolHeavy.brems_cool_hminus = 0.;
00445                         /* this is all done in opacity_addtotal - why do here too? */
00446                         for( long int i=rfield.ipEnergyBremsThin; i < limit; i++ )
00447                         {
00448                                 long int ion = 1;
00449 
00450                                 /* in all following CoolHeavy.lgFreeOn is flag set with 'no free-free' to
00451                                  * turn off brems heating and cooling */
00452                                 BremsThisEnergy = rfield.gff[ion][i]*rfield.widflx[i]*rfield.ContBoltz[i];
00453                                 /*ASSERT( BremsThisEnergy >= 0. );*/
00454                                 CoolHeavy.brems_cool_h += BremsThisEnergy;
00455 
00456                                 /* for case of hydrogen, add H- brems - OpacStack contains the ratio
00457                                  * of the H- to H brems cross section - multiply by this and H(1s) population */
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                         /* now do helium, both He+ and He++ */
00466                         CoolHeavy.brems_cool_he = 0.;
00467                         for( long int i=rfield.ipEnergyBremsThin; i < limit; i++ )
00468                         {
00469                                 long int nelem = ipHELIUM;
00470                                 /* eff. charge is ion, so first rfield.gff argument must be "ion".      */
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                         /* >>chng 05 jul 13, rewrite this for speed */
00480                         /* gaunt factors depend only on photon energy and ion charge, so do
00481                         * sum of ions here before entering into loop over photon energy */
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                                 /* now include the charge, density, and temperature */
00494                                 sumion[ion] *= POW2((realnum)ion);
00495                         }
00496 
00497                         /* add molecular ions */
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                         /* now find lowest and highest ion we need to consider - following loop is over
00515                          * full continuum and eats time
00516                          * >>chng 05 oct 19, bounds check had been on ion, rather than ion_lo and ion_hi, so
00517                          * array bounds were exceeded */
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                         /* heavy elements */
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                                 /* the total heating due to bremsstrahlung */
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 /*&& iteration > 1*/)
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         /* these two terms are both large, nearly canceling, near lte */
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         /*fprintf(ioQQQ,"DEBUG brems\t%.2f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
00574                 fnzone,
00575                 phycon.te,
00576                 CoolHeavy.brems_cool_net,
00577                 CoolHeavy.brems_cool_h ,
00578                 CoolHeavy.brems_cool_he ,
00579                 CoolHeavy.brems_cool_hminus,
00580                 CoolHeavy.brems_cool_metals ,
00581                 CoolHeavy.brems_heat_total);*/
00582 
00583         /* net free free brems cooling, count as cooling if positive */
00584         CoolAdd( "FF c" , 0, MAX2(0.,CoolHeavy.brems_cool_net) );
00585 
00586         /* now stuff into heating array if negative */
00587         thermal.heating[0][11] = MAX2(0.,-CoolHeavy.brems_cool_net);
00588 
00589         /* >>chng 96 oct 30, from HFFNet to just FreeFreeCool,
00590          * since HeatSum picks up CoolHeavy.brems_heat_total */
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         /* >>chng 02 jun 21, net cooling already includes this */
00596         /* end of brems cooling */
00597 
00598         /* heavy element recombination cooling, do not count hydrogenic since
00599          * already done above, also helium singlets have been done */
00600         /* >>chng 02 jul 21, put in charge dependent rec term */
00601         CoolHeavy.heavfb = 0.;
00602         for( long int nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
00603         {
00604                 if( dense.lgElmtOn[nelem] )
00605                 {
00606                         /* note that detailed iso seq atoms are done in iso_cool */
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                                 /* factor of 0.9 is roughly correct for nebular conditions, see
00612                                  * >>refer      H       rec cooling     LaMothe, J., & Ferland, G.J., 2001, PASP, 113, 165 */
00613                                 /* note that ionbal.RR_rate_coef_used is the rate coef, cm3 s-1, needs eden */
00614                                 /* >>chng 02 nov 07, move rec arrays around, this now has ONLY rad rec,
00615                                  * previously had included grain rec and three body */
00616                                 /* recombination cooling for iso-seq atoms are done in iso_cool */
00617                                 double one = dense.xIonDense[nelem][ion] * ionbal.RR_rate_coef_used[nelem][ion-1]*
00618                                         dense.eden * phycon.te * BOLTZMANN;
00619                                 /*fprintf(ioQQQ,"debugggfb\t%li\t%li\t%.3e\t%.3e\t%.3e\n", nelem, ion, one, 
00620                                         dense.xIonDense[nelem][ion] , ionbal.RR_rate_coef_used[nelem][ion]);*/
00621                                 CoolHeavy.heavfb += one;
00622                         }
00623                 }
00624         }
00625 
00626         /*fprintf(ioQQQ,"debuggg hvFB\t%i\t%.2f\t%.2e\t%.2e\n",iteration, fnzone,CoolHeavy.heavfb, dense.eden);*/
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         /* electron-electron brems, approx form from 
00634          * >>refer      ee      brems   Stepney and Guilbert, MNRAS 204, 1269 (1983)
00635          * ok for T<10**9 */
00636         CoolHeavy.eebrm = POW2(dense.eden*phycon.te*1.84e-21);
00637 
00638         /* >>chng 97 mar 12, added deriv */
00639         thermal.dCooldT += CoolHeavy.eebrm*thermal.halfte;
00640         CoolAdd("eeff",0,CoolHeavy.eebrm);
00641 
00642         /* add advective heating and cooling */
00643         /* this is cooling due to loss of matter from this region */
00644         CoolAdd("adve",0,dynamics.Cool() );
00645         /* >>chng02 dec 04, rm factor of 8 in front of dCooldT */
00646         thermal.dCooldT += dynamics.dCooldT();
00647         /* local heating due to matter moving into this location */
00648         thermal.heating[1][5] = dynamics.Heat();
00649         thermal.dHeatdT += dynamics.dHeatdT;
00650 
00651         /* total Compton cooling */
00652         CoolHeavy.tccool = rfield.cmcool*phycon.te;
00653         CoolAdd("Comp",0,CoolHeavy.tccool);
00654         thermal.dCooldT += rfield.cmcool;
00655 
00656         /* option for "extra" cooling, expressed as power-law in temperature, these
00657          * are set with the CEXTRA command */
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         /* cooling due to wind expansion, only for winds expansion cooling */
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                 // pdV work term
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         /* cyclotron cooling */
00697         /* coef is 4/3 /8pi /c * vtherm(elec) */
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         /* heavy element collisional ionization
00703          * derivative should be zero since increased coll ion rate
00704          * decreases neutral fraction by proportional amount */
00705         CoolAdd("Hvin",0,CoolHeavy.colmet);
00706 
00707         /* evaluate H 21 cm spin changing collisions */
00708         coolnum = thermal.ncltot;
00709         if( !fp_equal(phycon.te,TeEvalCS_21cm) )
00710         {
00711                 {
00712                         /* this prints table of rates at points given in original data paper */
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                 /*only evaluate T dependent part when Te changes, but update
00734                  * density part below since densities may constantly change */
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         /* H 21 cm emission/population,
00741         * cs will be sum of e cs and H cs converted from rate */
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         /* fine structure lines */
00749         if( !fp_equal(phycon.te,TeEvalCS) )
00750         {
00751                 /* H 21 cm done above, now loop over remaining lines to get collision strengths */
00752                 for( long int i=1; i < nHFLines; i++ )
00753                 {
00754                         cs = HyperfineCS( i );
00755                         /* now generate the collision strength and put into the line array */
00756                         PutCS(  cs , HFLines[i] );
00757                 }
00758                 TeEvalCS = phycon.te;
00759         }
00760 
00761         /* do level pops for H 21 cm which is a special case since Lya pumping in included */
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         /* find total cooling due to hyperfine structure lines */
00768         hyperfine.cooling_total = HFLines[0].Coll().cool();
00769 
00770         /* now do level pops for all except 21 cm  */
00771         for( long int i=1; i < nHFLines; i++ )
00772         {
00773                 /* remember current gas-phase abundance of this isotope */
00774                 realnum save = dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1];
00775 
00776                 /* bail if no abundance */
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                 /* set gas-phase abundance to total times isotope ratio */
00783                 dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1] *= 
00784                         hyperfine.HFLabundance[i];
00785 
00786                 /* use the collision strength generated above and find pops and cooling */
00787                 atom_level2( HFLines[i] );
00788 
00789                 /* put the correct gas-phase abundance back in the array */
00790                 dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1] = save;
00791 
00792                 /* find total cooling due to hyperfine structure lines */
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                                 // zero abundance of species if we are using Chianti for this ion
00810                                 if( dense.lgIonChiantiOn[nelem][ion] || dense.lgIonStoutOn[nelem][ion] )
00811                                         dense.xIonDense[nelem][ion] = 0.;
00812                         }
00813                 }
00814         }
00815 
00816         /* Carbon cooling */
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         /* Nitrogen cooling */
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         /* Oxygen cooling */
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         /* Neon cooling */
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         /* Magnesium cooling */
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         /* Sodium cooling */
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         /* Aluminum cooling */
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         /* Silicon cooling */
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         /* Phosphorus */
00883         coolnum = thermal.ncltot;
00884         CoolPhos();
00885         for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00886                 thermal.elementcool[ipPHOSPHORUS] += thermal.cooling[coolcal];
00887 
00888         /* Sulphur cooling */
00889         coolnum = thermal.ncltot;
00890         CoolSulf();
00891         for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00892                 thermal.elementcool[ipSULPHUR] += thermal.cooling[coolcal];
00893 
00894         /* Chlorine cooling */
00895         coolnum = thermal.ncltot;
00896         CoolChlo();
00897         for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00898                 thermal.elementcool[ipCHLORINE] += thermal.cooling[coolcal];
00899 
00900         /* Argon cooling */
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         /* Potasium cooling */
00909         coolnum = thermal.ncltot;
00910         CoolPota();
00911         for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00912                 thermal.elementcool[ipPOTASSIUM] += thermal.cooling[coolcal];
00913 
00914         /* Calcium cooling */
00915         coolnum = thermal.ncltot;
00916         CoolCalc();
00917         for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00918                 thermal.elementcool[ipCALCIUM] += thermal.cooling[coolcal];
00919 
00920         /* Scandium cooling */
00921         coolnum = thermal.ncltot;
00922         CoolScan();
00923         for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00924                 thermal.elementcool[ipSCANDIUM] += thermal.cooling[coolcal];
00925 
00926         /* Chromium cooling */
00927         coolnum = thermal.ncltot;
00928         CoolChro();
00929         for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00930                 thermal.elementcool[ipCHROMIUM] += thermal.cooling[coolcal];
00931         
00932 
00933         /* Iron cooling */
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         /* Cobalt cooling */
00942         coolnum = thermal.ncltot;
00943         CoolCoba();
00944         for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00945                 thermal.elementcool[ipCOBALT] += thermal.cooling[coolcal];
00946 
00947         /* Nickel cooling */
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         // reset abundances to original values, may have been set zero to protect against old cloudy lines
00956         if( atmdat.lgChiantiOn || atmdat.lgStoutOn)
00957         {
00958                 // this clause, first reset abundances set to zero when Chianti included
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         /* opacity project lines Dima Verner added with g-bar approximation */
00969         CoolDima();
00970 
00971         for( int coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
00972                 thermal.dima += thermal.cooling[coolcal];
00973 
00974         /* do external database lines */
00975         dBase_solve();
00976         
00977         /* Print number of levels for each species */
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         /* now add up all the coolants */
01009         CoolSum(tot);
01010         if( PRT_DERIV )
01011                 fprintf(ioQQQ,"DEBUG dCdT 14 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
01012 
01013         /* negative cooling */
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         /* bad derivative */
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         /* lgTSetOn true for constant temperature model */
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         /* possibility of getting empirical cooling derivative
01053          * normally false, set true with 'set numerical derivatives' command */
01054         if( NumDeriv.lgNumDeriv )
01055         {
01056                 if( ((nzone > 2 && nzone == nzSave) && ! fp_equal( oldtemp, phycon.te )) && nhit > 4 )
01057                 {
01058                         /* hnit is number of tries on this zone - use to stop numerical problems
01059                          * do not evaluate numerical deriv until well into solution */
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 /*fndneg search cooling array to find negative values */
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 /*fndstr search cooling stack to find strongest values */
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                         /* this is the wavelength of the coolant, 0 for a continuum*/
01131                         wl = thermal.collam[i];
01132                         /* make sure labels are all valid*/
01133                         /*>>chng 06 jun 06, bug fix, assert length was ==4, should be <=NCOLNT_LAB_LEN */
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         /* option for extensive printout of lines */
01148         if( trace.lgCoolTr )
01149         {
01150                 realnum ratio;
01151 
01152                 /* flag all significant coolants, first zero out the array */
01153                 coolpr(ioQQQ,(char*)thermal.chClntLab[0],1,0.,"ZERO");
01154 
01155                 /* push all coolants onto the stack */
01156                 for( i=0; i < thermal.ncltot; i++ )
01157                 {
01158                         /* usually positive, although can be neg for coolants that heats, 
01159                          * only do positive here */
01160                         ratio = (realnum)(thermal.cooling[i]/thermal.ctot);
01161                         if( ratio >= EPS )
01162                         {
01163                                 /*>>chng 99 jan 27, only cal when ratio is significant */
01164                                 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i], ratio,"DOIT");
01165                         }
01166                 }
01167 
01168                 /* complete the printout for positive coolants */
01169                 coolpr(ioQQQ,"DONE",1,0.,"DONE");
01170 
01171                 /* now do cooling that was counted as a heat source if significant */
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 }
/home66/gary/public_html/cloudy/c13_branch/source/cool_eval.cpp
