33 double GrainNumRelHydrSilicate ,
34 GrainNumRelHydrCarbonaceous ,
36 GrainMassRelHydrSilicate,
37 GrainMassRelHydrCarbonaceous,
46 for( i=0; i <
LIMELM; i++ )
48 if(
dense.lgElmtOn[i] )
59 fprintf(
ioQQQ,
" \n" );
67 for( i=0; i <
LIMELM; i++ )
79 fprintf(
ioQQQ,
" \n" );
82 GrainNumRelHydrSilicate = 0.;
83 GrainNumRelHydrCarbonaceous = 0;
84 GrainNumRelHydr_PAH = 0.;
85 GrainMassRelHydrSilicate = 0.;
86 GrainMassRelHydrCarbonaceous = 0;
87 GrainMassRelHydr_PAH = 0.;
89 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
95 realnum DensityNumberPerHydrogen =
96 (
gv.bin[nd]->IntVol/
gv.bin[nd]->AvVol)*
gv.bin[nd]->dstAbund /
99 realnum DensityMassPerHydrogen =
100 gv.bin[nd]->IntVol*
gv.bin[nd]->dustp[0]*
gv.bin[nd]->dstAbund/
107 GrainNumRelHydrCarbonaceous += DensityNumberPerHydrogen;
108 GrainMassRelHydrCarbonaceous += DensityMassPerHydrogen;
113 GrainNumRelHydrSilicate += DensityNumberPerHydrogen;
114 GrainMassRelHydrSilicate += DensityMassPerHydrogen;
121 GrainNumRelHydr_PAH += DensityNumberPerHydrogen;
122 GrainMassRelHydr_PAH += DensityMassPerHydrogen;
129 fprintf(
ioQQQ,
" Number of grains per hydrogen (scale=1) Mass of grains per hydrogen (scale=1)\n");
130 fprintf(
ioQQQ,
" Carbonaceous: %.3f Silicate: %.3f PAH: %.3f Carbonaceous: %.3f Silicate: %.3f PAH: %.3f\n\n" ,
131 log10(
MAX2( 1e-30, GrainNumRelHydrCarbonaceous ) ) ,
132 log10(
MAX2( 1e-30, GrainNumRelHydrSilicate ) ) ,
133 log10(
MAX2( 1e-30, GrainNumRelHydr_PAH ) ) ,
134 log10(
MAX2( 1e-30, GrainMassRelHydrCarbonaceous ) ) ,
135 log10(
MAX2( 1e-30, GrainMassRelHydrSilicate ) ) ,
136 log10(
MAX2( 1e-30, GrainMassRelHydr_PAH ) ) );
148 static bool lgFirstCall=
true;
149 static bool lgElOnOff[
LIMELM];
166 lgElOnOff[i] =
dense.lgElmtOn[i];
175 dense.lgElmtOn[i] = lgElOnOff[i] &&
dense.lgElmtOn[i];
192 if(
dense.lgDenFlucOn )
206 abund.ScaleElement[i]*fac);
214 if(
abund.lgAbunTabl[nelem] )
231 if(
abund.solar[nelem] <
dense.AbundanceLimit )
232 dense.lgElmtOn[nelem] =
false;
234 if(
dense.lgElmtOn[nelem] )
237 if(
dense.gas_phase[nelem] <= 0. )
239 fprintf(
ioQQQ,
" Abundances must be greater than zero. "
240 "Check entered abundance for element%3ld = %2.2s\n",
246 fprintf(
ioQQQ,
" Abundance for %s is %.2e, less than lower "
247 "limit of %.3e, so turning element off.\n",
249 dense.gas_phase[nelem],
251 dense.lgElmtOn[nelem] =
false;
255 fprintf(
ioQQQ,
" Abundance for %s is %.2e. This version of Cloudy does not "
256 "permit densities greater than %e cm-3.\n",
258 dense.gas_phase[nelem],
263 if( !
dense.lgElmtOn[nelem] )
266 dense.SetGasPhaseDensity( nelem, 0. );
270 dense.xIonDense[nelem][0] =
dense.gas_phase[nelem];
271 for(
long int ion=1; ion <
LIMELM+1; ion++ )
273 dense.xIonDense[nelem][ion] = 0.;
290 fprintf(
ioQQQ,
"\n >>> \n"
291 " >>> The simulation is going into possibly molecular gas but the carbon/oxygen abundance ratio is greater than unity.\n" );
292 fprintf(
ioQQQ,
" >>> Standard interstellar chemistry networks are designed for environments with C/O < 1.\n" );
293 fprintf(
ioQQQ,
" >>> The chemistry network may (or may not) collapse deep in molecular regions where CO is fully formed.\n" );
294 fprintf(
ioQQQ,
" >>> \n\n\n\n\n" );
300 realnum sumx , sumy , sumz = 0.;
305 fprintf(
ioQQQ,
"\n AbundancesSet sets following densities (cm^-3); \n" );
308 for( nelem=i*10; nelem < i*10+10; nelem++ )
313 sumz +=
dense.gas_phase[nelem]*
dense.AtomicWeight[nelem];
315 fprintf(
ioQQQ,
" \n" );
317 fprintf(
ioQQQ,
"\n AbundancesSet sets following abundances rel to H; \n" );
320 for( nelem=i*10; nelem < i*10+10; nelem++ )
325 fprintf(
ioQQQ,
" \n" );
327 fprintf(
ioQQQ,
" \n" );
328 fprintf(
ioQQQ,
" Gas-phase mass fractions, X:%.3e Y:%.3e Z:%.3e\n\n",
329 sumx/
SDIV(sumx+sumy+sumz) ,
330 sumy/
SDIV(sumx+sumy+sumz) ,
331 sumz/
SDIV(sumx+sumy+sumz) );
351 static long int nelem;
355 if( strcmp(chJob,
"initG") == 0 )
360 " Gas Phase Chemical Composition\n" );
362 else if( strcmp(chJob,
"initD") == 0 )
367 " Grain Chemical Composition\n" );
370 else if( strcmp(chJob,
"fill") == 0 )
373 abund_prt = log10( abund_prt );
375 sprintf( chAllLabels[nelem],
" %2.2s:%8.4f", chLabl, abund_prt );
379 fprintf(
ioQQQ,
" " );
382 fprintf(
ioQQQ,
"%13.13s", chAllLabels[i] );
384 fprintf(
ioQQQ,
"\n" );
397 else if( strcmp(chJob,
"fillp") == 0 )
400 abund_prt = log10( abund_prt );
403 sprintf( chAllLabels[nelem],
" %2.2s:%8.4f", chLabl, abund_prt );
407 fprintf(
ioQQQ,
" " );
410 fprintf(
ioQQQ,
"%13.13s", chAllLabels[i] );
412 fprintf(
ioQQQ,
"\n" );
424 else if( strcmp(chJob,
"flus") == 0 )
430 fprintf(
ioQQQ,
" " );
432 for(i=0; i < noffset; i++)
435 fprintf(
ioQQQ,
" " );
439 if( !(nelem%2) && nelem > 0)
442 for( i=0; i < nelem; i++ )
444 fprintf(
ioQQQ,
"%13.13s", chAllLabels[i] );
447 fprintf(
ioQQQ,
"\n" );
451 fprintf(
ioQQQ,
" PrtElem does not understand job=%4.4s\n",