test_sphor_MVREA_btdf_lambertian.c (8056B)
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 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 36 write_emitting_surface() 37 { 38 FILE* fp = fopen("emitting_surface.stl", "w"); 39 CHK(fp != NULL); 40 41 fprintf(fp, 42 "solid emitting_surface\n" 43 " facet normal 0.5 0.8660254037844386 0\n" 44 " outer loop\n" 45 " vertex -1.4433756729740643 -1.5 -1\n" 46 " vertex -0.5773502691896257 -2 1\n" 47 " vertex -0.5773502691896257 -2 -1\n" 48 " endloop\n" 49 " endfacet\n" 50 " facet normal 0.5 0.8660254037844386 0\n" 51 " outer loop\n" 52 " vertex -0.5773502691896257 -2 1\n" 53 " vertex -1.4433756729740643 -1.5 -1\n" 54 " vertex -1.4433756729740643 -1.5 1\n" 55 " endloop\n" 56 " endfacet\n" 57 "endsolid emitting_surface\n"); 58 59 CHK(fclose(fp) == 0); 60 } 61 62 static void 63 write_btdf_surface() 64 { 65 FILE* fp = fopen("btdf_surface.stl", "w"); 66 CHK(fp != NULL); 67 68 fprintf(fp, 69 "solid btdf_surface\n" 70 " facet normal 1 0 0\n" 71 " outer loop\n" 72 " vertex 0 -1 -1\n" 73 " vertex 0 1 1\n" 74 " vertex 0 -1 1\n" 75 " endloop\n" 76 " endfacet\n" 77 " facet normal 1 0 0\n" 78 " outer loop\n" 79 " vertex 0 -1 -1\n" 80 " vertex 0 1 -1\n" 81 " vertex 0 1 1\n" 82 " endloop\n" 83 " endfacet\n" 84 "endsolid btdf_surface\n"); 85 86 CHK(fclose(fp) == 0); 87 } 88 89 static void 90 write_absorbing_surface() 91 { 92 FILE* fp = fopen("absorbing_surface.stl", "w"); 93 CHK(fp != NULL); 94 95 fprintf(fp, 96 "solid absorbing_surface\n" 97 " facet normal 1 0 0\n" 98 " outer loop\n" 99 " vertex 1 -1 -1\n" 100 " vertex 1 1 1\n" 101 " vertex 1 -1 1\n" 102 " endloop\n" 103 " endfacet\n" 104 " facet normal 1 0 0\n" 105 " outer loop\n" 106 " vertex 1 -1 -1\n" 107 " vertex 1 1 -1\n" 108 " vertex 1 1 1\n" 109 " endloop\n" 110 " endfacet\n" 111 "endsolid absorbing_surface\n"); 112 113 CHK(fclose(fp) == 0); 114 } 115 116 static void 117 write_spectrum() 118 { 119 FILE* fp = fopen("emission_spectrum.dat", "w"); 120 CHK(fp != NULL); 121 122 fprintf(fp, "400 1\n"); 123 fprintf(fp, "401 1\n"); 124 125 CHK(fclose(fp) == 0); 126 } 127 128 static void 129 write_btdf() 130 { 131 FILE* fp = fopen("btdf.dat", "w"); 132 CHK(fp != NULL); 133 134 fprintf(fp, "400 1\n"); 135 fprintf(fp, "401 1\n"); 136 137 CHK(fclose(fp) == 0); 138 } 139 140 static void 141 write_bsdf_null() 142 { 143 FILE* fp = fopen("bsdf_null.dat", "w"); 144 CHK(fp != NULL); 145 146 fprintf(fp, "400 0\n"); 147 fprintf(fp, "401 0\n"); 148 149 CHK(fclose(fp) == 0); 150 } 151 152 static void 153 write_input_file(FILE* fp) 154 { 155 write_emitting_surface(); 156 write_btdf_surface(); 157 write_absorbing_surface(); 158 write_spectrum(); 159 write_btdf(); 160 write_bsdf_null(); 161 162 fprintf(fp, "surface: \"source\"\n"); 163 fprintf(fp, "\tgeometry: FRONT emitting_surface.stl\n"); 164 fprintf(fp, "\tsource:\n"); 165 fprintf(fp, "\tflux_density: 2 umol/m^2/s emission_spectrum.dat nm nm^-1\n"); 166 fprintf(fp, "\t\tdirection: COLLIM NORMAL\n"); 167 fprintf(fp, "\n"); 168 169 fprintf(fp, "surface: \"transmitter\"\n"); 170 fprintf(fp, "\tgeometry: BACK btdf_surface.stl\n"); 171 fprintf(fp, "\tbtdf: LAMBERT btdf.dat\n"); 172 fprintf(fp, "\n"); 173 174 fprintf(fp, "surface: \"absorber\"\n"); 175 fprintf(fp, "\tgeometry: BACK absorbing_surface.stl\n"); 176 fprintf(fp, "\tbrdf: SPECULAR bsdf_null.dat\n"); 177 fprintf(fp, "\tsensor:\n"); 178 fprintf(fp, "\t\tresponse_function: 1\n"); 179 180 CHK(fflush(fp) == 0); 181 } 182 183 int main(int argc, char** argv) 184 { 185 struct sphor* sphor = NULL; 186 struct sphor_create_args args = SPHOR_CREATE_ARGS_DEFAULT; 187 struct txtrdr* txtrdr = NULL; 188 struct str line; 189 char* token = NULL; 190 char* token_ptr = NULL; 191 char input_filename[] = "input"; 192 char output_filename[] = "output"; 193 double avg = 0; 194 double std_err = 0; 195 196 FILE* fp = NULL; 197 198 fp = fopen(input_filename, "w+"); 199 CHK(NULL != fp); 200 201 write_input_file(fp); 202 CHK(fclose(fp) == 0); 203 204 args.input_filename = input_filename; 205 args.output_filename = output_filename; 206 args.force = 1; 207 208 (void)argc; 209 (void)argv; 210 211 /* Test sphor_create with NULL logger and allocator */ 212 CHK(sphor_create(&args, &sphor) == RES_OK); 213 214 /* Run sphor with the config file */ 215 CHK(sphor_run(sphor) == RES_OK); 216 217 CHK(sphor_ref_put(sphor) == RES_OK); 218 219 fp = fopen(output_filename, "r"); 220 CHK(NULL != fp); 221 222 str_init(NULL, &line); 223 CHK(txtrdr_stream(NULL, fp, output_filename, '#', &txtrdr) == RES_OK); 224 225 /* Verify first level: total scene MVREA */ 226 CHK(txtrdr_read_line(txtrdr) == RES_OK); 227 CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK); 228 token = strtok_r(str_get(&line), ":", &token_ptr); 229 CHK(strcmp(token, "MVREA") == 0); 230 token = strtok_r(NULL, ":", &token_ptr); 231 CHK(cstr_to_double(token, &avg) == RES_OK); 232 CHK(cstr_to_double(token_ptr, &std_err) == RES_OK); 233 /* Compare to analytical solution */ 234 CHK(eq_eps(avg, 0, 2 * std_err)); 235 236 /* Verify first level: total scene losses */ 237 CHK(txtrdr_read_line(txtrdr) == RES_OK); 238 CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK); 239 token = strtok_r(str_get(&line), ":", &token_ptr); 240 CHK(strcmp(token, "LOSSES") == 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 against the analytical view factor. 245 * 246 * This computes the exact configuration (view) factor F12 between two 247 * coaxial, parallel square surfaces of side length 2 separated by a 248 * distance l = 1. 249 * 250 * The expression evaluated is: 251 * 252 * F12 = (1 / A1) int_{A1} int_{A2} ( l^2 / (pi r^4) ) dA2 dA1 253 * 254 * where: 255 * - A1 and A2 are the two square surfaces, 256 * - l is the axial separation, 257 * - r = ||p2 - p1|| is the distance between differential elements, 258 * - the integrand results from the general view factor formula 259 * (cos theta_1 cos theta_2) / (pi r^2), 260 * using cos theta_1 = cos theta_2 = l / r for parallel planes. 261 * 262 * The result is the density of diffuse power leaving the btdf surface 263 * that reaches the absorbing surface, weighted by the surface area ratio of 264 * the emitting surface and the absorbing surface. 265 */ 266 CHK(eq_eps(avg, 0.4152532835771472, 3 * std_err)); 267 268 /* Verify second level: per surface losses */ 269 CHK(txtrdr_read_line(txtrdr) == RES_OK); 270 CHK(str_set(&line, txtrdr_get_cline(txtrdr)) == RES_OK); 271 token = strtok_r(str_get(&line), ":", &token_ptr); 272 CHK(strcmp(token, "LOSSES") == 0); 273 token = strtok_r(NULL, ":", &token_ptr); 274 CHK(strcmp(token, "absorber") == 0); 275 token = strtok_r(NULL, ":", &token_ptr); 276 CHK(cstr_to_double(token, &avg) == RES_OK); 277 CHK(cstr_to_double(token_ptr, &std_err) == RES_OK); 278 CHK(eq_eps(avg, 0.4152532835771472, 3 * std_err)); 279 280 CHK(fclose(fp) == 0); 281 282 str_release(&line); 283 txtrdr_ref_put(txtrdr); 284 285 CHK(mem_allocated_size() == 0); 286 287 return 0; 288 }