Hot fixes to C06.02


This page lists all hot fixes to C06.02. In each case the hot fixes have been incorporated in the next patch release. That is, all hot fixes to C06.02c are included in C06.02d.


fixes to version C06.02c

An assert can be thrown at line 176 or 177 of cont_ffun.c. Remove these two asserts, which read as follows:

	ASSERT( y >= MIN2( rfield.tFluxLog[rfield.ipspec][i-1] , rfield.tFluxLog[rfield.ipspec][i] ));
	ASSERT( y <= MAX2( rfield.tFluxLog[rfield.ipspec][i-1] , rfield.tFluxLog[rfield.ipspec][i] ));

and replace with the following block

# 	ifndef NDEBUG
	{
	double ys1 = rfield.tFluxLog[rfield.ipspec][i-1];
	double ys2 = rfield.tFluxLog[rfield.ipspec][i];
	double dy = ys2-ys1;
	/* apply safety margin around end points */
	double margin = 3.*sqrt(POW2(ys1) + POW2(ys2-ys1))*FLT_EPSILON;
	ys1 -= sign(margin,dy);
	ys2 += sign(margin,dy);
	ASSERT( y >= MIN2( ys1, ys2 ) );
	ASSERT( y <= MAX2( ys1, ys2 ) );
	}
# 	endif

Thanks to astronerma@… for reporting the crash, and to Peter van Hoof for fixing it. 2007 Feb 01.


The power law command could result in a floating point exception when the vary option was included. This was caused by incrementing the number of continua before the command was fully parsed. To correct, edit the file parse_powerlawcontinuum.c and move the block of code between lines 92 and 99, which reads as follows:

	++rfield.nspec;
	if( rfield.nspec >= LIMSPC )
	{
		fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
		puts( "[Stop in ParsePowerlawContinuum]" );
		cdEXIT(EXIT_FAILURE);
	}

to after the block that ends at line 153. After the change the code will read as follows (this includes some code before and after the change):

		}
		++optimize.nparm;
	}

	/*>>chng 06 nov 10, BUGFIX, nspec was incremented before previous branch
	 * and so crashed with log10 0 since nspec was beyond set values
	 * caused fpe domain function error with log 0
	 * fpe caught by Pavel Abolmasov */
	++rfield.nspec;
	if( rfield.nspec >= LIMSPC )
	{
		fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
		puts( "[Stop in ParsePowerlawContinuum]" );
		cdEXIT(EXIT_FAILURE);
	}

#	ifdef DEBUG_FUN
	fputs( " <->ParsePowerlawContinuum()\n", debug_fp );
#	endif

Many thanks to Pavel Abolmasov for reporting this crash. This change is only needed if you want to vary the options on the power law command. 06 Nov 10.


fixes to version 06.02b

The table stars command could cause crash on certain types of Macs. Edit file stars_rebinatmosphere.c and change the block starting at line 89 so that it reads as follows:

for( j=0; j < nCont-1; j++ )
{
	/* >>chng 06 aug 11, on some systems (e.g., macbook pro) y/x can get evaluated as y*(1/x);
	 * this causes overflows if x is a denormalized number, hence we force a cast to double, PvH */
	double ratio_x = (double)StarEner[j+1]/(double)StarEner[j];
	double ratio_y = (double)StarFlux[j+1]/(double)StarFlux[j];
	StarPower[j] = (float)(log(ratio_y)/log(ratio_x));
}

This need only be done if the table stars commands do not work on your Mac. If does not affect the calculations other than preventing the crash on some Macs.

Many thanks to Dick Henry for uncovering the problem and to Peter van Hoof for fixing it. 2006 Aug 12.


The code could incorrectly detect level populations insanity and end with an assert. Some important quantities are calculated with two independent methods, their results compared, and insanity declared if the results differ. This is a basic concept in safe programming (discussed in this link). The predicted densities of trace stages of ionization could differ by a fraction of a percent when the ionization fraction fell below one part in 1015 due to numerical imprecision. The code detected this and threw an assert around line 138 of radius_increment.c. The fix is to not do this test for ions which have extremely small fractional abundances.

Edit the file radius_increment.c and change line 99 and add new lines so that they read as follows

	if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][nelem-ipISO]/dense.eden>conv.EdenErrorAllowed / 10. &&
	/*>>chng 06 jul 23, add test on inaccurate results with LotT solver */
	!(iso.chTypeAtomUsed[ipISO][nelem][0]=='L' && iso.xIonSimple[ipISO][nelem]<1e-10) )

Change line 123 and add two new lines of code:

	if( iso.Pop2Ion[ipISO][nelem][0] >
		dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk &&
		/*>>chng 06 jul 22, only do check if pops are within tolerance of double */
		iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 )

Change line 142 and add new lines so that they read as

	ASSERT( conv.lgSearch ||
	iso.Pop2Ion[ipISO][nelem][0] <=
	dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk ||
	/*>>chng 06 jul 22, only do check if pops are within tolerance of double */
	iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15);

Change line 148 and add two lines so that it reads as

	if( EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo*dense.xIonDense[nelem][nelem+1-ipISO] >
	dense.xIonDense[nelem][nelem-ipISO]*err_chk &&
	/*>>chng 06 jul 22, only do check if pops are within tolerance of double */
	iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15)

Finally change line 160 so that it reads as

	ASSERT( (conv.lgSearch || !conv.nPres2Ioniz) || 
	/*(conv.lgSearch || !conv.nPres2Ioniz) || */
	(EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo*dense.xIonDense[nelem][nelem+1-ipISO] <=
		dense.xIonDense[nelem][nelem-ipISO]*err_chk) ||
	/*>>chng 06 jul 22, only do check if pops are within tolerance of double */
	iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15);

This fix need only be made if the code throws an assert near line 138 of radius_increment.c. The correction will not affect results. Thanks to Joseph Meiring for reporting this problem. 2006 July 23


The mean H0 gas temperature that is weighted by the 21 cm opacity should be weighted over radius but was weighted with respect to volume. This is given with the label T(<nH/Tspin> in the final printout. To fix this edit file prt_final.c and change line 1515 from

	/* get harmonic mean HI gas kin temperature, as needed for 21 cm spin temperature */
	if( cdTemp( "21cm",0,&THI,"volume" ) )

to read as follows:

	/* get harmonic mean HI gas kin temperature, as needed for 21 cm spin temperature 
	 * >>chng 06 jul 21, this was over volume but hazy said radius - change to radius,
	 * bug discovered by Eric Pellegrini & Jack Baldwin  */
	if( cdTemp( "21cm",0,&THI,"radius" ) )

Thanks to Eric Pellegrini & Jack Baldwin for catching this problem. 2006 July 21.


fixes to version 06.02a

The UV FeII optical depths predicted by the large FeII model are not correct when gas becomes very cold. The large model trims down the number of levels it computes when gas become very cold and the populations of highly excited levels becomes very small. When this happened the code did not continue to update optical depths of lines between lower levels and these high levels that are not computed. As a result the optical depth printed for the transition is too small. The fix, which will be part of 06.02b, is to continue to update optical depths for such lines. 2006 Jun 25


Thermal instabilities, followed by a thrown assert, could occur in a cosmic-ray dominated gas. There are two fixes needed. Edit file helike_recom.c and add the following after line 2260

	/* >>chng 06 jun 06, NP, assert thrown in helike_recom at T == 1e10 K, just bump the 
	 * high temperature end slightly.  */
	TeRRCoef[N_HE_TE_RECOMB-1] += 0.01f;

The comes just before the statement that reads "if( !helike.lgNoRecombInterp )".

Second, edit file heat_sum.c and change the code at lines 445 to 450 from the old code

thermal.heating[1][6] += 
	ionbal.CosRayHeatRate*
	((dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1]) + 
	  dense.xIonDense[ipHELIUM][0])* Cosmic_ray_heat_eff;
/* >>chng 04 mar 24, do not multiply by further sec2prim since heat is 35 eV per primary */
/*	secondaries.sec2prim;*/

to the following:

{
	/* add on cosmic rays - most important in neutral or molecular gas, use
	 * difference between total H and H+ density as substitute for sum of
	 * H in atoms and all molecules, but only in limit where difference is
	 * significant */
	double hneut;
	if( (dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1])/
		dense.gas_phase[ipHYDROGEN]<0.001 )
	{
		/* limit where most H is ionized - simply use atomic and H2 */
		hneut = dense.xIonDense[ipHYDROGEN][0] + 2.*(hmi.Hmolec[ipMH2g]+hmi.Hmolec[ipMH2s]);
	}
	else
	{
		/* limit where neutral is significant, use different between total and H+ */
		hneut = dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1];
	}
	thermal.heating[1][6] += 
		ionbal.CosRayHeatRate*
		/*>>chng 06 jun 26, bug fix, this became unstable at very high ionization 
		 * where cr heating dominated */
		/*((dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1]) + */
		(hneut + dense.xIonDense[ipHELIUM][0])* Cosmic_ray_heat_eff;
	/* >>chng 04 mar 24, do not multiply by further sec2prim since heat is 35 eV per primary */
	/*	secondaries.sec2prim;*/
}

Many thanks to Jennifer West for reporting this bug. 2006 June 07.


The punch FeII lines command only outputs emission from the lower levels in an optically thick cloud. This is because the number of FeII levels is adjusted as the calculation progresses. The number of levels considered at the end of the calculation may be small if the gas becomes cold. To fix this edit atom_feii.c and change the variable FeII.nFeIILevel to FeII.nFeIILevelAlloc at lines 1900 and 1902. After the change the lines will read as follows:

	for( ipLo=0; ipLo < (FeII.nFeIILevelAlloc - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevelAlloc; ipHi++ )
		{

2006 May 20


The interpolated continuum was truncated at a certain cell number when very high resolution continua with very large numbers of points were used. To fix change two files. First edit edit parse_interp.c and change rfield.nupper to NCELL at line 260 and line 263. Before the fix the lines will read a

if( npairs + 2 < rfield.nupper )
{
    /* zero out remainder of array */
    for( i=npairs + 1; i < rfield.nupper; i++ )
    {
        rfield.tNuRyd[rfield.nspec][i] = 0.;
    }
}

after the change this will read as

if( npairs + 2 < NCELL )
{
    /* zero out remainder of array */
    /*>>chng 06 may 17, same as above */
    /*for( i=npairs + 1; i < rfield.nupper; i++ )*/
    for( i=npairs + 1; i < NCELL; i++ )
    {
        rfield.tNuRyd[rfield.nspec][i] = 0.;
    }
}

The second change is in file cont_ffun.c. Change rfield.nupper to NCELL-1 at line 160. Before the fix the lines will read as

/* loop starts at second array element, [1], since want to
 * find next continuum energy greater than desired point */
/*for( i=1; i < rfield.nupper; i++ )*/
i = 1;
while( i< rfield.nupper && rfield.tNuRyd[rfield.ipspec][i]>0. )
{
    if( xnu < rfield.tNuRyd[rfield.ipspec][i] )
    {

after the change these lines will read as

/*>>chng 06 may 17, upper limit of nupper fails when reading in very fine
 * continuum mesh such as that output by stellar atmosphere 
 * bug caught by Christophe Morisset */
/*while( i< rfield.nupper && rfield.tNuRyd[rfield.ipspec][i]>0. )*/
while( i< NCELL-1 && rfield.tNuRyd[rfield.ipspec][i]>0. )
{
	if( xnu < rfield.tNuRyd[rfield.ipspec][i] )
	{

Many thanks to Christophe Morisset for discovering this problem and to Swarna K Ghosh for helping clarify this fix. 2006 May 17


The abundances starburst command did not work with the vary option. If the starburst abundances are to be varied during an optimization run then it will be necessary to change the format used for the abundances scale facror. Edit routine abund_starburst.c and change line 317 from the following:

strcpy( optimize.chVarFmt[optimize.nparm], "ABUNDANCES STARBURST %9.5F LOG" );

to read as follows:

strcpy( optimize.chVarFmt[optimize.nparm], "ABUNDANCES STARBURST %f LOG" );

That is, change %9.5F to %f. 06 May 14


The FeII spectrum produced with the punch FeII continuum command is shifted by about 600 km/s. It was initially only intended to make rough plots of the UV continuum and imprecise physical constants were used. Additionally, vacuum rather than air wavelengths are used for all wavelengths. The first problem can be fixed by editing atom_feii.c and changing the code around line 1500 to read as follows:

		/* convert to line energy in ergs */
		/*ahi = 1.99e-8f/wl1;
		alo = 1.99e-8f/wl2;*/
		/*>>chng 06 apr 10, from above to below, imprecise converstion
		 * factor for Angstroms into ergs, caused shift in spectrum by
		 * about 600 km/s.  Bug caught by Yumihiko Tsuzuki */
		ahi = 1.9864e-8f/wl1;
		alo = 1.9864e-8f/wl2;

The next bug fix roll up with include this fix and also use air wavelengths for wl > 2000A This is a very small effect and so not done here as a hot fix. Many thanks to Yumihiko Tsuzuki for reporting this bug. 2006 April 10.


fixes to version 06.02

The atom H2 dissociation collisions command turned off all collisions within the large H2 molecule. Edit routine parse_atomh2.c and add an "else" before the second if within the collisions branch at line 142. The changed code will read as follows:

	/* option to turn collisional dissociation off or on 
	 * >>chng 06 mar 01, had been simply if - so all collisions were turned off
	 * when dissociation collisions turned off - 
	 * due to bucket else at end */
	else if( lgMatch("ORTH",chCard ) && lgMatch("PARA",chCard ) )
	{
		/* option to turn ortho - para collisions with particles off */

2006 March 1


Approximations to Nussbaumer & Storey DR rates could diverge for very low temperatures. The approximations in their published papers were not mean to be used outside certain temperature limits. These rates are still the default but are being replaced by the Badnell values. One test on the range of validity of the NS rates was left out, and the rates could strongly diverge for very low temperatures. Edit routine ion_recom.c near line 214 and add a test on the value of tefac. After the change the code will read as follows (code before and after the change is included for reference):

	if( ff[ion] != 0. && nelem!=ipIRON )
	{
		/* >>chng 06 feb 14, O+3 ff[ion] is very negative, as a result exp goes to +inf
		 * at very low tempertatures.  This is a error in the Nussbaumer & Storey fits to DR.
		 * do not use them is tefac = ff[ion] / t4 is very negative */
		if( tefac > -5. )
		{
			factor = (((aa[ion]*t4m1+bb[ion])*t4m1+cc[ion])*t4m1+dd[ion])* sexp(tefac);
			DielRecomRateCoef_LowT[ion] = ionbal.DielSupprs[1][ion]*fac2*
				MAX2(0.,factor );
		}
		else
		{
			DielRecomRateCoef_LowT[ion] = 0.;
		}
	}

2006 Feb 27


The wavelength used to identify the continuum bands (given in bands_continuum.dat and entered into the emission line printout) gave the wavelength in Angstroms when it should have been in microns. Edit routine cont_createpointers.c and add the following code to line1809 (code before and after the change is also shown)

			continuum.ContBandWavelength[k] = (float)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			/* >>chng 06 feb 21, multiply by 1e4 to convert miron wavelength into Angstroms,
			 * which is assumed by the code.  before this correction the band centroid 
			 * wavelenth was given in the output incorrectly listed as Angstroms.
			 * results were correct just label was wrong */
			continuum.ContBandWavelength[k] *= 1e4f;

			/* these are short and long wave limits, which are high and
			 * low energy limits - these are now wl in microns */

2006 Feb 27


FeII may fail with insanity under some conditions when more than one iteration is performed. To fix this edit routine atom_feii.c and add two lines, defining and setting the variable nFeII_old, and another line testing its value. After these changes the code between lines 1203 and 1222 will read as follows (snippet begins and ends with neighboring lines of code to place in context):

	long int ipTemp,
		ipTempp1 , 
		ipLoFe2, 
		ipHiFe2;
	static long int nFeII_old = -1;

#	ifdef DEBUG_FUN
	fputs( "<+>FeIICollRatesBoltzmann()\n", debug_fp );
#	endif

	/* return if temperature has not changed */
	/* >>chng 06 feb 14, add test on number of levels changing */
	/*lint -e777 test floats for equality*/
	if( phycon.te == OldTemp && FeII.nFeIILevel== nFeII_old )
	{
#		ifdef DEBUG_FUN
		fputs( " <->FeIICollRatesBoltzmann()\n", debug_fp );
#		endif

		return;
	}
	/*lint +e777 test floats for equality*/
	OldTemp = phycon.te;
	nFeII_old = FeII.nFeIILevel;

	for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ )

2006 Feb 27


The code will generate a bogus warning when the density is set by an external function, such as fabden or dlaw. It will complain that the density is changing by too much, even when it is constant. The warning is benign and can simply be ignored. To fix the code, edit routine pressure_change.c and change the code near line 572 to read as follows:

		pressure.PresTotlCorrect = pressure.PresTotlCurr*pressure_change_factor;
		/* >>chng 06 feb 21, from above to below.  we want to look at change,
		 * not the value.  Bug found by Valentina Luridiana */
		/*if( pressure.PresTotlCorrect > 3. || pressure.PresTotlCorrect< 1./3 )*/
		if( pressure_change_factor > 3. || pressure_change_factor< 1./3 )
		{
			static int lgWARN2BIG=FALSE;

Many thanks to Valentina Luridiana for finding this problem. 2006 Feb 23.


Return to HotFixes main page

Return to the StepByStep instructions

Return to  nublado.org