diff -ur c07.02.01/source/cont_setintensity.cpp c07.02.01.new/source/cont_setintensity.cpp
--- c07.02.01/source/cont_setintensity.cpp 2007-07-09 04:28:29.000000000 +0200
+++ c07.02.01.new/source/cont_setintensity.cpp 2008-02-14 12:27:26.000000000 +0100
@@ -199,6 +199,13 @@
/*fprintf(ioQQQ,"DEBUG %4li %e %e %.3e %.3e\n",
i , rfield.anu[i] , (amean2/amean) , amean2 , amean );*/
rfield.anu[i] = (float)(amean2/amean);
+ if( i )
+ {
+ /* prevent roundoff from allowing i cell to lie below i-1
+ * * cell when continuum mesh is very fine */
+ rfield.anu[i] = MAX2( rfield.anu[i] , rfield.anu[i-1]*(1.f+FLT_EPSILON) );
+ ASSERT( rfield.anu[i] > rfield.anu[i-1] );
+ }
rfield.anu2[i] = (float)(amean3/amean);
}
diff -ur c07.02.01/source/grid_xspec.cpp c07.02.01.new/source/grid_xspec.cpp
--- c07.02.01/source/grid_xspec.cpp 2007-07-09 04:28:29.000000000 +0200
+++ c07.02.01.new/source/grid_xspec.cpp 2008-02-14 12:34:06.000000000 +0100
@@ -72,7 +72,7 @@
sprintf( grid.paramNames[i], "%s%ld", "PARAM", i+1 );
/* Method is 0 for linear, 1 for logarithmic, for now all linear. */
- grid.paramMethods[i] = 1;
+ grid.paramMethods[i] = 0;
/* Initial */
grid.paramRange[i][0] = xc[i]+grid.paramIncrements[i]*(grid.numParamValues[i]-1.f)/2.f;
/* Delta */
diff -ur c07.02.01/source/iso_create.cpp c07.02.01.new/source/iso_create.cpp
--- c07.02.01/source/iso_create.cpp 2007-07-09 04:28:29.000000000 +0200
+++ c07.02.01.new/source/iso_create.cpp 2008-02-14 12:31:27.000000000 +0100
@@ -364,7 +364,7 @@
}
else
{
- jj = jj;
+ jj = j;
}
hydro.strkar[n][j] = (float)pow(((float)ii*jj),1.2f);
hydro.pestrk[n][j] = 0.;
diff -ur c07.02.01/source/opacity.h c07.02.01.new/source/opacity.h
--- c07.02.01/source/opacity.h 2007-07-09 04:28:29.000000000 +0200
+++ c07.02.01.new/source/opacity.h 2008-02-14 12:40:08.000000000 +0100
@@ -172,7 +172,7 @@
/** another case b option, turn off background opacities, no Pdest */
bool lgCaseB_no_pdest;
- /** optical depths to compton and H- */
+ /** optical depths to Compton and H- */
float telec,
thmin;
@@ -203,9 +203,6 @@
* depths have been predicted yet */
bool lgTauOutOn;
- /** which opacity to punch out */
- char chOpcTyp[5];
-
/** total number of opacity cells used in opacity stack
* in OpacityCreateAll used as a counter to remember where
* next opacity goes */
@@ -217,7 +214,7 @@
* in OpacityCreateAll
*/
- /** ipRayScat opacity pointer for rayleigh scattering*/
+ /** ipRayScat opacity pointer for Rayleigh scattering*/
long int ipRayScat,
/** iopcom compton scatterin, total recoil*/
diff -ur c07.02.01/source/parse_grain.cpp c07.02.01.new/source/parse_grain.cpp
--- c07.02.01/source/parse_grain.cpp 2007-07-09 04:28:29.000000000 +0200
+++ c07.02.01.new/source/parse_grain.cpp 2008-02-14 12:31:02.000000000 +0100
@@ -364,8 +364,8 @@
/* find out how long chOption is - it was declared as a pointer, then
* set equal to a local string */
length = strlen( chOption );
- if( length == 0 )
- TotalInsanity();
+ /*if( length == 0 )
+ TotalInsanity();*/
/* make sure we have enough room to store the string */
if( length > static_cast<size_t>(FILENAME_PATH_LENGTH_2) )
diff -ur c07.02.01/source/parse_punch.cpp c07.02.01.new/source/parse_punch.cpp
--- c07.02.01/source/parse_punch.cpp 2007-07-09 04:28:29.000000000 +0200
+++ c07.02.01.new/source/parse_punch.cpp 2008-02-14 12:40:08.000000000 +0100
@@ -239,7 +239,7 @@
{
/* DoPunch will call punch_opacity to parse the subcommands
* punch total opacity */
- strcpy( opac.chOpcTyp, "TOTL" );
+ strcpy( punch.chOpcTyp[punch.npunch], "TOTL" );
sprintf( chHeader,
"#nu\tTot opac\tAbs opac\tScat opac\tAlbedo\telem\n" );
}
@@ -247,24 +247,25 @@
else if( nMatch("FIGU",chCard) )
{
/* do figure for hazy */
- strcpy( opac.chOpcTyp, "FIGU" );
+ strcpy( punch.chOpcTyp[punch.npunch], "FIGU" );
sprintf( chHeader,
"#nu, H, He, tot opac\n" );
}
else if( nMatch("FINE",chCard) )
{
- /* punch out the fine opacity array with only lines */
+ /* punch the fine opacity array */
rfield.lgPunchOpacityFine = true;
- strcpy( opac.chOpcTyp, "FINE" );
+ strcpy( punch.chOpcTyp[punch.npunch], "FINE" );
/* check for keyword UNITS on line, then scan wavelength or energy units if present,
- * units are copied into punch.chConPunEnr */
+ * units are copied into punch.chConPunEnr */
ChkUnits(chCard);
sprintf( chHeader,
"#nu\topac\n" );
ipFFmt = 5;
- /* range option - important since so much data */
+ /* range option - important since so much data - usually want to
+ * only give portion of the continuum */
if( nMatch("RANGE",chCard) )
{
/* get lower and upper range, must be in Ryd */
@@ -296,6 +297,7 @@
/* default is to average together ten */
if( lgEOL )
punch.punarg[punch.npunch][2] = 10;
+
if( punch.punarg[punch.npunch][2] < 1 )
{
fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n");
@@ -316,7 +318,7 @@
else if( nMatch("BREM",chCard) )
{
/* punch bremsstrahlung opacity */
- strcpy( opac.chOpcTyp, "BREM" );
+ strcpy( punch.chOpcTyp[punch.npunch], "BREM" );
sprintf( chHeader,
"#nu\tbrem opac\n" );
}
@@ -327,7 +329,7 @@
strcpy( punch.chPunch[punch.npunch], "OPAC" );
/* punch subshell cross sections */
- strcpy( opac.chOpcTyp, "SHEL" );
+ strcpy( punch.chOpcTyp[punch.npunch], "SHEL" );
/* this is element */
ipFFmt = 3;
@@ -370,7 +372,7 @@
}
/* copy string over */
- strcpy( opac.chOpcTyp, elementnames.chElementNameShort[nelem] );
+ strcpy( punch.chOpcTyp[punch.npunch], elementnames.chElementNameShort[nelem] );
}
else
{
@@ -437,7 +439,7 @@
else if( nMatch("OPAC",chCard) )
{
/* create table for appendix in AGN */
- strcpy( opac.chOpcTyp, " AGN" );
+ strcpy( punch.chOpcTyp[punch.npunch], " AGN" );
strcpy( punch.chPunch[punch.npunch], "OPAC" );
}
diff -ur c07.02.01/source/parse_set.cpp c07.02.01.new/source/parse_set.cpp
--- c07.02.01/source/parse_set.cpp 2007-07-09 04:28:29.000000000 +0200
+++ c07.02.01.new/source/parse_set.cpp 2008-02-14 12:33:46.000000000 +0100
@@ -109,7 +109,7 @@
{
if( nMatch( " NEW" , chCard ) )
{
- dense.lgAsChoose[nelem][ion] = true;
+ dense.lgAsChoose[nelem][ion-1] = true;
}
else if( nMatch( " OLD" , chCard ) )
{
diff -ur c07.02.01/source/prt_lines_lv1_li_ne.cpp c07.02.01.new/source/prt_lines_lv1_li_ne.cpp
--- c07.02.01/source/prt_lines_lv1_li_ne.cpp 2007-07-09 04:28:29.000000000 +0200
+++ c07.02.01.new/source/prt_lines_lv1_li_ne.cpp 2008-02-14 12:28:26.000000000 +0100
@@ -912,6 +912,8 @@
/* recombination and specific pump for OII 4638.86-4696.35 (8 lines) */
rec = GetLineRec(82, 4651 );
+ PntForLine(4651.,"O 2",&ipnt);
+ lindst(rec,4651,"O 2r",ipnt,'i',true );
/* convert UV pump rate to intensity with branching ratio and hnu */
/* recombination part of O II 4651 line */
linadd(rec,4651,"O 2r",'i' );
diff -ur c07.02.01/source/punch.h c07.02.01.new/source/punch.h
--- c07.02.01/source/punch.h 2007-07-09 04:28:29.000000000 +0200
+++ c07.02.01.new/source/punch.h 2008-02-14 12:40:08.000000000 +0100
@@ -147,6 +147,9 @@
/**chPunch - what is it we want to punch? set in GetPunch, used in DoPunch */
char chPunch[LIMPUN][5];
+ /** which opacity to punch out */
+ char chOpcTyp[LIMPUN][5];
+
/** this flag tells us whether to punch results of a grid to separate files
* for each grid point or all to the same file. Different for different
* punch commands */
diff -ur c07.02.01/source/punch_opacity.cpp c07.02.01.new/source/punch_opacity.cpp
--- c07.02.01/source/punch_opacity.cpp 2007-07-09 04:28:29.000000000 +0200
+++ c07.02.01.new/source/punch_opacity.cpp 2008-02-14 12:40:08.000000000 +0100
@@ -49,7 +49,7 @@
/* punch total opacity in any element, punch opacity command
* ioPUN is ioPUN unit number, ipPun is pointer within stack of punches */
- if( strcmp(opac.chOpcTyp,"TOTL") == 0 )
+ if( strcmp(punch.chOpcTyp[ipPun],"TOTL") == 0 )
/* total opacity */
{
for( j=0; j < rfield.nflux; j++ )
@@ -66,8 +66,8 @@
fprintf( ioPUN, "\n" );
}
- else if( strcmp(opac.chOpcTyp,"BREM") == 0 )
- /* brems opacity */
+ else if( strcmp(punch.chOpcTyp[ipPun],"BREM") == 0 )
+ /* bremsstrahlung opacity */
{
for( j=0; j < rfield.nflux; j++ )
{
@@ -80,7 +80,7 @@
}
/* subshell photo cross sections */
- else if( strcmp(opac.chOpcTyp,"SHEL") == 0 )
+ else if( strcmp(punch.chOpcTyp[ipPun],"SHEL") == 0 )
{
nelem = (long)punch.punarg[ipPun][0];
ion = (long)punch.punarg[ipPun][1];
@@ -95,7 +95,7 @@
}
}
- else if( strcmp(opac.chOpcTyp,"FINE") == 0 )
+ else if( strcmp(punch.chOpcTyp[ipPun],"FINE") == 0 )
{
/* the fine opacity array */
float sum;
@@ -119,7 +119,9 @@
{
nu_hi = rfield.nfine;
}
+
nskip = (long)punch.punarg[ipPun][2];
+ ASSERT( nskip > 0 );
for( i=j; i<nu_hi; i+=nskip )
{
@@ -131,12 +133,12 @@
sum /= nskip;
if( sum > 0. )
fprintf(ioPUN,"%.5e\t%.3e\n",
- AnuUnit( rfield.fine_anu[i+k/2] ), sum );
+ AnuUnit( rfield.fine_anu[i+k/2] ), sum );
}
}
/* figure for hazy */
- else if( strcmp(opac.chOpcTyp,"FIGU") == 0 )
+ else if( strcmp(punch.chOpcTyp[ipPun],"FIGU") == 0 )
{
nelem = 0;
for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1; i < (iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] - 1); i++ )
@@ -166,7 +168,7 @@
}
/* photoionization data table for AGN */
- else if( strcmp(opac.chOpcTyp," AGN") == 0 )
+ else if( strcmp(punch.chOpcTyp[ipPun]," AGN") == 0 )
{
long int
ipop,
@@ -234,7 +236,7 @@
}
/* hydrogen */
- else if( strcmp(opac.chOpcTyp,"HYDR") == 0 )
+ else if( strcmp(punch.chOpcTyp[ipPun],"HYDR") == 0 )
{
nelem = ipHYDROGEN;
/* zero out the opacity arrays */
@@ -278,7 +280,7 @@
}
/* helium */
- else if( strcmp(opac.chOpcTyp,"HELI") == 0 )
+ else if( strcmp(punch.chOpcTyp[ipPun],"HELI") == 0 )
{
/* atomic helium first, HELI1.opc */
nelem = ipHELIUM;
@@ -356,12 +358,12 @@
else
{
- /* check for hydroge through zinc, nelem is atomic number on the c scale */
+ /* check for hydrogen through zinc, nelem is atomic number on the c scale */
nelem = -1;
i = 0;
while( i < LIMELM )
{
- if( strcmp(opac.chOpcTyp,elementnames.chElementNameShort[i]) == 0 )
+ if( strcmp(punch.chOpcTyp[ipPun],elementnames.chElementNameShort[i]) == 0 )
{
nelem = i;
break;
@@ -373,7 +375,7 @@
if( nelem < 0 )
{
fprintf( ioQQQ, " Unidentified opacity key=%4.4s\n",
- opac.chOpcTyp );
+ punch.chOpcTyp[ipPun] );
puts( "[Stop in punopac]" );
cdEXIT(EXIT_FAILURE);
}
diff -ur c07.02.01/source/zero.cpp c07.02.01.new/source/zero.cpp
--- c07.02.01/source/zero.cpp 2007-07-09 04:28:29.000000000 +0200
+++ c07.02.01.new/source/zero.cpp 2008-02-14 12:40:08.000000000 +0100
@@ -1011,10 +1011,8 @@
/* another case b option, turn off background opacities, no Pdest */
opac.lgCaseB_no_pdest = false;
- strcpy( opac.chOpcTyp, " " );
-
/* smallest allowed line and Lya optical depths, reset with
- * caseb command */
+ * Case B command */
opac.taumin = 1e-20f;
opac.tlamin = 1e-20f;