Skip to content

Commit

Permalink
Bug511 ff (sirocco-rt#512)
Browse files Browse the repository at this point in the history 
* Removed emrng structure, which was causing problems in restarts of macro-atom models 

* Modifications to one_continuum to assure that  the wavelength array created for
cdf_gen_from_array does not have two points at the end of the array with the same
wavelength.  See sirocco-rt#511.
  • Loading branch information
@kslong
kslong authored May 6, 2019
1 parent 718f3bf commit 9ebbcc7
Show file tree
Hide file tree
Showing 8 changed files with 66 additions and 68 deletions.
61 changes: 36 additions & 25 deletions source/cdf.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,7 +123,7 @@ cdf_gen_from_func (cdf, func, xmin, xmax, njumps, jump)
if (xmax <= xmin)
{
Error ("pdf_gen_from_func: xmin %g <= xmax %g\n", xmin, xmax);
Exit (0);
Exit (1);
}
if (njumps > 0)
{
Expand All @@ -132,7 +132,7 @@ cdf_gen_from_func (cdf, func, xmin, xmax, njumps, jump)
if (jump[j] <= jump[j - 1])
{
Error ("pdf_gen_from_func: jump[%d]=%g <=jump[%d]=%g out of order\n", j, jump[j], j - 1, jump[j - 1]);
Exit (0);
Exit (1);
}
}
njump_min = 0;
Expand Down Expand Up @@ -293,7 +293,7 @@ gen_array_from_func (func, xmin, xmax, pdfsteps)
if ((pdf_array = calloc (sizeof (x), pdfsteps)) == NULL)
{
Error ("pdf: Could not allocate space for pdf_array\n");
Exit (0);
Exit (1);
}
init_pdf = 1;
pdf_steps_current = pdfsteps;
Expand Down Expand Up @@ -369,13 +369,13 @@ int pdf_n;
/**
* @brief Generate a cumulative distribution function (cdf) from an array of values
*
* @param [in] cdf Pointer to the cdf that will be populated here
* @param [in] x Locations of points in supplied data - usually wavelength or freuency
* @param [in] y Data - e.g. cross section or emissivity
* @param [in] n_xy Number of data points
* @param [in] xmin Lower limit to required cdf
* @param [in] xmin Upper limit to required cdf
* @return echeck the reutrn value from cdf_check
* @param [in] cdf Pointer to the cdf that will be populated here
* @param [in] x Locations of points in supplied data - usually wavelength or freuency
* @param [in] y Data - e.g. cross section or emissivity
* @param [in] n_xy Number of data points
* @param [in] xmin Lower limit to required cdf
* @param [in] xmin Upper limit to required cdf
* @return echeck the return value from cdf_check
*
* @details
* Takes two arrays of some values of a function y at locations
Expand All @@ -387,10 +387,12 @@ int pdf_n;
*
* cdf_gen_from_array does not have the concept of jumps. All of
* this needs to be taken care of by the routine that generates
* the array.
* the array. This means the x's should be monotonic, withoug
* duplicates, and all the values of y should be 0 or positive.
* The program exits if these conditions are not fulfilled.
*
* In constructing the cdf the routine collapses regions of the
* pdf that have zero probability of occuring
* pdf that have zero probability of occuring.
*
**********************************************************/

Expand All @@ -403,7 +405,7 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax)
{
int allzero;
int nmin, nmax, cdf_n;
int n;
int n, nn;
double sum;
int echeck, zcheck = 0;

Expand All @@ -414,14 +416,14 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax)
if (n_xy > NCDF) //If the supplied data array is larger than the output array - fixed as NCDF
{
Error ("cdf_gen_from_array: supplied data %i is larger than default aray size %i - increase NCDF\n", n_xy, NCDF);
Exit (0);
Exit (1);
}


if (xmax < xmin) //This must be a mistake, the limits are reversed
{
Error ("cdf_gen_from_array: xmin %g <= xmax %g\n", xmin, xmax);
Exit (0);
Exit (1);
}

/*We are now going to crawl down the input array, checking for mistakes.
Expand All @@ -435,13 +437,22 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax)
{
if (x[n] <= x[n - 1]) //One x-point is less than the one above - this would cause problems later
{
for (nn = 0; nn < n_xy; nn++)
{
Log ("cdf_gen_from_array: %5d %11.6e %11.6e\n", nn, x[nn], y[nn]);
}
Error ("cdf_gen_from_array: input x not in ascending order at element %5d/%5d %11.6e %11.6e\n", n, n_xy, x[n - 1], x[n]);
Exit (0);

Exit (1);
}
if (y[n] < 0) //A negative point!
{
for (nn = 0; nn < n_xy; nn++)
{
Log ("cdf_gen_from_array: %5d %11.6e %11.6e\n", nn, x[nn], y[nn]);
}
Error ("cdf_gen_from_array: probability density %g < 0 at element %d\n", y[n], n);
Exit (0);
Exit (1);
}
else if (y[n] > 0) //At least one point is positive
{
Expand All @@ -451,13 +462,13 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax)



/* ksl - Our requrements are as follows. We want the cdf to have boundaries given by xmin and xmax, unless
/* ksl - Our requirements are as follows. We want the cdf to have boundaries given by xmin and xmax, unless
* the array does not encompass this, in which case we will only include what the array contains. We also
* want to suppress regions at the end of the cdf where the pdf was 0
*/


/*Firstly, find the point in the input array that matches the required start point
/*First, find the point in the input array that matches the required start point
we step up through the input data until we reach a point where the x-point is greater
than xmin, or we reach the end of the array */

Expand All @@ -477,7 +488,7 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax)
if (nmin == n_xy)
{
Error ("cdf_gen_from_array: xmin greater than all array values\n");
Exit (0);
Exit (1);
}

/* if xmin is equal to one of the values in the x array then nmin will point to that value, but otherwise it
Expand All @@ -501,7 +512,7 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax)
if (nmax == nmin)
{
Error ("cdf_gen_from_array: nmin and nmax are identical which is not desirable\n");
Exit (0);
Exit (1);
}

if (nmin == 0 && xmin < x[0]) /*We are requesting a CDF that starts below where we have data
Expand Down Expand Up @@ -552,7 +563,7 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax)
if (nmax == nmin)
{
Error ("cdf_gen_from_array - modified nmax=nmin followin interpolation at ends\n");
Exit (0);
Exit (1);
}


Expand Down Expand Up @@ -610,7 +621,7 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax)
{
Error ("cdf_gen_from_array: error %d on cdf_check - check CDF_err.diag\n", echeck);
cdf_to_file (cdf, "CDF_err.diag"); //output the CDF to a file
Exit (0);
Exit (1);
}
// cdf_to_file(cdf,"foo.diag"); //output the CDF to a file

Expand Down Expand Up @@ -729,7 +740,7 @@ cdf_limit (cdf, xmin, xmax)
if (cdf->y[cdf->ncdf] != 1.0)
{
Error ("cdf_limit: cdf not defined!)");
Exit (0);
Exit (1);
}
if (xmin >= cdf->x[cdf->ncdf])
{
Expand All @@ -738,7 +749,7 @@ cdf_limit (cdf, xmin, xmax)
if (xmax <= cdf->x[0])
{
Error ("cdf_limit: xmax %g < cdf->x[0] %g\n", xmax, cdf->x[0]);
Exit (0);
Exit (1);
}

/* Set the limits for the minimum */
Expand Down
15 changes: 11 additions & 4 deletions source/continuum.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,8 @@ one_continuum (spectype, t, g, freqmin, freqmax)
lambdamax = C * 1e8 / freqmin;
nwave = 0;

/* if the first wavelength in the model is below the wavelength range in the simulation,
/* Create the first element of the array, precisely at lambdamin if that is possible
Specifically, if the first wavelength in the model is below the wavelength range in the simulation,
interpolate on the model flux to get the flux at lambdamin. copy relevant wavelengths and
fluxes to w_local and f_local */
if (comp[spectype].xmod.w[0] < lambdamin && lambdamin < comp[spectype].xmod.w[comp[spectype].nwaves - 1])
Expand All @@ -150,19 +151,25 @@ one_continuum (spectype, t, g, freqmin, freqmax)
nwave++;
}

/* loop over rest of model wavelengths and fluxes and copy to w_local and f_local */
/* loop over rest of model wavelengths and fluxes and copy to w_local and f_local.
This does not include the end points
*/

for (n = 0; n < comp[spectype].nwaves; n++)
{
if (comp[spectype].xmod.w[n] > lambdamin && comp[spectype].xmod.w[n] <= lambdamax)
//OLD if (comp[spectype].xmod.w[n] > lambdamin && comp[spectype].xmod.w[n] <= lambdamax)
if (comp[spectype].xmod.w[n] > lambdamin && comp[spectype].xmod.w[n] < lambdamax)
{
w_local[nwave] = comp[spectype].xmod.w[n];
f_local[nwave] = comp[spectype].xmod.f[n];
nwave++;
}
}

/* now check if upper bound is beyond lambdamax, and if so, interpolate to get appropriate flux
/* No add a point at lambdamax. Specicxally check if upper bound is beyond lambdamax, and if so,
interpolate to get appropriate flux
at lambda max. copy to w_local and f_local */

if (comp[spectype].xmod.w[0] < lambdamax && lambdamax < comp[spectype].xmod.w[comp[spectype].nwaves - 1])
{
w_local[nwave] = lambdamax;
Expand Down
10 changes: 5 additions & 5 deletions source/matom.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -814,8 +814,8 @@ kpkt (p, nres, escape, mode)
else
{
/* in spectral cycles, use the frequency range of the final spectrum */
freqmin = em_rnge.fmin;
freqmax = em_rnge.fmax;
freqmin = geo.sfmin;
freqmax = geo.sfmax;
}


Expand Down Expand Up @@ -1521,7 +1521,7 @@ emit_matom (w, p, nres, upper)
line_ptr = &line[config[uplvl].bbd_jump[n]];
/* Since we are only interested in making an r-packet here we can (a) ignore collisional
deactivation and (b) ignore lines outside the frequency range of interest. */
if ((line_ptr->freq > em_rnge.fmin) && (line_ptr->freq < em_rnge.fmax)) // correct range
if ((line_ptr->freq > geo.sfmin) && (line_ptr->freq < geo.sfmax)) // correct range
{
bb_cont = (a21 (line_ptr) * p_escape (line_ptr, xplasma));

Expand All @@ -1545,7 +1545,7 @@ emit_matom (w, p, nres, upper)

/* If the edge is above the frequency range we are interested in then we need not consider this
bf process. */
if (cont_ptr->freq[0] < em_rnge.fmax) //means that it may contribute
if (cont_ptr->freq[0] < geo.sfmax) //means that it may contribute
{
sp_rec_rate = alpha_sp (cont_ptr, xplasma, 0);
eprbs[m] = sp_rec_rate * ne * (config[uplvl].ex - config[phot_top[config[uplvl].bfd_jump[n]].nlev].ex); //energy difference
Expand Down Expand Up @@ -1687,7 +1687,7 @@ matom_emit_in_line_prob (WindPtr one, struct lines *line_ptr_emit)

if (eprbs_line == 0.0)
{
Error ("matom_emit_in_line_prob: Line frequency %g lies outside spectral range %g-%g!\n", line_ptr->freq, em_rnge.fmin, em_rnge.fmax);
Error ("matom_emit_in_line_prob: Line frequency %g lies outside spectral range %g-%g!\n", line_ptr->freq, geo.sfmin, geo.sfmax);
return (-1.0);
}

Expand Down
5 changes: 2 additions & 3 deletions source/models.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,7 @@
/* 1405 JM -- Increased PDF array for use with disk14 models- also removed duplication of ncomps */
/* 1409 JM -- Increased LINELEN to 160 */

#define NWAVES 28000
// #define NWAVES 2000
#define NWAVES 28000 // The maximum number of wavelength bins in models
#define NDIM 10 // The maximum number of free parameters
#define NCOMPS 10 //The maximum number of separate components
#define NPARS 10 //The maximum number of parameters in a component (not all variable)
Expand Down Expand Up @@ -48,7 +47,7 @@ struct ModSum
char name[LINELEN];
int npars; // Number of parameters that describe the model
int nmods, modstart, modstop; //number of models of this type and the index into the model struct
double min[NPARS]; // The minimum and maximum for each "free" paratmenter of the model
double min[NPARS]; // The minimum and maximum for each "free" parameter of the model
double max[NPARS];
int nwaves; //All models in each comp should have same wavelengths;
struct Model xmod; //The current intepolated model of this type
Expand Down
8 changes: 4 additions & 4 deletions source/photo_gen_matom.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -364,7 +364,7 @@ get_matom_f (mode)
}


if (ppp.freq > em_rnge.fmin && ppp.freq < em_rnge.fmax && ppp.istat != P_ADIABATIC)
if (ppp.freq > geo.sfmin && ppp.freq < geo.sfmax && ppp.istat != P_ADIABATIC)
{
if (which_out == 1)
{
Expand Down Expand Up @@ -567,8 +567,8 @@ photo_gen_kpkt (p, weight, photstart, nphot)
else
{
/* we are in the spectral cycles, so use all the required frequency range */
fmin = em_rnge.fmin;
fmax = em_rnge.fmax;
fmin = geo.sfmin;
fmax = geo.sfmax;
/* we only want k->r processes */
kpkt_mode = KPKT_MODE_CONTINUUM;
}
Expand Down Expand Up @@ -790,7 +790,7 @@ photo_gen_matom (p, weight, photstart, nphot)

test = pp.freq;

while (test > em_rnge.fmax || test < em_rnge.fmin)
while (test > geo.sfmax || test < geo.sfmin)
{

/* Call routine that will select an emission process for the
Expand Down
2 changes: 1 addition & 1 deletion source/python.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -651,7 +651,7 @@ main (argc, argv)
* radiation and then value given to return a spectrum type. The output is not the same
* number as one inputs. It' s not obvious that this is a good idea. */

if (geo.pcycles > 0 && geo.run_type!=RUN_TYPE_RESTART)
if (geo.pcycles > 0 && geo.run_type != RUN_TYPE_RESTART)
{

rdpar_comment ("Parameters defining the spectra seen by observers\n");
Expand Down
Loading

0 comments on commit 9ebbcc7

Please sign in to comment.