test_sphor_MVREA_analytical1.c (8620B)
1 /* Copyright (C) 2024-2026 Centre National de la Recherche Scientifique 2 * Copyright (C) 2024-2026 Clermont Auvergne INP 3 * Copyright (C) 2024-2026 INSA Lyon 4 * Copyright (C) 2024-2026 Institut Mines Télécom Albi-Carmaux 5 * Copyright (C) 2024-2026 Institut National Polytechnique de Toulouse 6 * Copyright (C) 2024-2026 |Méso|Star> (contact@meso-star.com) 7 * Copyright (C) 2024-2026 PhotonLyX (info@photonlyx.com) 8 * Copyright (C) 2024-2026 Université de Lorraine 9 * Copyright (C) 2024-2026 Université Paul Sabatier 10 * Copyright (C) 2024-2026 Université Toulouse - Jean Jaurès 11 * 12 * This program is free software: you can redistribute it and/or modify 13 * it under the terms of the GNU General Public License as published by 14 * the Free Software Foundation, either version 3 of the License, or 15 * (at your option) any later version. 16 * 17 * This program is distributed in the hope that it will be useful, 18 * but WITHOUT ANY WARRANTY; without even the implied warranty of 19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 20 * GNU General Public License for more details. 21 * 22 * You should have received a copy of the GNU General Public License 23 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 24 #define _POSIX_C_SOURCE 200112L /* for strtok_r support */ 25 26 #include "sphor.h" 27 28 #include <rsys/cstr.h> 29 #include <rsys/logger.h> 30 #include <rsys/str.h> 31 #include <rsys/text_reader.h> 32 33 #include <stdio.h> 34 35 static void write_stl( 36 const char* filename, 37 const float verts[][3], 38 const unsigned tris[][3], 39 unsigned ntris) 40 { 41 size_t i = 0; 42 FILE* fp = fopen(filename, "w"); 43 CHK(fp != NULL); 44 45 fprintf(fp, "solid\n"); 46 FOR_EACH(i, 0, ntris) { 47 fprintf(fp, "facet normal 0 0 0\n"); 48 fprintf(fp, " outer loop\n"); 49 fprintf(fp, " vertex %f %f %f\n", SPLIT3(verts[tris[i][0]])); 50 fprintf(fp, " vertex %f %f %f\n", SPLIT3(verts[tris[i][1]])); 51 fprintf(fp, " vertex %f %f %f\n", SPLIT3(verts[tris[i][2]])); 52 fprintf(fp, " endloop\n"); 53 fprintf(fp, "endfacet\n"); 54 } 55 fprintf(fp, "endsolid\n"); 56 57 CHK(fclose(fp) == 0); 58 } 59 60 static void 61 write_cube() 62 { 63 64 const float cube_verts[][3] = { 65 {0.f, 0.f, 0.f}, 66 {1.f, 0.f, 0.f}, 67 {0.f, 1.f, 0.f}, 68 {1.f, 1.f, 0.f}, 69 {0.f, 0.f, 1.f}, 70 {1.f, 0.f, 1.f}, 71 {0.f, 1.f, 1.f}, 72 {1.f, 1.f, 1.f} 73 }; 74 75 /* Front faces are CW. Normals point outside the cube */ 76 const unsigned front_tris[][3] = { 77 {0, 2, 1}, {1, 2, 3} 78 }; 79 80 const unsigned back_tris[][3] = { 81 {4, 5, 6}, {6, 5, 7} 82 }; 83 84 const unsigned cube_tris[][3] = { 85 {0, 4, 2}, {2, 4, 6}, /* Left */ 86 {3, 7, 1}, {1, 7, 5}, /* Right */ 87 {2, 6, 3}, {3, 6, 7}, /* Top */ 88 {0, 1, 4}, {4, 1, 5} /* Bottom */ 89 }; 90 91 const unsigned front_ntris = sizeof(front_tris) / sizeof(front_tris[0]); 92 93 const unsigned back_ntris = sizeof(back_tris) / sizeof(back_tris[0]); 94 95 const unsigned walls_ntris = sizeof(cube_tris) / sizeof(cube_tris[0]); 96 97 write_stl("cube_front.stl", cube_verts, front_tris, front_ntris); 98 99 write_stl("cube_back.stl", cube_verts, back_tris, back_ntris); 100 101 write_stl("cube.stl", cube_verts, cube_tris, walls_ntris); 102 } 103 104 static void 105 write_spectral_file(const char* filename) 106 { 107 FILE* fp = NULL; 108 109 fp = fopen(filename, "w"); 110 CHK(NULL != fp); 111 112 fprintf(fp, "450 1\n"); 113 fprintf(fp, "550 1\n"); 114 CHK(fclose(fp) == 0); 115 } 116 117 static void 118 write_reflectivity_file(const char* filename) 119 { 120 FILE* fp = NULL; 121 122 fp = fopen(filename, "w"); 123 CHK(NULL != fp); 124 125 fprintf(fp, "450 0.5\n"); 126 fprintf(fp, "550 0.5\n"); 127 CHK(fclose(fp) == 0); 128 } 129 130 static void 131 write_input_file(FILE* fp) 132 { 133 const char* spectral_filename = "spec.txt"; 134 const char* reflectivity_filename = "rho.txt"; 135 136 write_cube(); 137 write_spectral_file(spectral_filename); 138 write_reflectivity_file(reflectivity_filename); 139 140 fprintf(fp, "surface: \"source\"\n"); 141 fprintf(fp, "\tgeometry: BACK cube_front.stl\n"); 142 fprintf(fp, "\tsource:\n"); 143 fprintf(fp, "\tflux_density: 1 umol/m^2/s spec.txt nm nm^-1\n"); 144 fprintf(fp, "\tdirection: COLLIM NORMAL\n"); 145 fprintf(fp, "\n"); 146 fprintf(fp, "volume: \"reactional_volume\"\n"); 147 fprintf(fp, "\tgeometry: BACK cube.stl\n"); 148 fprintf(fp, "\tgeometry: BACK cube_front.stl\n"); 149 fprintf(fp, "\tgeometry: BACK cube_back.stl\n"); 150 fprintf(fp, "\tprop_rad: \"sigma_a\" \n"); 151 fprintf(fp, "\tscatterer:\n"); 152 fprintf(fp, "\t\tconcentration: 1 mol/m^3\n"); 153 fprintf(fp, "\t\tcross_sections:\n"); 154 fprintf(fp, "\t\t\tabs_cross_sec: spec.txt nm m^2/mol\n"); 155 fprintf(fp, "\tsensor:\n"); 156 fprintf(fp, "\t\tresponse_function: 1\n"); 157 fprintf(fp, "surface: \"mirror\"\n"); 158 fprintf(fp, "\tgeometry: BACK cube_back.stl\n"); 159 fprintf(fp, "\tbrdf: SPECULAR rho.txt\n"); 160 fprintf(fp, "\tsensor:\n"); 161 fprintf(fp, "\t\tresponse_function: 1\n"); 162 CHK(fflush(fp) == 0); 163 } 164 165 int 166 main(int argc, char** argv) 167 { 168 struct sphor* sphor = NULL; 169 struct sphor_create_args args = SPHOR_CREATE_ARGS_DEFAULT; 170 struct txtrdr* txtrdr = NULL; 171 struct str line; 172 char* token = NULL; 173 char* token_ptr = NULL; 174 char input_filename[] = "input"; 175 char output_filename[] = "output"; 176 double avg = 0; 177 double std_err = 0; 178 179 FILE* fp = NULL; 180 181 fp = fopen(input_filename, "w+"); 182 CHK(NULL != fp); 183 184 write_input_file(fp); 185 CHK(fclose(fp) == 0); 186 187 args.input_filename = input_filename; 188 args.output_filename = output_filename; 189 args.force = 1; 190 191 (void)argc; 192 (void)argv; 193 194 /* Test sphor_create with NULL logger and allocator */ 195 CHK(sphor_create(&args, &sphor) == RES_OK); 196 197 /* Run sphor with the config file */ 198 CHK(sphor_run(sphor) == RES_OK); 199 200 CHK(sphor_ref_put(sphor) == RES_OK); 201 202 fp = fopen(output_filename, "r"); 203 CHK(NULL != fp); 204 205 str_init(NULL, &line); 206 CHK(txtrdr_stream(NULL, fp, output_filename, '#', &txtrdr) == RES_OK); 207 208 /* Verify first level: total scene MVREA */ 209 CHK(txtrdr_read_line(txtrdr) == RES_OK); 210 CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK); 211 token = strtok_r(str_get(&line), ":", &token_ptr); 212 CHK(strcmp(token, "MVREA") == 0); 213 token = strtok_r(NULL, ":", &token_ptr); 214 CHK(cstr_to_double(token, &avg) == RES_OK); 215 CHK(cstr_to_double(token_ptr, &std_err) == RES_OK); 216 /* Compare to analytical solution */ 217 CHK(eq_eps(avg, 1 - 0.5 * (exp(-1) + exp(-2)), 2 * std_err)); 218 219 /* Verify second level: per volume MVREA */ 220 CHK(txtrdr_read_line(txtrdr) == RES_OK); 221 CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK); 222 token = strtok_r(str_get(&line), ":", &token_ptr); 223 CHK(strcmp(token, "MVREA") == 0); 224 token = strtok_r(NULL, ":", &token_ptr); 225 CHK(strcmp(token, "reactional_volume") == 0); 226 token = strtok_r(NULL, ":", &token_ptr); 227 CHK(cstr_to_double(token, &avg) == RES_OK); 228 CHK(cstr_to_double(token_ptr, &std_err) == RES_OK); 229 /* Compare to analytical solution */ 230 CHK(eq_eps(avg, 1 - 0.5 * (exp(-1) + exp(-2)), 2 * std_err)); 231 232 /* Verify third level: per prop_rad MVREA */ 233 CHK(txtrdr_read_line(txtrdr) == RES_OK); 234 CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK); 235 token = strtok_r(str_get(&line), ":", &token_ptr); 236 CHK(strcmp(token, "MVREA") == 0); 237 token = strtok_r(NULL, ":", &token_ptr); 238 CHK(strcmp(token, "reactional_volume") == 0); 239 token = strtok_r(NULL, ":", &token_ptr); 240 CHK(strcmp(token, "sigma_a") == 0); 241 token = strtok_r(NULL, ":", &token_ptr); 242 CHK(cstr_to_double(token, &avg) == RES_OK); 243 CHK(cstr_to_double(token_ptr, &std_err) == RES_OK); 244 /* Compare to analytical solution */ 245 CHK(eq_eps(avg, 1 - 0.5 * (exp(-1) + exp(-2)), 2 * std_err)); 246 247 /* Verify first level: total scene losses */ 248 CHK(txtrdr_read_line(txtrdr) == RES_OK); 249 CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK); 250 token = strtok_r(str_get(&line), ":", &token_ptr); 251 CHK(strcmp(token, "LOSSES") == 0); 252 token = strtok_r(NULL, ":", &token_ptr); 253 CHK(cstr_to_double(token, &avg) == RES_OK); 254 CHK(cstr_to_double(token_ptr, &std_err) == RES_OK); 255 /* Compare to analytical solution */ 256 CHK(eq_eps(avg, 0.5 * exp(-1), 2 * std_err)); 257 258 /* Verify second level: per surface losses */ 259 CHK(txtrdr_read_line(txtrdr) == RES_OK); 260 CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK); 261 token = strtok_r(str_get(&line), ":", &token_ptr); 262 CHK(strcmp(token, "LOSSES") == 0); 263 token = strtok_r(NULL, ":", &token_ptr); 264 CHK(strcmp(token, "mirror") == 0); 265 token = strtok_r(NULL, ":", &token_ptr); 266 CHK(cstr_to_double(token, &avg) == RES_OK); 267 CHK(cstr_to_double(token_ptr, &std_err) == RES_OK); 268 /* Compare to analytical solution */ 269 CHK(eq_eps(avg, 0.5 * exp(-1), 2 * std_err)); 270 271 CHK(fclose(fp) == 0); 272 273 str_release(&line); 274 txtrdr_ref_put(txtrdr); 275 276 CHK(mem_allocated_size() == 0); 277 return 0; 278 }