commit c36c173300a96a912eefcd1858e5a16211824f44
parent e8b44e9a3d396694530044a7d8e8b003771fbde7
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 24 Nov 2023 10:55:29 +0100
Changing the wavelength unit of the Planck distribution
The unit used to be nanometers, whereas the rest of the API uses meters.
This commit unifies the units: all wavelengths are now defined in
meters.
Diffstat:
2 files changed, 28 insertions(+), 38 deletions(-)
diff --git a/src/sbb.h b/src/sbb.h
@@ -46,7 +46,7 @@ struct mem_allocator;
struct sbb_ran_planck_create_args {
struct mem_allocator* allocator; /* NULL <=> use default allocator */
struct logger* logger; /* NULL <=> user default logger */
- double range[2]; /* Spectral range [nm] */
+ double range[2]; /* Spectral range [m] */
double ref_temperature; /* Reference temperature [K] */
/* Indication of the number of bands used to discretize the CDF. These bands
@@ -59,7 +59,7 @@ struct sbb_ran_planck_create_args {
#define SBB_RAN_PLANCK_CREATE_ARGS_DEFAULT__ {\
NULL, /* Allocator */\
NULL, /* Logger */\
- {9600, 11500}, /* Spectral range [nm] */\
+ {9.6e-6, 11.5e-6}, /* Spectral range [m] */\
300, /* Reference temperature [K] */\
1900, /* #bands used to speed up sampling, actually 1 band per nanometers */\
0 /* Verbosity */\
@@ -139,12 +139,12 @@ SBB_API res_T
sbb_ran_planck_ref_put
(struct sbb_ran_planck* planck);
-SBB_API double /* Wavelength [nm] */
+SBB_API double /* Wavelength [m] */
sbb_ran_planck_sample
(struct sbb_ran_planck* planck,
const double r0, /* Random number in [0, 1[ */
const double r1, /* Random number in [0, 1[ */
- double* pdf); /* [nm^-1]. May be NULL */
+ double* pdf); /* [m^-1]. May be NULL */
END_DECLS
diff --git a/src/sbb_ran_planck.c b/src/sbb_ran_planck.c
@@ -28,8 +28,8 @@
struct sbb_ran_planck {
struct darray_double pdf;
struct darray_double cdf;
- double range[2]; /* Boundaries of the spectral integration interval [nm] */
- double band_len; /* Length of a band [nm] */
+ double range[2]; /* Boundaries of the spectral integration interval [m] */
+ double band_len; /* Length of a band [m] */
double ref_temperature; /* Reference temperature [K] */
size_t nbands; /* #bands */
@@ -57,7 +57,7 @@ check_ranst_planck_create_args(struct sbb_ran_planck_create_args* args)
|| args->range[0] > args->range[1]) {
if(args->verbose) {
logger_print(logger, LOG_ERROR, MSG_ERROR_PREFIX
- "Invalid spectral range for the Planck distribution: [%g,%g] nm\n",
+ "Invalid spectral range for the Planck distribution: [%g,%g] m\n",
args->range[0], args->range[1]);
}
return RES_BAD_ARG;
@@ -113,10 +113,6 @@ setup_cdf(struct sbb_ran_planck* planck)
ASSERT(lo > planck->range[0] || eq_eps(lo, planck->range[0], 1.e-6));
ASSERT(lo < planck->range[1] || eq_eps(lo, planck->range[1], 1.e-6));
- /* Convert from nanometer to meter */
- lo *= 1.e-9;
- hi *= 1.e-9;
-
/* Compute the probability of the current band */
pdf[i] = sbb_blackbody_fraction(lo, hi, planck->ref_temperature);
@@ -152,7 +148,7 @@ static double
sample_continue
(const struct sbb_ran_planck* planck,
const double r, /* In [0, 1[ */
- const double range[2], /* [nm]*/
+ const double range[2], /* [m]*/
double* pdf) /* May be NULL */
{
/* Numerical parameters of the dichotomy search */
@@ -164,11 +160,10 @@ sample_continue
double bf = 0;
double bf_prev = 0;
double bf_min_max = 0;
- double lambda_m = 0;
- double lambda_m_prev = 0;
- double lambda_m_min = 0;
- double lambda_m_max = 0;
- double range_m[2] = {0, 0};
+ double lambda = 0;
+ double lambda_prev = 0;
+ double lambda_min = 0;
+ double lambda_max = 0;
size_t i;
/* Check precondition */
@@ -176,42 +171,38 @@ sample_continue
ASSERT(range[0] < range[1] && range[0] >= 0 && range[1] >= 0);
ASSERT(0 <= r && r < 1);
- /* Convert the wavelength range in meters */
- range_m[0] = range[0] * 1.0e-9;
- range_m[1] = range[1] * 1.0e-9;
-
/* Setup the dichotomy search */
- lambda_m_min = range_m[0];
- lambda_m_max = range_m[1];
+ lambda_min = range[0];
+ lambda_max = range[1];
bf_min_max = sbb_blackbody_fraction
- (range_m[0], range_m[1], planck->ref_temperature);
+ (range[0], range[1], planck->ref_temperature);
/* Numerically search the lambda corresponding to the submitted canonical
* number */
FOR_EACH(i, 0, MAX_ITER) {
double r_test;
- lambda_m = (lambda_m_min + lambda_m_max) * 0.5;
+ lambda = (lambda_min + lambda_max) * 0.5;
bf = sbb_blackbody_fraction
- (range_m[0], lambda_m, planck->ref_temperature);
+ (range[0], lambda, planck->ref_temperature);
r_test = bf / bf_min_max;
if(r_test < r) {
- lambda_m_min = lambda_m;
+ lambda_min = lambda;
} else {
- lambda_m_max = lambda_m;
+ lambda_max = lambda;
}
- if(fabs(lambda_m_prev - lambda_m) < EPSILON_LAMBDA_M
+ if(fabs(lambda_prev - lambda) < EPSILON_LAMBDA_M
|| fabs(bf_prev - bf) < EPSILON_BF)
break;
- lambda_m_prev = lambda_m;
+ lambda_prev = lambda;
bf_prev = bf;
}
if(i >= MAX_ITER && planck->verbose) {
logger_print(planck->logger, LOG_WARNING, MSG_WARNING_PREFIX
"Could not sample a wavelength wrt the Planck distribution "
- "(spectral range = [%g, %g] nm; reference temperature = %g K)\n",
+ "(spectral range = [%g, %g] m; reference temperature = %g K)\n",
range[0], range[1], planck->ref_temperature);
}
@@ -219,14 +210,13 @@ sample_continue
const double Tref = planck->ref_temperature; /* K */
/* W/m²/sr/m */
- const double B_lambda = sbb_planck(lambda_m, lambda_m, Tref);
- const double B_mean = sbb_planck(range_m[0], range_m[1], Tref);
+ const double B_lambda = sbb_planck(lambda, lambda, Tref);
+ const double B_mean = sbb_planck(range[0], range[1], Tref);
- *pdf = B_lambda / (B_mean * (range_m[1]-range_m[0]));
- *pdf *= 1.e-9; /* Transform from m^-1 to nm^-1 */
+ *pdf = B_lambda / (B_mean * (range[1]-range[0]));
}
- return lambda_m * 1.e9; /* Convert in nanometers */
+ return lambda;
}
static FINLINE int
@@ -372,12 +362,12 @@ sbb_ran_planck_ref_put(struct sbb_ran_planck* planck)
return RES_OK;
}
-double /* Wavelength [nm] */
+double /* Wavelength [m] */
sbb_ran_planck_sample
(struct sbb_ran_planck* planck,
const double r0, /* Random number in [0, 1[ */
const double r1, /* Random number in [0, 1[ */
- double* pdf) /* [nm^-1]. May be NULL */
+ double* pdf) /* [m^-1]. May be NULL */
{
ASSERT(planck && r0 >= 0 && r1 >= 0 && r0 < 1 && r1 < 1);
