test_atrri_load.c (10610B)
1 /* Copyright (C) 2022, 2023 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2021 Centre National de la Recherche Scientifique 3 * 4 * This program is free software: you can redistribute it and/or modify 5 * it under the terms of the GNU General Public License as published by 6 * the Free Software Foundation, either version 3 of the License, or 7 * (at your option) any later version. 8 * 9 * This program is distributed in the hope that it will be useful, 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12 * GNU General Public License for more details. 13 * 14 * You should have received a copy of the GNU General Public License 15 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 16 17 #define _POSIX_C_SOURCE 200112L /* nextafter support */ 18 19 #include "atrri.h" 20 21 #include <rsys/math.h> 22 #include <rsys/mem_allocator.h> 23 24 #include <float.h> 25 #include <math.h> 26 27 /******************************************************************************* 28 * Helper functions 29 ******************************************************************************/ 30 static INLINE double 31 rand_canonic() 32 { 33 return (double)rand() / (double)((size_t)RAND_MAX + 1); 34 } 35 36 static void 37 check_desc 38 (const struct atrri_desc* desc, 39 const double refract_ids[], 40 const size_t nrefract_ids) 41 { 42 size_t i; 43 44 CHK(desc && refract_ids && nrefract_ids); 45 CHK(desc->nindices == nrefract_ids); 46 FOR_EACH(i, 0, nrefract_ids) { 47 if(i != 0) { 48 CHK(desc->indices[i].wavelength > desc->indices[i-1].wavelength); 49 } 50 CHK(desc->indices[i].wavelength == refract_ids[i*3+0]); 51 CHK(desc->indices[i].n == refract_ids[i*3+1]); 52 CHK(desc->indices[i].kappa == refract_ids[i*3+2]); 53 } 54 } 55 56 static void 57 fetch_refractive_index 58 (const struct atrri_desc* desc, 59 const double wavelength, 60 struct atrri_refractive_index* index) 61 { 62 size_t i; 63 CHK(desc && wavelength >= 0 && index); 64 65 /* Naive linear search the first refractive index whose associated wavelength 66 * is greater than or equal to the submitted wavelength */ 67 FOR_EACH(i, 0, desc->nindices) { 68 if(desc->indices[i].wavelength >= wavelength) break; 69 } 70 71 if(i >= desc->nindices) { 72 *index = desc->indices[desc->nindices-1]; 73 } else if(i == 0) { 74 *index = desc->indices[0]; 75 } else { 76 const struct atrri_refractive_index* next = desc->indices + i; 77 const struct atrri_refractive_index* prev = desc->indices + i - 1; 78 const double len = next->wavelength - prev->wavelength; 79 const double u = CLAMP((wavelength - prev->wavelength) / len, 0, 1); 80 81 CHK(wavelength > prev->wavelength); 82 CHK(wavelength <= next->wavelength); 83 84 index->wavelength = wavelength; 85 index->n = u * prev->n + (1.0-u) * next->n; 86 index->kappa = u * prev->kappa + (1.0-u) * next->kappa; 87 } 88 } 89 90 static void 91 test_load(struct atrri* atrri) 92 { 93 const double refract_ids[] = { /* Wavelength n kappa */ 94 300.000011921, 1.740000010, 0.469999999, 95 310.000002384, 1.741990328, 0.469004273, 96 319.999992847, 1.743917465, 0.468010038, 97 329.999983311, 1.745785475, 0.467014253, 98 339.999973774, 1.748130083, 0.465993971, 99 349.999964237, 1.750000000, 0.464985400, 100 359.999954700, 1.750000000, 0.463959038, 101 369.999945164, 1.750000000, 0.462961257, 102 379.999935627, 1.750000000, 0.461990029, 103 389.999926090, 1.750000000, 0.460981935, 104 399.999916553, 1.750000000, 0.459957659, 105 409.999907017, 1.750000000, 0.458958119, 106 419.999897480, 1.750000000, 0.457983255, 107 429.999887943, 1.750000000, 0.456973791, 108 439.999878407, 1.750000000, 0.455950528, 109 449.999868870, 1.750000000, 0.454950303, 110 459.999859333, 1.750000000, 0.453971595, 111 469.999849796, 1.750000000, 0.452998221, 112 479.999840260, 1.750000000, 0.451974779, 113 489.999830723, 1.750000000, 0.450974643, 114 499.999821186, 1.750000000, 0.449385107, 115 509.999811649, 1.750000000, 0.447702825, 116 519.999802113, 1.750000000, 0.445923656, 117 529.999792576, 1.750000000, 0.443915308, 118 539.999783039, 1.750000000, 0.441953331, 119 549.999773502, 1.750000000, 0.440498799, 120 559.999763966, 1.750000000, 0.439258754, 121 569.999754429, 1.750000000, 0.438044012, 122 579.999744892, 1.750000000, 0.437014997, 123 589.999735355, 1.750000000, 0.436012030, 124 599.999725819, 1.750000000, 0.435097486, 125 609.999716282, 1.750000000, 0.435064077, 126 619.999706745, 1.750000000, 0.435031146, 127 629.999697208, 1.750000000, 0.434985548, 128 639.999687672, 1.750000000, 0.434602797, 129 649.999678135, 1.750000000, 0.434226334, 130 659.999668598, 1.750000000, 0.433855951, 131 669.999659061, 1.750000000, 0.432947546, 132 679.999649525, 1.750000000, 0.431957632, 133 689.999639988, 1.750000000, 0.430984408, 134 699.999630451, 1.750000000, 0.430428028, 135 709.999620914, 1.750000000, 0.430284500, 136 719.999611378, 1.750000000, 0.430143028, 137 729.999601841, 1.750000000, 0.430003524, 138 739.999592304, 1.750000000, 0.430000007, 139 749.999582767, 1.750000000, 0.430000007, 140 759.999573231, 1.750000000, 0.430000007, 141 769.999563694, 1.750000000, 0.430021703, 142 779.999554157, 1.750000000, 0.430100024, 143 789.999544621, 1.750000000, 0.430177301, 144 799.999535084, 1.750000000, 0.430253685 145 }; 146 const size_t nrefract_ids = sizeof(refract_ids) / (sizeof(double)*3); 147 struct atrri_desc desc = ATRRI_DESC_NULL; 148 struct atrri_refractive_index index = ATRRI_REFRACTIVE_INDEX_NULL; 149 size_t irefract_id = 0; 150 const char* filename = "test_file.atrri"; 151 FILE* fp = NULL; 152 double wlen; 153 size_t i; 154 155 fp = fopen(filename, "w+"); 156 CHK(fp); 157 158 /* Write the input file */ 159 fprintf(fp, "# Wavelength (in nm) n kappa\n"); 160 FOR_EACH(irefract_id, 0, nrefract_ids) { 161 fprintf(fp, "%.9f %.9f %.9f\n", 162 refract_ids[irefract_id*3+0], 163 refract_ids[irefract_id*3+1], 164 refract_ids[irefract_id*3+2]); 165 } 166 167 rewind(fp); 168 CHK(atrri_load_stream(NULL, fp, filename) == RES_BAD_ARG); 169 CHK(atrri_load_stream(atrri, NULL, filename) == RES_BAD_ARG); 170 CHK(atrri_load_stream(atrri, fp, NULL) == RES_OK); 171 CHK(atrri_get_desc(NULL, &desc) == RES_BAD_ARG); 172 CHK(atrri_get_desc(atrri, NULL) == RES_BAD_ARG); 173 CHK(atrri_get_desc(atrri, &desc) == RES_OK); 174 check_desc(&desc, refract_ids, nrefract_ids); 175 176 rewind(fp); 177 CHK(atrri_load_stream(atrri, fp, filename) == RES_OK); 178 check_desc(&desc, refract_ids, nrefract_ids); 179 180 wlen = desc.indices[0].wavelength; 181 CHK(atrri_fetch_refractive_index(NULL, wlen, &index) == RES_BAD_ARG); 182 CHK(atrri_fetch_refractive_index(atrri, -1, &index) == RES_BAD_ARG); 183 CHK(atrri_fetch_refractive_index(atrri, wlen, NULL) == RES_BAD_ARG); 184 CHK(atrri_fetch_refractive_index(atrri, wlen, &index) == RES_OK); 185 CHK(index.wavelength == desc.indices[0].wavelength); 186 CHK(index.n == desc.indices[0].n); 187 CHK(index.kappa == desc.indices[0].kappa); 188 189 /* Test clamp to lower bound */ 190 i = 0; 191 wlen = nextafter(wlen, 0); 192 CHK(atrri_fetch_refractive_index(atrri, wlen, &index) == RES_OK); 193 CHK(index.wavelength == desc.indices[i].wavelength); 194 CHK(index.n == desc.indices[i].n); 195 CHK(index.kappa == desc.indices[i].kappa); 196 197 /* Test clamp to upper bound */ 198 i = desc.nindices-1; 199 wlen = nextafter(desc.indices[i].wavelength, DBL_MAX); 200 CHK(atrri_fetch_refractive_index(atrri, wlen, &index) == RES_OK); 201 CHK(index.wavelength == desc.indices[i].wavelength); 202 CHK(index.n == desc.indices[i].n); 203 CHK(index.kappa == desc.indices[i].kappa); 204 205 /* Test random fetches */ 206 FOR_EACH(i, 0, 10000) { 207 struct atrri_refractive_index index_ref = ATRRI_REFRACTIVE_INDEX_NULL; 208 const double r = rand_canonic(); 209 wlen = 210 r * (desc.indices[desc.nindices-1].wavelength - desc.indices[0].wavelength) 211 + desc.indices[0].wavelength; 212 213 CHK(atrri_fetch_refractive_index(atrri, wlen, &index) == RES_OK); 214 215 fetch_refractive_index(&desc, wlen, &index_ref); 216 CHK(index_ref.wavelength == index.wavelength); 217 CHK(eq_eps(index_ref.n, index.n, 1.e-6)); 218 CHK(eq_eps(index_ref.kappa, index.kappa, 1.e-6)); 219 } 220 221 fclose(fp); 222 } 223 224 static void 225 test_load_failure(struct atrri* atrri) 226 { 227 FILE* fp = NULL; 228 229 fp = tmpfile(); 230 CHK(fp); 231 fprintf(fp, "# Empty file\n"); 232 rewind(fp); 233 CHK(atrri_load_stream(atrri, fp, NULL) == RES_BAD_ARG); 234 fclose(fp); 235 236 fp = tmpfile(); 237 fprintf(fp, "-1 10 10\n"); 238 rewind(fp); 239 CHK(atrri_load_stream(atrri, fp, NULL) == RES_BAD_ARG); 240 fclose(fp); 241 242 fp = tmpfile(); 243 fprintf(fp, "1 -1 1\n"); 244 rewind(fp); 245 CHK(atrri_load_stream(atrri, fp, NULL) == RES_BAD_ARG); 246 fclose(fp); 247 248 fp = tmpfile(); 249 fprintf(fp, "1 1 -1\n"); 250 rewind(fp); 251 CHK(atrri_load_stream(atrri, fp, NULL) == RES_BAD_ARG); 252 fclose(fp); 253 254 fp = tmpfile(); 255 fprintf(fp, "1 1 1\n"); 256 fprintf(fp, "3 3 3\n"); 257 fprintf(fp, "2 2 2\n"); 258 rewind(fp); 259 CHK(atrri_load_stream(atrri, fp, NULL) == RES_BAD_ARG); 260 fclose(fp); 261 262 fp = tmpfile(); 263 fprintf(fp, "1, 1 1\n"); 264 rewind(fp); 265 CHK(atrri_load_stream(atrri, fp, NULL) == RES_BAD_ARG); 266 fclose(fp); 267 268 fp = tmpfile(); 269 fprintf(fp, "1 1; 1\n"); 270 rewind(fp); 271 CHK(atrri_load_stream(atrri, fp, NULL) == RES_BAD_ARG); 272 fclose(fp); 273 274 fp = tmpfile(); 275 fprintf(fp, "1 1 1.0.1\n"); 276 rewind(fp); 277 CHK(atrri_load_stream(atrri, fp, NULL) == RES_BAD_ARG); 278 fclose(fp); 279 280 fp = tmpfile(); 281 fprintf(fp, "1 1 1 hello world !\n"); 282 rewind(fp); 283 CHK(atrri_load_stream(atrri, fp, NULL) == RES_OK); /* Warning */ 284 fclose(fp); 285 } 286 287 static void 288 test_load_files(struct atrri* atrri, int argc, char** argv) 289 { 290 int iarg; 291 CHK(atrri); 292 FOR_EACH(iarg, 1, argc) { 293 struct atrri_desc desc = ATRRI_DESC_NULL; 294 size_t i; 295 296 printf("Load %s\n", argv[iarg]); 297 CHK(atrri_load(atrri, argv[iarg]) == RES_OK); 298 CHK(atrri_get_desc(atrri, &desc) == RES_OK); 299 CHK(desc.nindices != 0); 300 301 FOR_EACH(i, 0, desc.nindices) { 302 CHK(desc.indices[i].wavelength == desc.indices[i].wavelength); /* !NaN */ 303 CHK(desc.indices[i].n == desc.indices[i].n); /* !NaN */ 304 CHK(desc.indices[i].kappa == desc.indices[i].kappa); /* !NaN */ 305 CHK(desc.indices[i].wavelength >= 0); 306 CHK(desc.indices[i].n >= 0); 307 CHK(desc.indices[i].kappa >= 0); 308 } 309 } 310 } 311 312 313 /******************************************************************************* 314 * Main function 315 ******************************************************************************/ 316 int 317 main(int argc, char** argv) 318 { 319 struct atrri* atrri = NULL; 320 (void)argc, (void)argv; 321 322 CHK(atrri_create(NULL, &mem_default_allocator, 1, &atrri) == RES_OK); 323 324 test_load(atrri); 325 test_load_failure(atrri); 326 test_load_files(atrri, argc, argv); 327 328 CHK(atrri_ref_put(atrri) == RES_OK); 329 CHK(mem_allocated_size() == 0); 330 return 0; 331 }