39 double cosmic_ray_ionization_rate ,
40 pair_production_ionization_rate ,
41 fast_neutron_ionization_rate , SecExcitLyaRate;
44 double SeconIoniz_iso[
NISO] ,
61 static long int nhit=0,
63 double photo_heat_2lev_atoms,
64 photo_heat_ISO_atoms ,
70 static double oldheat=0.,
74 realnum SaveOxygen1 , SaveCarbon1;
78 double Cosmic_ray_heat_eff ,
108 realnum xValenceDonorDensity, ElectronFraction;
117 long ipFullShell = -1, ipNextFull = 0;
118 xValenceDonorDensity = 0.;
122 xValenceDonorDensity +=
123 (nelem-ipFullShell)*
dense.gas_phase[nelem];
124 if (nelem == ipNoble[ipNextFull])
134 xValenceDonorDensity),1e-4));
143 realnum xNeutralParticleDensity = 0.;
144 for(
long nelem=0; nelem <
LIMELM; nelem++ )
146 xNeutralParticleDensity +=
dense.xIonDense[nelem][0];
149 xNeutralParticleDensity +=
deut.xIonDense[0];
156 xNeutralParticleDensity += (
realnum)
mole.species[i].den;
160 enum {DEBUG_LOC=
false};
163 fprintf(
ioQQQ,
" xIonDense xNeutralParticleDensity tot\t%.3e",xNeutralParticleDensity);
164 for(
long nelem=0; nelem <
LIMELM; nelem++ )
166 fprintf(
ioQQQ,
"\t%.2e",
dense.xIonDense[nelem][0]);
171 xValenceDonorDensity = (xNeutralParticleDensity+
dense.EdenTrue);
173 xValenceDonorDensity),1e-4));
179 realnum xBoundValenceDensity = (1.0f-ElectronFraction)*xValenceDonorDensity;
186 Cosmic_ray_heat_eff = 0.95;
187 Cosmic_ray_sec2prim = 0.05f;
213 sec2prim_par_1 = -(1.251f + 195.932f * ElectronFraction);
214 sec2prim_par_2 = 1.f + 46.814f * ElectronFraction - 44.969f *
215 ElectronFraction * ElectronFraction;
217 Cosmic_ray_sec2prim = (exp(sec2prim_par_1/
SDIV( sec2prim_par_2)));
234 " csupra H0 %.2e HeatSum eval sec ion effic, ElectronFraction = %.3e HeatEfficPrimary %.3e SecIon2PrimaryErg %.3e\n",
268 Cosmic_ray_heat_eff = - 8.189 - 18.394 * ElectronFraction - 6.608 * ElectronFraction * ElectronFraction * log(ElectronFraction)
269 + 8.322 * exp( ElectronFraction ) + 4.961 * sqrt(ElectronFraction);
279 photo_heat_2lev_atoms = 0.;
280 photo_heat_ISO_atoms = 0.;
308 SeconIoniz_iso[ipISO] = 0.;
309 SeconExcit_iso[ipISO] = 0.;
314 for(
long nelem=0; nelem<
LIMELM; ++nelem)
316 secmetsave[nelem] = 0.;
317 if(
dense.lgElmtOn[nelem] )
328 for( ion=
dense.IonLow[nelem]; ion < limit; ion++ )
333 for( ns=0; ns <
Heavy.nsShells[nelem][ion]; ns++ )
336 HeatingLo +=
ionbal.PhotoRate_Shell[nelem][ion][ns][1];
337 HeatingHi +=
ionbal.PhotoRate_Shell[nelem][ion][ns][2];
342 thermal.heating[nelem][ion] =
dense.xIonDense[nelem][ion]*
343 (HeatingLo + HeatingHi*
secondaries.HeatEfficPrimary);
345 photo_heat_UTA +=
dense.xIonDense[nelem][ion]*
346 ionbal.UTA_heat_rate[nelem][ion];
349 photo_heat_2lev_atoms +=
thermal.heating[nelem][ion];
358 secmetsave[nelem] += HeatingHi*
secondaries.SecIon2PrimaryErg*
dense.xIonDense[nelem][ion];
361 smetla += HeatingHi*
secondaries.SecExcitLya2PrimaryErg*
dense.xIonDense[nelem][ion];
363 secmet += secmetsave[nelem];
366 limit =
MAX2( limit,
dense.IonLow[nelem] );
367 for( ion=
MAX2(0,limit); ion <
dense.IonHigh[nelem]; ion++ )
375 for( ns=0; ns <
Heavy.nsShells[nelem][ion]; ns++ )
378 HeatingLo +=
ionbal.PhotoRate_Shell[nelem][ion][ns][1];
379 HeatingHi +=
ionbal.PhotoRate_Shell[nelem][ion][ns][2];
383 thermal.heating[nelem][ion] =
iso_sp[ipISO][nelem].st[0].Pop()*
384 (HeatingLo + HeatingHi*
secondaries.HeatEfficPrimary);
387 photo_heat_UTA +=
dense.xIonDense[nelem][ion]*
388 ionbal.UTA_heat_rate[nelem][ion];
391 photo_heat_ISO_atoms +=
thermal.heating[nelem][ion];
398 SeconIoniz_iso[ipISO] += HeatingHi*
secondaries.SecIon2PrimaryErg*
399 iso_sp[ipISO][nelem].st[0].Pop()/
400 xBoundValenceDensity;
403 SeconExcit_iso[ipISO] += HeatingHi*
secondaries.SecExcitLya2PrimaryErg*
404 iso_sp[ipISO][nelem].st[0].Pop()/
405 xBoundValenceDensity;
407 ASSERT( SeconIoniz_iso[ipISO]>=0. &&
408 SeconExcit_iso[ipISO]>=0.);
414 for( ion=0; ion<
dense.IonLow[nelem]; ion++ )
418 for( ion=
dense.IonHigh[nelem]+1; ion<nelem+1; ion++ )
427 long int ipsavemax=-1;
430 if( secmetsave[nelem] > savemax )
432 savemax = secmetsave[nelem];
437 " HeatSum secmet largest contributor element %li frac of total %.3e, total %.3e\n",
439 savemax/
SDIV(secmet),
454 secmet /= xBoundValenceDensity;
455 smetla /= xBoundValenceDensity;
469 ionbal.CompRecoilHeatLocal = 0.;
470 for(
long nelem=0; nelem<
LIMELM; ++nelem )
472 for( ion=0; ion<nelem+1; ++ion )
474 ionbal.CompRecoilHeatLocal +=
480 ionbal.CompRecoilHeatLocal +=
501 thermal.heating[0][7] = photo_heat_UTA;
509 thermal.heating[0][18] += (*diatom)->Abund() *
510 ( (*diatom)->photo_heat_soft + (*diatom)->photo_heat_hard*
secondaries.HeatEfficPrimary);
542 ionbal.CosRayHeatNeutralParticles*
543 xBoundValenceDensity * Cosmic_ray_heat_eff;
550 for(
long nelem=0; nelem<
LIMELM; ++nelem)
556 for( i=nelem+1; i <
LIMELM; i++ )
558 OtherHeat +=
thermal.heating[nelem][i];
562 thermal.htot = OtherHeat + photo_heat_2lev_atoms + photo_heat_ISO_atoms;
571 " Total heating is <0; is this model collisionally ionized? zone is %li\n",
577 " Total heating is 0; is the density small? zone is %li\n",
586 fprintf(
ioQQQ,
" HeatSum gets negative line heating,%10.2e%10.2e this is insane.\n",
589 fprintf(
ioQQQ,
" this is zone%4ld\n",
nzone );
628 pair_production_ionization_rate =
633 ionbal.PairProducPhotoRate[2]*
secondaries.SecExcitLya2PrimaryErg*1.3333*DensityRatio;
636 fast_neutron_ionization_rate =
646 ionbal.CosRayHeatNeutralParticles*
secondaries.SecExcitLya2PrimaryErg*1.3333*DensityRatio;
651 pair_production_ionization_rate = 0.;
652 SecExcitLyaRate = 0.;
653 fast_neutron_ionization_rate = 0.;
658 cosmic_ray_ionization_rate =
663 ionbal.CosRayIonRate*Cosmic_ray_sec2prim +
677 for( ion=0; ion<nelem+1; ++ion )
686 double csupra = cosmic_ray_ionization_rate +
687 pair_production_ionization_rate +
688 fast_neutron_ionization_rate +
697 for( ion=0; ion<nelem+1; ++ion )
711 fprintf(
ioQQQ,
"x12 %ld %15.8g %15.8g %15.8g\n",
nzone,
713 secondaries.SecExcitLya2PrimaryErg,ElectronFraction);
719 " HeatSum return CSUPRA %9.2e SECCMP %6.3f SecHI %6.3f SECHE %6.3f SECMET %6.3f efrac= %9.2e \n",
721 (cosmic_ray_ionization_rate+pair_production_ionization_rate+fast_neutron_ionization_rate)/
secondaries.csupra[
ipHYDROGEN][0],
737 thermal.dHeatdT = -0.7*(photo_heat_2lev_atoms+photo_heat_ISO_atoms+
738 photo_heat_UTA)/
phycon.te;
744 fprintf(
ioQQQ,
"DEBUG dhdt 0 %.3e %.3e %.3e \n",
746 photo_heat_2lev_atoms,
747 photo_heat_ISO_atoms);
782 if(
hmi.HeatH2Dexc_used > 0. )
819 for(
long nelem=ipISO; nelem<
LIMELM; ++nelem)
827 if(
FeII.Fe2_large_heat > 0. )
852 if(
nzone != nzSave )
866 for( i=0; i <
LIMELM; i++ )
868 for( j=0; j <
LIMELM; j++ )
876 ipnt =
MIN2((
long)NDIM-1,ipnt+1);
884 for( i=0; i <
thermal.ncltot; i++ )
891 ipnt =
MIN2((
long)NDIM-1,ipnt+1);
896 " HeatSum HTOT %.3e Te:%.3e dH/dT%.2e other %.2e photo 2lev %.2e photo iso %.2e\n",
903 photo_heat_2lev_atoms,
904 photo_heat_ISO_atoms);
906 fprintf(
ioQQQ,
" " );
907 for( i=0; i < ipnt; i++ )
910 fprintf(
ioQQQ,
" [%ld][%ld]%6.3f",
916 if( !(i%8) && i>0 && i!=(ipnt-1) )
918 fprintf(
ioQQQ,
"\n " );
921 fprintf(
ioQQQ,
"\n" );