commit 634507724af9edfb9ed5cc5ca08eaf3c731ae54a
parent 34cfdb97f8eb4e6ee8245ec04f67197cb29f4051
Author: Eduardo Fontana Lazzari <edufonlaz@gmail.com>
Date: Thu, 29 Jan 2026 10:29:48 +0100
Add analytical test case to the MVREA algorithm
Previous tests checked only the behavior of the library, without
actually running a simulation. This commit adds a new test using a
simple configuration and compares it to the analytical solution. It
verifies whether the analytical result lies within the 95% confidence
interval of the simulated value and also verifies the conformity of the
output file with the star-phor-output file format specification.
Diffstat:
3 files changed, 275 insertions(+), 4 deletions(-)
diff --git a/.gitignore b/.gitignore
@@ -10,7 +10,8 @@ spec.txt
*.stl
tags
input
+output
star-phor.1
star-phor
-test_sphor
-test_sphor_lib
+test_*
+!test*.[ch]
diff --git a/Makefile b/Makefile
@@ -168,7 +168,8 @@ uninstall:
# Tests
################################################################################
TEST_SRC =\
- src/test_sphor_lib.c
+ src/test_sphor_lib.c\
+ src/test_sphor_MVREA_analytical1.c
TEST_OBJ = $(TEST_SRC:.c=.o)
TEST_DEP = $(TEST_SRC:.c=.d)
TEST_TGT = $(TEST_SRC:.c=.t)
@@ -215,9 +216,11 @@ $(TEST_OBJ): config.mk sphor-local.pc
$(CC) $(CFLAGS_TEST) -c $(@:.o=.c) -o $@
test_sphor_lib\
+test_sphor_MVREA_analytical1\
: config.mk sphor-local.pc $(LIBNAME)
$(CC) $(CFLAGS_TEST) -o $@ src/$@.o $(LDFLAGS_TEST)
clean_test:
- rm -f $(TEST_DEP) $(TEST_OBJ) $(TEST_TGT) input cube.stl spec.txt
+ rm -f $(TEST_DEP) $(TEST_OBJ) $(TEST_TGT) \
+ input output cube.stl cube_back.stl cube_front.stl spec.txt
for i in $(TEST_SRC); do rm -f "$$(basename "$${i}" ".c")"; done
diff --git a/src/test_sphor_MVREA_analytical1.c b/src/test_sphor_MVREA_analytical1.c
@@ -0,0 +1,267 @@
+/* Copyright (C) 2024-2025 Centre National de la Recherche Scientifique
+ * Copyright (C) 2024-2025 Clermont Auvergne INP
+ * Copyright (C) 2024-2025 INSA Lyon
+ * Copyright (C) 2024-2025 Institut Mines Télécom Albi-Carmaux
+ * Copyright (C) 2024-2025 Institut National Polytechnique de Toulouse
+ * Copyright (C) 2024-2025 |Méso|Star> (contact@meso-star.com)
+ * Copyright (C) 2024-2025 PhotonLyX (info@photonlyx.com)
+ * Copyright (C) 2024-2025 Université de Lorraine
+ * Copyright (C) 2024-2025 Université Paul Sabatier
+ * Copyright (C) 2024-2025 Université Toulouse - Jean Jaurès
+ *
+ * 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/>. */
+#define _POSIX_C_SOURCE 200112L /* for strtok_r support */
+
+#include "sphor.h"
+
+#include <rsys/cstr.h>
+#include <rsys/logger.h>
+#include <rsys/str.h>
+#include <rsys/text_reader.h>
+
+#include <math.h>
+#include <stdio.h>
+
+static void write_stl(
+ const char* filename,
+ const float verts[][3],
+ const unsigned tris[][3],
+ unsigned ntris)
+{
+ size_t i = 0;
+ FILE* fp = fopen(filename, "w");
+ CHK(fp != NULL);
+
+ fprintf(fp, "solid\n");
+ FOR_EACH(i, 0, ntris) {
+ fprintf(fp, "facet normal 0 0 0\n");
+ fprintf(fp, " outer loop\n");
+ fprintf(fp, " vertex %f %f %f\n", SPLIT3(verts[tris[i][0]]));
+ fprintf(fp, " vertex %f %f %f\n", SPLIT3(verts[tris[i][1]]));
+ fprintf(fp, " vertex %f %f %f\n", SPLIT3(verts[tris[i][2]]));
+ fprintf(fp, " endloop\n");
+ fprintf(fp, "endfacet\n");
+ }
+ fprintf(fp, "endsolid\n");
+
+ CHK(fclose(fp) == 0);
+}
+
+static void
+write_cube()
+{
+
+ const float cube_verts[][3] = {
+ {0.f, 0.f, 0.f},
+ {1.f, 0.f, 0.f},
+ {0.f, 1.f, 0.f},
+ {1.f, 1.f, 0.f},
+ {0.f, 0.f, 1.f},
+ {1.f, 0.f, 1.f},
+ {0.f, 1.f, 1.f},
+ {1.f, 1.f, 1.f}
+ };
+
+ /* Front faces are CW. Normals point outside the cube */
+ const unsigned front_tris[][3] = {
+ {0, 2, 1}, {1, 2, 3}
+ };
+
+ const unsigned back_tris[][3] = {
+ {4, 5, 6}, {6, 5, 7}
+ };
+
+ const unsigned cube_tris[][3] = {
+ {0, 4, 2}, {2, 4, 6}, /* Left */
+ {3, 7, 1}, {1, 7, 5}, /* Right */
+ {2, 6, 3}, {3, 6, 7}, /* Top */
+ {0, 1, 4}, {4, 1, 5} /* Bottom */
+ };
+
+ const unsigned front_ntris = sizeof(front_tris) / sizeof(front_tris[0]);
+
+ const unsigned back_ntris = sizeof(back_tris) / sizeof(back_tris[0]);
+
+ const unsigned walls_ntris = sizeof(cube_tris) / sizeof(cube_tris[0]);
+
+ write_stl("cube_front.stl", cube_verts, front_tris, front_ntris);
+
+ write_stl("cube_back.stl", cube_verts, back_tris, back_ntris);
+
+ write_stl("cube.stl", cube_verts, cube_tris, walls_ntris);
+}
+
+static void
+write_spectral_file(const char* filename)
+{
+ FILE* fp = NULL;
+
+ fp = fopen(filename, "w");
+ CHK(NULL != fp);
+
+ fprintf(fp, "450 1\n");
+ fprintf(fp, "550 1\n");
+ CHK(fclose(fp) == 0);
+}
+
+
+static void
+write_input_file(FILE* fp)
+{
+ const char* spectral_filename = "spec.txt";
+
+ write_cube();
+ write_spectral_file(spectral_filename);
+
+ fprintf(fp, "surface: \"source\"\n");
+ fprintf(fp, "\tgeometry: BACK cube_front.stl\n");
+ fprintf(fp, "\tsource:\n");
+ fprintf(fp, "\tflux_density: 1 umol/m^2/s spec.txt nm nm^-1\n");
+ fprintf(fp, "\tdirection: COLLIM NORMAL\n");
+ fprintf(fp, "\n");
+ fprintf(fp, "volume: \"reactional_volume\"\n");
+ fprintf(fp, "\tgeometry: BACK cube.stl\n");
+ fprintf(fp, "\tgeometry: BACK cube_front.stl\n");
+ fprintf(fp, "\tgeometry: BACK cube_back.stl\n");
+ fprintf(fp, "\tprop_rad: \"sigma_a\" SCATTERER\n");
+ fprintf(fp, "\t\tconcentration: 1 mol/m^3\n");
+ fprintf(fp, "\t\tcross_sections:\n");
+ fprintf(fp, "\t\t\tabs_cross_sec: spec.txt nm m^2/mol\n");
+ fprintf(fp, "\tsensor:\n");
+ fprintf(fp, "\t\tresponse_function: 1\n");
+ fprintf(fp, "surface: \"mirror\"\n");
+ fprintf(fp, "\tgeometry: BACK cube_back.stl\n");
+ fprintf(fp, "\tbrdf: SPECULAR 0.5\n");
+ fprintf(fp, "\tsensor:\n");
+ fprintf(fp, "\t\tresponse_function: 1\n");
+ CHK(fflush(fp) == 0);
+}
+
+int
+main(int argc, char** argv)
+{
+ struct sphor* sphor = NULL;
+ struct sphor_create_args args = SPHOR_CREATE_ARGS_DEFAULT;
+ struct txtrdr* txtrdr = NULL;
+ struct str line;
+ char* token = NULL;
+ char* token_ptr = NULL;
+ char input_filename[] = "input";
+ char output_filename[] = "output";
+ double avg = 0;
+ double std = 0;
+ double nsamples = 0;
+
+ FILE* fp = NULL;
+
+ fp = fopen(input_filename, "w+");
+ CHK(NULL != fp);
+
+ write_input_file(fp);
+ CHK(fclose(fp) == 0);
+
+ args.input_filename = input_filename;
+ args.output_filename = output_filename;
+ args.force = 1;
+
+ (void)argc;
+ (void)argv;
+
+ /* Test sphor_create with NULL logger and allocator */
+ CHK(sphor_create(&args, &sphor) == RES_OK);
+
+ /* Run sphor with the config file */
+ CHK(sphor_run(sphor) == RES_OK);
+
+ CHK(sphor_ref_put(sphor) == RES_OK);
+
+ fp = fopen(output_filename, "r");
+ CHK(NULL != fp);
+
+ str_init(NULL, &line);
+ CHK(txtrdr_stream(NULL, fp, output_filename, '#', &txtrdr) == RES_OK);
+
+ nsamples = (double)args.samples;
+
+ /* Verify first level: total scene MVREA */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "MVREA") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare to analytical solution */
+ CHK(eq_eps(avg, 1 - 0.5 * (exp(-1) + exp(-2)), 2 * std / sqrt(nsamples)));
+
+ /* Verify second level: per volume MVREA */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "MVREA") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(strcmp(token, "reactional_volume") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare to analytical solution */
+ CHK(eq_eps(avg, 1 - 0.5 * (exp(-1) + exp(-2)), 2 * std / sqrt(nsamples)));
+
+ /* Verify third level: per prop_rad MVREA */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "MVREA") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(strcmp(token, "reactional_volume") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(strcmp(token, "sigma_a") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare to analytical solution */
+ CHK(eq_eps(avg, 1 - 0.5 * (exp(-1) + exp(-2)), 2 * std / sqrt(nsamples)));
+
+ /* Verify first level: total scene losses */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "LOSSES") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare to analytical solution */
+ CHK(eq_eps(avg, 0.5 * exp(-1), 2 * std / sqrt(nsamples)));
+
+ /* Verify second level: per surface losses */
+ CHK(txtrdr_read_line(txtrdr) == RES_OK);
+ CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK);
+ token = strtok_r(str_get(&line), ":", &token_ptr);
+ CHK(strcmp(token, "LOSSES") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(strcmp(token, "mirror") == 0);
+ token = strtok_r(NULL, ":", &token_ptr);
+ CHK(cstr_to_double(token, &avg) == RES_OK);
+ CHK(cstr_to_double(token_ptr, &std) == RES_OK);
+ /* Compare to analytical solution */
+ CHK(eq_eps(avg, 0.5 * exp(-1), 2 * std / sqrt(nsamples)));
+
+ CHK(fclose(fp) == 0);
+
+ str_release(&line);
+ txtrdr_ref_put(txtrdr);
+
+ CHK(mem_allocated_size() == 0);
+ return 0;
+}
