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