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