35 double cosmic_ray_ionization_rate ,
36 pair_production_ionization_rate ,
37 fast_neutron_ionization_rate , SecExcitLyaRate;
40 double SeconIoniz_iso[
NISO] ,
57 static long int nhit=0,
59 double photo_heat_2lev_atoms,
60 photo_heat_ISO_atoms ,
66 static double oldheat=0.,
70 realnum SaveOxygen1 , SaveCarbon1;
73 realnum xNeutralParticleDensity;
78 double Cosmic_ray_heat_eff ,
94 xNeutralParticleDensity = 0.;
95 for( nelem=0; nelem <
LIMELM; nelem++ )
104 if(
COmole[i]->n_nuclei > 1)
112 enum {DEBUG_LOC=
false};
115 fprintf(
ioQQQ,
" xIonDense xNeutralParticleDensity tot\t%.3e",xNeutralParticleDensity);
116 for( nelem=0; nelem <
LIMELM; nelem++ )
143 static realnum OldElectronFraction = 0,
144 OlderElectronFraction = 0;
147 OldElectronFraction = 0;
148 OlderElectronFraction = 0;
150 if( (ElectronFraction-OldElectronFraction)*
151 (OldElectronFraction-OlderElectronFraction) < 0. )
152 ElectronFraction = ( ElectronFraction+OldElectronFraction)/2.f;
154 OlderElectronFraction = OldElectronFraction;
155 OldElectronFraction = ElectronFraction;
162 Cosmic_ray_heat_eff = 0.95;
163 Cosmic_ray_sec2prim = 0.05f;
189 sec2prim_par_1 = -(1.251f + 195.932f * ElectronFraction);
190 sec2prim_par_2 = 1.f + 46.814f * ElectronFraction - 44.969f *
191 ElectronFraction * ElectronFraction;
193 Cosmic_ray_sec2prim = (exp(sec2prim_par_1/
SDIV( sec2prim_par_2)));
210 " csupra H0 %.2e HeatSum eval sec ion effic, ElectronFraction = %.3e HeatEfficPrimary %.3e SecIon2PrimaryErg %.3e\n",
244 Cosmic_ray_heat_eff = - 8.189 - 18.394 * ElectronFraction - 6.608 * ElectronFraction * ElectronFraction * log(ElectronFraction)
245 + 8.322 * exp( ElectronFraction ) + 4.961 * sqrt(ElectronFraction);
255 photo_heat_2lev_atoms = 0.;
256 photo_heat_ISO_atoms = 0.;
276 SeconIoniz_iso[ipISO] = 0.;
277 SeconExcit_iso[ipISO] = 0.;
282 for( nelem=0; nelem<
LIMELM; ++nelem)
284 secmetsave[nelem] = 0.;
331 secmet += secmetsave[nelem];
368 xNeutralParticleDensity;
373 xNeutralParticleDensity;
375 ASSERT( SeconIoniz_iso[ipISO]>=0. &&
376 SeconExcit_iso[ipISO]>=0.);
395 long int ipsavemax=-1;
398 if( secmetsave[nelem] > savemax )
400 savemax = secmetsave[nelem];
405 " HeatSum secmet largest contributor element %li frac of total %.3e, total %.3e\n",
407 savemax/
SDIV(secmet),
422 secmet /= xNeutralParticleDensity;
423 smetla /= xNeutralParticleDensity;
438 for( nelem=0; nelem<
LIMELM; ++nelem )
440 for( ion=0; ion<nelem+1; ++ion )
507 xNeutralParticleDensity * Cosmic_ray_heat_eff;
514 for( nelem=0; nelem<
LIMELM; ++nelem)
520 for( i=nelem+1; i <
LIMELM; i++ )
526 thermal.
htot = OtherHeat + photo_heat_2lev_atoms + photo_heat_ISO_atoms;
535 " Total heating is <0; is this model collisionally ionized? zone is %li\n",
541 " Total heating is 0; is the density small? zone is %li\n",
550 fprintf(
ioQQQ,
" HeatSum gets negative line heating,%10.2e%10.2e this is insane.\n",
553 fprintf(
ioQQQ,
" this is zone%4ld\n",
nzone );
574 pair_production_ionization_rate =
582 fast_neutron_ionization_rate =
597 pair_production_ionization_rate = 0.;
598 SecExcitLyaRate = 0.;
599 fast_neutron_ionization_rate = 0.;
604 cosmic_ray_ionization_rate =
623 for( ion=0; ion<nelem+1; ++ion )
633 double facold , facnew;
645 facnew = 1. - facold;
647 (cosmic_ray_ionization_rate +
648 pair_production_ionization_rate +
649 fast_neutron_ionization_rate +
658 for( ion=0; ion<nelem+1; ++ion )
676 " HeatSum return CSUPRA %9.2e SECCMP %6.3f SecHI %6.3f SECHE %6.3f SECMET %6.3f efrac= %9.2e \n",
694 thermal.
dHeatdT = -0.7*(photo_heat_2lev_atoms+photo_heat_ISO_atoms+
701 fprintf(
ioQQQ,
"DEBUG dhdt 0 %.3e %.3e %.3e \n",
703 photo_heat_2lev_atoms,
704 photo_heat_ISO_atoms);
776 for( nelem=ipISO; nelem<
LIMELM; ++nelem)
809 if(
nzone != nzSave )
823 for( i=0; i <
LIMELM; i++ )
825 for( j=0; j <
LIMELM; j++ )
833 ipnt =
MIN2((
long)NDIM-1,ipnt+1);
848 ipnt =
MIN2((
long)NDIM-1,ipnt+1);
853 " HeatSum HTOT %.3e Te:%.3e dH/dT%.2e other %.2e photo 2lev %.2e photo iso %.2e\n",
860 photo_heat_2lev_atoms,
861 photo_heat_ISO_atoms);
863 fprintf(
ioQQQ,
" " );
864 for( i=0; i < ipnt; i++ )
867 fprintf(
ioQQQ,
" [%ld][%ld]%6.3f",
873 if( !(i%8) && i>0 && i!=(ipnt-1) )
875 fprintf(
ioQQQ,
"\n " );
878 fprintf(
ioQQQ,
"\n" );
891 for( i=0; i <
LIMELM; i++ )
893 for( j=0; j <
LIMELM; j++ )