00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "physconst.h"
00007 #include "thermal.h"
00008 #include "heavy.h"
00009 #include "trace.h"
00010 #include "secondaries.h"
00011 #include "conv.h"
00012 #include "called.h"
00013 #include "coolheavy.h"
00014 #include "iso.h"
00015 #include "mole.h"
00016 #include "h2.h"
00017 #include "hmi.h"
00018 #include "dense.h"
00019 #include "ionbal.h"
00020 #include "phycon.h"
00021 #include "numderiv.h"
00022 #include "atomfeii.h"
00023 #include "cooling.h"
00024 #include "grainvar.h"
00025 #include "taulines.h"
00026 #include "deuterium.h"
00027
00028
00029 static const double FAINT_HEAT = 0.02;
00030
00031 static const bool PRT_DERIV = false;
00032
00033 void HeatSum( void )
00034 {
00035
00036 const int NDIM = 40;
00037
00038
00039 double cosmic_ray_ionization_rate ,
00040 pair_production_ionization_rate ,
00041 fast_neutron_ionization_rate , SecExcitLyaRate;
00042
00043
00044 double SeconIoniz_iso[NISO] ,
00045 SeconExcit_iso[NISO];
00046
00047 long int i,
00048 ion,
00049 ipnt,
00050 ipsave[NDIM],
00051 j,
00052 jpsave[NDIM],
00053 limit;
00054
00055 double HeatingLo ,
00056 HeatingHi ,
00057 secmet ,
00058 smetla;
00059 long ipISO,
00060 ns;
00061 static long int nhit=0,
00062 nzSave=0;
00063 double photo_heat_2lev_atoms,
00064 photo_heat_ISO_atoms ,
00065 photo_heat_UTA ,
00066 OtherHeat ,
00067 deriv,
00068 oldfac,
00069 save[NDIM];
00070 static double oldheat=0.,
00071 oldtemp=0.;
00072 double secmetsave[LIMELM];
00073
00074 realnum SaveOxygen1 , SaveCarbon1;
00075
00076
00077
00078 double Cosmic_ray_heat_eff ,
00079 Cosmic_ray_sec2prim;
00080 realnum sec2prim_par_1;
00081 realnum sec2prim_par_2;
00082
00083 DEBUG_ENTRY( "HeatSum()" );
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108 realnum xValenceDonorDensity, ElectronFraction;
00109
00110 if (1)
00111 {
00112 long ipNoble[] = {
00113 ipHELIUM, ipNEON, ipARGON, ipKRYPTON
00114
00115 };
00116
00117 long ipFullShell = -1, ipNextFull = 0;
00118 xValenceDonorDensity = 0.;
00119 for (long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem)
00120 {
00121
00122 xValenceDonorDensity +=
00123 (nelem-ipFullShell)*dense.gas_phase[nelem];
00124 if (nelem == ipNoble[ipNextFull])
00125 {
00126 ipFullShell = nelem;
00127 ++ipNextFull;
00128 }
00129 }
00130
00131
00132
00133 ElectronFraction = (realnum)(MAX2(MIN2(1.0,dense.eden/
00134 xValenceDonorDensity),1e-4));
00135 }
00136 else
00137 {
00138
00139
00140
00141
00143 realnum xNeutralParticleDensity = 0.;
00144 for( long nelem=0; nelem < LIMELM; nelem++ )
00145 {
00146 xNeutralParticleDensity += dense.xIonDense[nelem][0];
00147 }
00148
00149 xNeutralParticleDensity += deut.xIonDense[0];
00150
00151
00152 for( i=0; i < mole_global.num_calc; i++ )
00153 {
00154
00155 if(mole.species[i].location == NULL && mole_global.list[i]->charge == 0 && mole_global.list[i]->parentLabel.empty())
00156 xNeutralParticleDensity += (realnum)mole.species[i].den;
00157 }
00158
00159 {
00160 enum {DEBUG_LOC=false};
00161 if( DEBUG_LOC )
00162 {
00163 fprintf(ioQQQ," xIonDense xNeutralParticleDensity tot\t%.3e",xNeutralParticleDensity);
00164 for( long nelem=0; nelem < LIMELM; nelem++ )
00165 {
00166 fprintf(ioQQQ,"\t%.2e",dense.xIonDense[nelem][0]);
00167 }
00168 fprintf(ioQQQ,"\n");
00169 }
00170 }
00171 xValenceDonorDensity = (xNeutralParticleDensity+dense.EdenTrue);
00172 ElectronFraction = (realnum)(MAX2(MIN2(1.0,dense.EdenTrue/
00173 xValenceDonorDensity),1e-4));
00174 }
00175
00176
00177
00178
00179 realnum xBoundValenceDensity = (1.0f-ElectronFraction)*xValenceDonorDensity;
00180
00181 if( secondaries.lgSecOFF )
00182 {
00183 secondaries.HeatEfficPrimary = 1.;
00184 secondaries.SecIon2PrimaryErg = 0.;
00185 secondaries.SecExcitLya2PrimaryErg = 0.;
00186 Cosmic_ray_heat_eff = 0.95;
00187 Cosmic_ray_sec2prim = 0.05f;
00188 }
00189
00190
00191 else
00192 {
00193
00194
00195
00196
00197
00198
00199 secondaries.HeatEfficPrimary = 0.9971f*(1.f - pow(1.f - pow(ElectronFraction,(realnum)0.2663f),(realnum)1.3163f));
00200
00201
00202
00203
00204
00205 secondaries.SecIon2PrimaryErg = 0.3908f*pow(1.f - pow(ElectronFraction,(realnum)0.4092f),(realnum)1.7592f);
00206
00207
00208 secondaries.SecIon2PrimaryErg /= ((realnum)EN1RYD);
00209
00210
00211
00212
00213 sec2prim_par_1 = -(1.251f + 195.932f * ElectronFraction);
00214 sec2prim_par_2 = 1.f + 46.814f * ElectronFraction - 44.969f *
00215 ElectronFraction * ElectronFraction;
00216
00217 Cosmic_ray_sec2prim = (exp(sec2prim_par_1/SDIV( sec2prim_par_2)));
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229 secondaries.SecExcitLya2PrimaryErg = 0.4766f*pow(1.f - pow(ElectronFraction,(realnum)0.2735f),(realnum)1.5221f)/((realnum)EN1RYD);
00230
00231 if( (trace.lgTrace && trace.lgSecIon) )
00232 {
00233 fprintf( ioQQQ,
00234 " csupra H0 %.2e HeatSum eval sec ion effic, ElectronFraction = %.3e HeatEfficPrimary %.3e SecIon2PrimaryErg %.3e\n",
00235 secondaries.csupra[ipHYDROGEN][0],
00236 ElectronFraction,
00237 secondaries.HeatEfficPrimary ,
00238 secondaries.SecIon2PrimaryErg );
00239 }
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268 Cosmic_ray_heat_eff = - 8.189 - 18.394 * ElectronFraction - 6.608 * ElectronFraction * ElectronFraction * log(ElectronFraction)
00269 + 8.322 * exp( ElectronFraction ) + 4.961 * sqrt(ElectronFraction);
00270 }
00271
00272
00273
00274
00275
00276
00277
00278
00279 photo_heat_2lev_atoms = 0.;
00280 photo_heat_ISO_atoms = 0.;
00281 photo_heat_UTA = 0.;
00282
00283
00284
00285
00286 SaveOxygen1 = dense.xIonDense[ipOXYGEN][0];
00287 SaveCarbon1 = dense.xIonDense[ipCARBON][0];
00288
00289
00290 fixit();
00291
00292
00293
00294
00295
00296
00297
00298 dense.xIonDense[ipOXYGEN][0] += (realnum) findspecieslocal("CO")->den;
00299 dense.xIonDense[ipCARBON][0] += (realnum) findspecieslocal("CO")->den;
00300
00301
00302 CoolHeavy.colmet = 0.;
00303
00304 secmet = 0.;
00305 smetla = 0.;
00306 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00307 {
00308 SeconIoniz_iso[ipISO] = 0.;
00309 SeconExcit_iso[ipISO] = 0.;
00310 }
00311
00312
00313
00314 for( long nelem=0; nelem<LIMELM; ++nelem)
00315 {
00316 secmetsave[nelem] = 0.;
00317 if( dense.lgElmtOn[nelem] )
00318 {
00319
00320
00321
00322
00323
00324
00325 limit = MIN2( dense.IonHigh[nelem] , nelem+1-NISO );
00326
00327
00328 for( ion=dense.IonLow[nelem]; ion < limit; ion++ )
00329 {
00330
00331 HeatingLo = 0.;
00332 HeatingHi = 0.;
00333 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
00334 {
00335
00336 HeatingLo += ionbal.PhotoRate_Shell[nelem][ion][ns][1];
00337 HeatingHi += ionbal.PhotoRate_Shell[nelem][ion][ns][2];
00338 }
00339
00340
00341
00342 thermal.heating[nelem][ion] = dense.xIonDense[nelem][ion]*
00343 (HeatingLo + HeatingHi*secondaries.HeatEfficPrimary);
00344
00345 photo_heat_UTA += dense.xIonDense[nelem][ion]*
00346 ionbal.UTA_heat_rate[nelem][ion];
00347
00348
00349 photo_heat_2lev_atoms += thermal.heating[nelem][ion];
00350
00351
00352
00353
00354
00355 CoolHeavy.colmet += dense.xIonDense[nelem][ion]*ionbal.CollIonRate_Ground[nelem][ion][1];
00356
00357
00358 secmetsave[nelem] += HeatingHi*secondaries.SecIon2PrimaryErg* dense.xIonDense[nelem][ion];
00359
00360
00361 smetla += HeatingHi*secondaries.SecExcitLya2PrimaryErg* dense.xIonDense[nelem][ion];
00362 }
00363 secmet += secmetsave[nelem];
00364
00365
00366 limit = MAX2( limit, dense.IonLow[nelem] );
00367 for( ion=MAX2(0,limit); ion < dense.IonHigh[nelem]; ion++ )
00368 {
00369
00370 ipISO = nelem-ion;
00371
00372 HeatingLo = 0.;
00373 HeatingHi = 0.;
00374
00375 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
00376 {
00377
00378 HeatingLo += ionbal.PhotoRate_Shell[nelem][ion][ns][1];
00379 HeatingHi += ionbal.PhotoRate_Shell[nelem][ion][ns][2];
00380 }
00381
00382
00383 thermal.heating[nelem][ion] = iso_sp[ipISO][nelem].st[0].Pop()*
00384 (HeatingLo + HeatingHi*secondaries.HeatEfficPrimary);
00385
00386
00387 photo_heat_UTA += dense.xIonDense[nelem][ion]*
00388 ionbal.UTA_heat_rate[nelem][ion];
00389
00390
00391 photo_heat_ISO_atoms += thermal.heating[nelem][ion];
00392
00393 if( secondaries.SecIon2PrimaryErg > SMALLFLOAT && xBoundValenceDensity > SMALLFLOAT )
00394 {
00395
00396
00397
00398 SeconIoniz_iso[ipISO] += HeatingHi*secondaries.SecIon2PrimaryErg*
00399 iso_sp[ipISO][nelem].st[0].Pop()/
00400 xBoundValenceDensity;
00401
00402
00403 SeconExcit_iso[ipISO] += HeatingHi*secondaries.SecExcitLya2PrimaryErg*
00404 iso_sp[ipISO][nelem].st[0].Pop()/
00405 xBoundValenceDensity;
00406
00407 ASSERT( SeconIoniz_iso[ipISO]>=0. &&
00408 SeconExcit_iso[ipISO]>=0.);
00409 }
00410 }
00411
00412
00413
00414 for( ion=0; ion<dense.IonLow[nelem]; ion++ )
00415 {
00416 ASSERT( thermal.heating[nelem][ion] == 0. );
00417 }
00418 for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ion++ )
00419 {
00420 ASSERT( thermal.heating[nelem][ion] == 0. );
00421 }
00422 }
00423 }
00424 if( trace.lgTrace && trace.lgSecIon )
00425 {
00426 double savemax=0.;
00427 long int ipsavemax=-1;
00428 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00429 {
00430 if( secmetsave[nelem] > savemax )
00431 {
00432 savemax = secmetsave[nelem];
00433 ipsavemax = nelem;
00434 }
00435 }
00436 fprintf( ioQQQ,
00437 " HeatSum secmet largest contributor element %li frac of total %.3e, total %.3e\n",
00438 ipsavemax,
00439 savemax/SDIV(secmet),
00440 secmet);
00441 }
00442
00443
00444 dense.xIonDense[ipOXYGEN][0] = SaveOxygen1;
00445 dense.xIonDense[ipCARBON][0] = SaveCarbon1;
00446
00447
00448
00449
00450
00451 if( secondaries.SecIon2PrimaryErg > SMALLFLOAT && xBoundValenceDensity > SMALLFLOAT )
00452 {
00453
00454 secmet /= xBoundValenceDensity;
00455 smetla /= xBoundValenceDensity;
00456 }
00457 else
00458 {
00459
00460 secmet = 0.;
00461 smetla = 0.;
00462 }
00463
00464
00465
00466
00467
00468
00469 ionbal.CompRecoilHeatLocal = 0.;
00470 for( long nelem=0; nelem<LIMELM; ++nelem )
00471 {
00472 for( ion=0; ion<nelem+1; ++ion )
00473 {
00474 ionbal.CompRecoilHeatLocal +=
00475 ionbal.CompRecoilHeatRate[nelem][ion]*secondaries.HeatEfficPrimary*dense.xIonDense[nelem][ion];
00476 }
00477 }
00478
00479
00480 ionbal.CompRecoilHeatLocal +=
00481 2.*ionbal.CompRecoilHeatRate[ipHYDROGEN][0]*secondaries.HeatEfficPrimary*hmi.H2_total;
00482
00483
00484
00485
00486 thermal.heating[0][24] = thermal.char_tran_heat;
00487
00488
00489 thermal.heating[0][21] =
00490 ionbal.PairProducPhotoRate[2]*secondaries.HeatEfficPrimary*(dense.gas_phase[ipHYDROGEN] + 4.F*dense.gas_phase[ipHELIUM]);
00491
00492
00493
00494 thermal.heating[0][20] =
00495 ionbal.xNeutronHeatRate*secondaries.HeatEfficPrimary*dense.gas_phase[ipHYDROGEN];
00496
00497
00498 thermal.heating[0][20] += ionbal.ExtraHeatRate;
00499
00500
00501 thermal.heating[0][7] = photo_heat_UTA;
00502
00503
00504
00505 thermal.heating[0][18] = 0.;
00506 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00507 {
00508
00509 thermal.heating[0][18] += (*diatom)->Abund() *
00510 ( (*diatom)->photo_heat_soft + (*diatom)->photo_heat_hard*secondaries.HeatEfficPrimary);
00511 }
00514
00515
00516
00517 thermal.heating[0][26] = (findspecieslocal("H2+")->den+findspecieslocal("H3+")->den) *
00518 (ionbal.PhotoRate_Shell[ipHYDROGEN][0][0][1] +
00519 ionbal.PhotoRate_Shell[ipHYDROGEN][0][0][2]*secondaries.HeatEfficPrimary);
00520
00521
00522
00523
00524
00525 # if 0
00526 double hneut;
00527 if( (dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1])/
00528 dense.gas_phase[ipHYDROGEN]<0.001 )
00529 {
00530
00531 hneut = dense.xIonDense[ipHYDROGEN][0] + 2.*(findspecieslocal("H2")->den+findspecieslocal("H2*")->den);
00532 }
00533 else
00534 {
00535
00536 hneut = dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1];
00537 }
00538 # endif
00539
00540
00541 thermal.heating[1][6] =
00542 ionbal.CosRayHeatNeutralParticles*
00543 xBoundValenceDensity * Cosmic_ray_heat_eff;
00544
00545
00546 thermal.heating[1][6] += ionbal.CosRayHeatThermalElectrons * dense.eden;
00547
00548
00549 OtherHeat = 0.;
00550 for( long nelem=0; nelem<LIMELM; ++nelem)
00551 {
00552
00553
00554
00555
00556 for( i=nelem+1; i < LIMELM; i++ )
00557 {
00558 OtherHeat += thermal.heating[nelem][i];
00559 }
00560 }
00561
00562 thermal.htot = OtherHeat + photo_heat_2lev_atoms + photo_heat_ISO_atoms;
00563
00564
00565 if( called.lgTalk && !conv.lgSearch )
00566 {
00567
00568 if( thermal.htot < 0. && !thermal.lgTemperatureConstant)
00569 {
00570 fprintf( ioQQQ,
00571 " Total heating is <0; is this model collisionally ionized? zone is %li\n",
00572 nzone );
00573 }
00574 else if( thermal.htot == 0. )
00575 {
00576 fprintf( ioQQQ,
00577 " Total heating is 0; is the density small? zone is %li\n",
00578 nzone);
00579 }
00580 }
00581
00582
00583
00584 if( thermal.heatl/thermal.ctot < -1e-15 )
00585 {
00586 fprintf( ioQQQ, " HeatSum gets negative line heating,%10.2e%10.2e this is insane.\n",
00587 thermal.heatl, thermal.ctot );
00588
00589 fprintf( ioQQQ, " this is zone%4ld\n", nzone );
00590 ShowMe();
00591 cdEXIT(EXIT_FAILURE);
00592 }
00593
00594 fixit();
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607 if( secondaries.SecIon2PrimaryErg > SMALLFLOAT && xBoundValenceDensity > SMALLFLOAT )
00608 {
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625 realnum DensityRatio = (dense.gas_phase[ipHYDROGEN] +
00626 4.F*dense.gas_phase[ipHELIUM])/xBoundValenceDensity;
00627
00628 pair_production_ionization_rate =
00629 ionbal.PairProducPhotoRate[2]*secondaries.SecIon2PrimaryErg*DensityRatio;
00630
00631
00632 SecExcitLyaRate =
00633 ionbal.PairProducPhotoRate[2]*secondaries.SecExcitLya2PrimaryErg*1.3333*DensityRatio;
00634
00635
00636 fast_neutron_ionization_rate =
00637 ionbal.xNeutronHeatRate*secondaries.SecIon2PrimaryErg*DensityRatio;
00638
00639
00640 SecExcitLyaRate +=
00641 ionbal.xNeutronHeatRate*secondaries.SecExcitLya2PrimaryErg*1.3333*DensityRatio;
00642
00643
00644 SecExcitLyaRate +=
00645
00646 ionbal.CosRayHeatNeutralParticles*secondaries.SecExcitLya2PrimaryErg*1.3333*DensityRatio;
00647 }
00648 else
00649 {
00650
00651 pair_production_ionization_rate = 0.;
00652 SecExcitLyaRate = 0.;
00653 fast_neutron_ionization_rate = 0.;
00654 }
00655
00656
00657
00658 cosmic_ray_ionization_rate =
00659
00660
00661
00662
00663 ionbal.CosRayIonRate*Cosmic_ray_sec2prim +
00664
00665 ionbal.CosRayHeatNeutralParticles*secondaries.SecIon2PrimaryErg;
00666
00667
00668
00669
00670
00671
00672
00673 if( secondaries.lgCSetOn )
00674 {
00675 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00676 {
00677 for( ion=0; ion<nelem+1; ++ion )
00678 {
00679 secondaries.csupra[nelem][ion] = secondaries.SetCsupra*secondaries.csupra_effic[nelem][ion];
00680 }
00681 }
00682 secondaries.x12tot = secondaries.SetCsupra;
00683 }
00684 else
00685 {
00686 double csupra = cosmic_ray_ionization_rate +
00687 pair_production_ionization_rate +
00688 fast_neutron_ionization_rate +
00689 SeconIoniz_iso[ipH_LIKE] +
00690 SeconIoniz_iso[ipHE_LIKE] +
00691 secmet;
00692
00694
00695 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00696 {
00697 for( ion=0; ion<nelem+1; ++ion )
00698 {
00699 secondaries.csupra[nelem][ion] = (realnum)csupra*secondaries.csupra_effic[nelem][ion];
00700 }
00701 }
00702
00703 fixit();
00704
00705 secondaries.x12tot = SecExcitLyaRate +
00706 SeconExcit_iso[ipH_LIKE] +
00707 SeconExcit_iso[ipHE_LIKE] +
00708 smetla;
00709
00710 if (0)
00711 fprintf(ioQQQ,"x12 %ld %15.8g %15.8g %15.8g\n",nzone,
00712 secondaries.x12tot,
00713 secondaries.SecExcitLya2PrimaryErg,ElectronFraction);
00714 }
00715
00716 if( trace.lgTrace && secondaries.csupra[ipHYDROGEN][0] > 0. )
00717 {
00718 fprintf( ioQQQ,
00719 " HeatSum return CSUPRA %9.2e SECCMP %6.3f SecHI %6.3f SECHE %6.3f SECMET %6.3f efrac= %9.2e \n",
00720 secondaries.csupra[ipHYDROGEN][0],
00721 (cosmic_ray_ionization_rate+pair_production_ionization_rate+fast_neutron_ionization_rate)/secondaries.csupra[ipHYDROGEN][0],
00722 SeconIoniz_iso[ipH_LIKE] / secondaries.csupra[ipHYDROGEN][0],
00723 SeconIoniz_iso[ipHE_LIKE]/secondaries.csupra[ipHYDROGEN][0],
00724 secmet/secondaries.csupra[ipHYDROGEN][0] ,
00725 ElectronFraction );
00726 }
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737 thermal.dHeatdT = -0.7*(photo_heat_2lev_atoms+photo_heat_ISO_atoms+
00738 photo_heat_UTA)/phycon.te;
00739
00740
00741
00742 thermal.dHeatdT *= dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN];
00743 if( PRT_DERIV )
00744 fprintf(ioQQQ,"DEBUG dhdt 0 %.3e %.3e %.3e \n",
00745 thermal.dHeatdT,
00746 photo_heat_2lev_atoms,
00747 photo_heat_ISO_atoms);
00748
00749
00750
00751 oldfac = -0.7/phycon.te;
00752
00753
00754 thermal.dHeatdT += thermal.heating[0][9]*oldfac;
00755
00756
00757 thermal.dHeatdT += thermal.heating[0][10]*oldfac;
00758
00759
00760 thermal.dHeatdT += thermal.heating[0][22]*oldfac;
00761 if( PRT_DERIV )
00762 fprintf(ioQQQ,"DEBUG dhdt 1 %.3e \n", thermal.dHeatdT);
00763
00764
00765
00766
00767
00768 thermal.dHeatdT += CoolHeavy.brems_heat_total*(-1.5/phycon.te);
00769
00770
00771
00772 thermal.dHeatdT += gv.dHeatdT;
00773
00774
00775 thermal.dHeatdT += thermal.heating[1][2]*oldfac;
00776 if( PRT_DERIV )
00777 fprintf(ioQQQ,"DEBUG dhdt 2 %.3e \n", thermal.dHeatdT);
00778
00779
00780
00781
00782 if( hmi.HeatH2Dexc_used > 0. )
00783 thermal.dHeatdT += hmi.deriv_HeatH2Dexc_used;
00784
00785
00786
00787 thermal.dHeatdT += thermal.heating[0][15]*0.7/phycon.te;
00788
00789
00790 thermal.dHeatdT += thermal.heating[0][16]*oldfac;
00791
00792
00793
00794
00795 thermal.dHeatdT += thermal.heating[0][1]*oldfac;
00796 if( PRT_DERIV )
00797 fprintf(ioQQQ,"DEBUG dhdt 3 %.3e \n", thermal.dHeatdT);
00798
00799
00800 thermal.dHeatdT += thermal.heating[0][3]*oldfac;
00801
00802
00803
00804
00805
00806 thermal.dHeatdT += thermal.heating[0][19]*oldfac;
00807
00808
00809 thermal.dHeatdT += thermal.heating[0][20]*oldfac;
00810
00811
00812 thermal.dHeatdT += thermal.heating[0][21]*oldfac;
00813 if( PRT_DERIV )
00814 fprintf(ioQQQ,"DEBUG dhdt 4 %.3e \n", thermal.dHeatdT);
00815
00816
00817 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00818 {
00819 for( long nelem=ipISO; nelem<LIMELM; ++nelem)
00820 {
00821 thermal.dHeatdT += iso_sp[ipISO][nelem].dLTot;
00822 }
00823 }
00824
00825
00826
00827 if( FeII.Fe2_large_heat > 0. )
00828 {
00829 thermal.dHeatdT += FeII.ddT_Fe2_large_cool;
00830 }
00831 if( PRT_DERIV )
00832 fprintf(ioQQQ,"DEBUG dhdt 5 %.3e \n", thermal.dHeatdT);
00833
00834
00835
00836 if( NumDeriv.lgNumDeriv )
00837 {
00838 if( ((nzone > 2 && nzone == nzSave) &&
00839 ! fp_equal( oldtemp, phycon.te )) && nhit > 4 )
00840 {
00841
00842
00843 deriv = (oldheat - thermal.htot)/(oldtemp - phycon.te);
00844 thermal.dHeatdT = deriv;
00845 }
00846 else
00847 {
00848 deriv = thermal.dHeatdT;
00849 }
00850
00851
00852 if( nzone != nzSave )
00853 {
00854 nhit = 0;
00855 }
00856 nzSave = nzone;
00857 nhit += 1;
00858 oldheat = thermal.htot;
00859 oldtemp = phycon.te;
00860 }
00861
00862 if( trace.lgTrace && trace.lgHeatBug )
00863 {
00864 ipnt = 0;
00865
00866 for( i=0; i < LIMELM; i++ )
00867 {
00868 for( j=0; j < LIMELM; j++ )
00869 {
00870 if( thermal.heating[i][j]/thermal.htot > FAINT_HEAT )
00871 {
00872 ipsave[ipnt] = i;
00873 jpsave[ipnt] = j;
00874 save[ipnt] = thermal.heating[i][j]/thermal.htot;
00875
00876 ipnt = MIN2((long)NDIM-1,ipnt+1);
00877 }
00878 }
00879 }
00880
00881
00882
00883
00884 for( i=0; i < thermal.ncltot; i++ )
00885 {
00886 if( thermal.heatnt[i]/thermal.htot > FAINT_HEAT )
00887 {
00888 ipsave[ipnt] = -1;
00889 jpsave[ipnt] = i;
00890 save[ipnt] = thermal.heatnt[i]/thermal.htot;
00891 ipnt = MIN2((long)NDIM-1,ipnt+1);
00892 }
00893 }
00894
00895 fprintf( ioQQQ,
00896 " HeatSum HTOT %.3e Te:%.3e dH/dT%.2e other %.2e photo 2lev %.2e photo iso %.2e\n",
00897 thermal.htot,
00898 phycon.te,
00899 thermal.dHeatdT ,
00900
00901
00902 OtherHeat ,
00903 photo_heat_2lev_atoms,
00904 photo_heat_ISO_atoms);
00905
00906 fprintf( ioQQQ, " " );
00907 for( i=0; i < ipnt; i++ )
00908 {
00909
00910 fprintf( ioQQQ, " [%ld][%ld]%6.3f",
00911 ipsave[i],
00912 jpsave[i],
00913 save[i] );
00914
00915
00916 if( !(i%8) && i>0 && i!=(ipnt-1) )
00917 {
00918 fprintf( ioQQQ, "\n " );
00919 }
00920 }
00921 fprintf( ioQQQ, "\n" );
00922 }
00923 return;
00924 }
00925
00926
00927
00928 void HeatZero( void )
00929 {
00930 long int i , j;
00931
00932 DEBUG_ENTRY( "HeatZero()" );
00933
00934 for( i=0; i < LIMELM; i++ )
00935 {
00936 for( j=0; j < LIMELM; j++ )
00937 {
00938 thermal.heating[i][j] = 0.;
00939 }
00940 }
00941 return;
00942 }