star-blackbody

Handling the Planck function
git clone git://git.meso-star.fr/star-blackbody.git
Log | Files | Refs | README | LICENSE

commit e9cd371c77d3d4ea7ce393940b2a0e0cafb55916
parent 94534dda6586253b4e08a4a7971b35c11e2281ad
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 24 Nov 2023 11:20:31 +0100

Add non regression tests

Diffstat:
MMakefile | 8+++++---
Asrc/test_sbb_misc.c | 64++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/test_sbb_ran_planck.c | 137+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 files changed, 206 insertions(+), 3 deletions(-)

diff --git a/Makefile b/Makefile @@ -113,7 +113,7 @@ lint: ################################################################################ # Tests ################################################################################ -TEST_SRC = +TEST_SRC = src/test_sbb_misc.c src/test_sbb_ran_planck.c TEST_OBJ = $(TEST_SRC:.c=.o) TEST_DEP = $(TEST_SRC:.c=.d) @@ -142,5 +142,7 @@ $(TEST_DEP): config.mk sbb-local.pc $(TEST_OBJ): config.mk sbb-local.pc $(CC) $(CFLAGS_EXE) $(SBB_CFLAGS) $(RSYS_CFLAGS) -c $(@:.o=.c) -o $@ -#TODO_EXE: config.mk sbb-local.pc $(LIBNAME) -# $(CC) $(CFLAGS_EXE) -o $@ src/$@.o $(LDFLAGS_EXE) $(SBB_LIBS) $(RSYS_LIBS) +test_sbb_misc \ +test_sbb_ran_planck \ +: config.mk sbb-local.pc $(LIBNAME) + $(CC) $(CFLAGS_EXE) -o $@ src/$@.o $(LDFLAGS_EXE) $(SBB_LIBS) $(RSYS_LIBS) -lm diff --git a/src/test_sbb_misc.c b/src/test_sbb_misc.c @@ -0,0 +1,64 @@ +/* Copyright (C) 2018-2023 |Méso|Star> (contact@meso-star.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#include "sbb.h" +#include <rsys/math.h> + +int +main(int argc, char** argv) +{ + const double lambda_lo = 3.6e-6; /* [m] */ + const double lambda_hi = 12.3e-6; /* [m] */ + const double Tref = 550; /* [K] */ + const double delta_lambda = lambda_hi - lambda_lo; + + double L; /* Radiance */ + double d; /* Temporary double */ + double T; /* Temperature */ + double ref; /* Reference value to be verified */ + (void)argc, (void)argv; + + d = sbb_blackbody_fraction(lambda_lo, lambda_hi, Tref); + ref = 0.73033200108935725; + CHK(eq_eps(d, ref, ref*1e-3)); + + d = sbb_planck_monochromatic(lambda_lo, Tref); + ref = 137.68759834775730e6; + CHK(eq_eps(d, 137.68759834775730e6, ref*1e-5)); + + d = sbb_planck_interval(lambda_lo, lambda_hi, Tref); + ref = 1206.0731201171875/delta_lambda; + CHK(eq_eps(d, ref, ref*1e-3)); + + L = 1206.0731201171875/delta_lambda; /* [W/m^2/sr/m ] */ + CHK(sbb_brightness_temperature(lambda_hi, lambda_lo, L, &T) == RES_BAD_ARG); + CHK(sbb_brightness_temperature(-lambda_lo, lambda_hi,-L, &T) == RES_BAD_ARG); + CHK(sbb_brightness_temperature(-lambda_lo,-lambda_hi,-L, &T) == RES_BAD_ARG); + CHK(sbb_brightness_temperature(lambda_lo, lambda_hi,-L, &T) == RES_BAD_ARG); + CHK(sbb_brightness_temperature(lambda_lo, lambda_hi, L, NULL) == RES_BAD_ARG); + CHK(sbb_brightness_temperature(lambda_lo, lambda_hi, L, &T) == RES_OK); + CHK(eq_eps(T, Tref, Tref*1.e-4)); + + L = 1206.0731201171875; /* [W/m^2/sr/m] */ + CHK(sbb_radiance_temperature(lambda_hi, lambda_lo, L, &T) == RES_BAD_ARG); + CHK(sbb_radiance_temperature(-lambda_lo, lambda_hi, L, &T) == RES_BAD_ARG); + CHK(sbb_radiance_temperature(-lambda_lo,-lambda_hi, L, &T) == RES_BAD_ARG); + CHK(sbb_radiance_temperature(lambda_lo, lambda_hi,-L, &T) == RES_BAD_ARG); + CHK(sbb_radiance_temperature(lambda_lo, lambda_hi, L, NULL) == RES_BAD_ARG); + CHK(sbb_radiance_temperature(lambda_lo, lambda_hi, L, &T) == RES_OK); + CHK(eq_eps(T, Tref, Tref*1e-4)); + + return 0; +} diff --git a/src/test_sbb_ran_planck.c b/src/test_sbb_ran_planck.c @@ -0,0 +1,137 @@ +/* Copyright (C) 2018-2023 |Méso|Star> (contact@meso-star.com) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#include "sbb.h" +#include <rsys/math.h> +#include <math.h> + +static double +rand_canonic(void) +{ + return (double)rand() / (double)((size_t)RAND_MAX+1); +} + +static void +check_planck(void) +{ + struct sbb_ran_planck_create_args args = SBB_RAN_PLANCK_CREATE_ARGS_DEFAULT; + struct sbb_ran_planck* planck = NULL; + + args.range[0] = 3.6e-6; + args.range[1] = 12.3e-6; + args.ref_temperature = 500; + args.nbands = 0; + + CHK(sbb_ran_planck_create(NULL, &planck) == RES_BAD_ARG); + CHK(sbb_ran_planck_create(&args, NULL) == RES_BAD_ARG); + CHK(sbb_ran_planck_create(&args, &planck) == RES_OK); + + CHK(sbb_ran_planck_ref_get(NULL) == RES_BAD_ARG); + CHK(sbb_ran_planck_ref_get(planck) == RES_OK); + CHK(sbb_ran_planck_ref_put(NULL) == RES_BAD_ARG); + CHK(sbb_ran_planck_ref_put(planck) == RES_OK); + CHK(sbb_ran_planck_ref_put(planck) == RES_OK); + + args.range[0] = 1; + args.range[1] = 0; + CHK(sbb_ran_planck_create(&args, &planck) == RES_BAD_ARG); + args.range[0] = 0; + args.range[1] = 1; + args.ref_temperature = -100; + CHK(sbb_ran_planck_create(&args, &planck) == RES_BAD_ARG); + args.range[0] = 3.6e-6; + args.range[1] =-12.3e-6; + args.ref_temperature = 300; + CHK(sbb_ran_planck_create(&args, &planck) == RES_BAD_ARG); +} + +static void +planck_integration + (const double lambda_lo, /* [m] */ + const double lambda_hi, /* [m] */ + const double Tref, /* [K] */ + const double T0, /* [K] */ + const size_t nbands) /* > 0 <=> Speed up Planck sampling */ +{ + struct sbb_ran_planck_create_args args = SBB_RAN_PLANCK_CREATE_ARGS_DEFAULT; + struct sbb_ran_planck* planck = NULL; + + const double delta_lambda = lambda_hi - lambda_lo; /* [m] */ + + /* Variables of the Monte Carlo integration */ + const size_t N = 10000; /* #realisations */ + size_t irealisation = 0; + double sum = 0; /* Sum of weights */ + double sum2 = 0; /* Sum of weights squared */ + double E = 0; /* Expected value */ + double V = 0; /* Variance */ + double SE = 0; /* Standard Error */ + + double ref = 0; + + /* Setup the planck Distribution */ + args.range[0] = lambda_lo; + args.range[1] = lambda_hi; + args.ref_temperature = Tref; + args.nbands = nbands; + CHK(sbb_ran_planck_create(&args, &planck) == RES_OK); + + FOR_EACH(irealisation, 0, N) { + const double r0 = rand_canonic(); + const double r1 = rand_canonic(); + double lambda = 0; + double pdf = 0; + double w = 0; + + lambda = sbb_ran_planck_sample(planck, r0, r1, &pdf); + w = sbb_planck_monochromatic(lambda, T0) / pdf; + sum += w; + sum2 += w*w; + } + + E = sum / (double)N; + V = MMAX(sum2 / (double)N - E*E, 0); + SE = sqrt(V/(double)N); + + ref = sbb_planck(lambda_lo, lambda_hi, T0) * delta_lambda; + + printf("planck(%g m, %g m, %g K) = %g ~ %g +/- %g [W/m^2/sr] (Tref = %g K)\n", + lambda_lo, lambda_hi, T0, ref, E, SE, Tref); + CHK(eq_eps(ref, E, 3*SE)); + + CHK(sbb_ran_planck_ref_put(planck) == RES_OK); +} + +int +main(int argc, char** argv) +{ + /* Input parameters */ + const double lambda_lo = 3.6e-6; /* [m] */ + const double lambda_hi = 12.3e-6; /* [m] */ + const double delta_lambda = lambda_hi - lambda_lo; + size_t nbands = 0; + (void)argc, (void)argv; + + check_planck(); + + nbands = (size_t)(delta_lambda * 1e9); + + planck_integration(lambda_lo, lambda_hi, 500, 550, 0); + planck_integration(lambda_lo, lambda_hi, 500, 550, nbands); + planck_integration(lambda_lo, lambda_hi, 300, 550, 0); + planck_integration(lambda_lo, lambda_hi, 300, 550, nbands); + planck_integration(lambda_lo, lambda_hi, 300, 550, nbands/4); + return 0; +}