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