test_htrdr_planets_source.c (11691B)
1 /* Copyright (C) 2018-2019, 2022-2025 Centre National de la Recherche Scientifique 2 * Copyright (C) 2020-2022 Institut Mines Télécom Albi-Carmaux 3 * Copyright (C) 2022-2025 Institut Pierre-Simon Laplace 4 * Copyright (C) 2022-2025 Institut de Physique du Globe de Paris 5 * Copyright (C) 2018-2025 |Méso|Star> (contact@meso-star.com) 6 * Copyright (C) 2022-2025 Observatoire de Paris 7 * Copyright (C) 2022-2025 Université de Reims Champagne-Ardenne 8 * Copyright (C) 2022-2025 Université de Versaille Saint-Quentin 9 * Copyright (C) 2018-2019, 2022-2025 Université Paul Sabatier 10 * 11 * This program is free software: you can redistribute it and/or modify 12 * it under the terms of the GNU General Public License as published by 13 * the Free Software Foundation, either version 3 of the License, or 14 * (at your option) any later version. 15 * 16 * This program is distributed in the hope that it will be useful, 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 * GNU General Public License for more details. 20 * 21 * You should have received a copy of the GNU General Public License 22 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 23 24 #include "planets/htrdr_planets_args.h" 25 #include "planets/htrdr_planets_source.h" 26 27 #include "core/htrdr.h" 28 #include "core/htrdr_ran_wlen_discrete.h" 29 30 #include <rsys/math.h> 31 #include <rsys/mem_allocator.h> 32 33 #include <stdio.h> 34 35 static void 36 write_per_wlen_radiances 37 (FILE* fp, 38 const size_t pagesize, 39 const size_t size, 40 const size_t szelmt, 41 const size_t alelmt) 42 { 43 const char byte = 0; 44 size_t i; 45 46 CHK(fp); 47 48 /* Header */ 49 CHK(fwrite(&pagesize, sizeof(pagesize), 1, fp) == 1); 50 CHK(fwrite(&size, sizeof(size), 1, fp) == 1); 51 CHK(fwrite(&szelmt, sizeof(szelmt), 1, fp) == 1); 52 CHK(fwrite(&alelmt, sizeof(alelmt), 1, fp) == 1); 53 54 /* Padding */ 55 CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize), SEEK_SET) == 0); 56 57 FOR_EACH(i, 0, size) { 58 const double w = (double)i; 59 const double L = (double)(100 + i); 60 61 CHK(fwrite(&w, sizeof(w), 1, fp) == 1); 62 CHK(fwrite(&L, sizeof(L), 1, fp) == 1); 63 } 64 65 /* Padding. Write one char to position the EOF indicator */ 66 CHK(fseek(fp, (long)ALIGN_SIZE((size_t)ftell(fp), pagesize)-1, SEEK_SET) == 0); 67 CHK(fwrite(&byte, sizeof(byte), 1, fp) == 1); 68 69 CHK(fflush(fp) == 0); 70 } 71 72 static void 73 test_spectrum(struct htrdr* htrdr) 74 { 75 struct htrdr_planets_source_args source_args = HTRDR_PLANETS_SOURCE_ARGS_NULL; 76 struct htrdr_planets_source_spectrum spectrum; 77 struct htrdr_planets_source* source = NULL; 78 79 FILE* fp = NULL; 80 char rnrl_filename[] = "rnrl.bin"; 81 double range[2]; 82 double w, L; 83 84 CHK(fp = fopen(rnrl_filename, "w")); 85 write_per_wlen_radiances(fp, 4096, 10, 16, 16); 86 CHK(fclose(fp) == 0); 87 88 source_args.rnrl_filename = rnrl_filename; 89 source_args.longitude = 0; 90 source_args.latitude = 0; 91 source_args.distance = 0; 92 source_args.radius = 1e8; 93 source_args.temperature = -1; 94 CHK(htrdr_planets_source_create(htrdr, &source_args, &source) == RES_OK); 95 CHK(htrdr_planets_source_does_radiance_vary_spectrally(source) == 1); 96 CHK(htrdr_planets_source_get_spectral_range(source, range) == RES_OK); 97 CHK(range[0] == 0); 98 CHK(range[1] == 9); 99 100 range[0] = 0; range[1] = 10; 101 CHK(htrdr_planets_source_get_spectrum(source, range, &spectrum) == RES_BAD_ARG); 102 103 range[0] = 1; range[1] = 3; 104 CHK(htrdr_planets_source_get_spectrum(source, range, &spectrum) == RES_OK); 105 CHK(spectrum.source == source); 106 CHK(spectrum.range[0] == 1); 107 CHK(spectrum.range[1] == 3); 108 CHK(spectrum.size == 3); 109 110 htrdr_planets_source_spectrum_at(&spectrum, 0, &w, &L); 111 CHK(w == 1 && L == 101); 112 htrdr_planets_source_spectrum_at(&spectrum, 1, &w, &L); 113 CHK(w == 2 && L == 102); 114 htrdr_planets_source_spectrum_at(&spectrum, 2, &w, &L); 115 CHK(w == 3 && L == 103); 116 117 range[0] = 1.7; range[1] = 1.95; 118 CHK(htrdr_planets_source_get_spectrum(source, range, &spectrum) == RES_OK); 119 CHK(spectrum.source == source); 120 CHK(spectrum.range[0] = 1.7); 121 CHK(spectrum.range[1] = 1.95); 122 CHK(spectrum.size == 2); 123 htrdr_planets_source_spectrum_at(&spectrum, 0, &w, &L); 124 CHK(w == 1.7 && eq_eps(L, 101.7, 1.e-6)); 125 htrdr_planets_source_spectrum_at(&spectrum, 1, &w, &L); 126 CHK(w == 1.95 && eq_eps(L, 101.95, 1.e-6)); 127 128 range[0] = 2; range[1] = 2.01; 129 CHK(htrdr_planets_source_get_spectrum(source, range, &spectrum) == RES_OK); 130 CHK(spectrum.size == 2); 131 htrdr_planets_source_spectrum_at(&spectrum, 0, &w, &L); 132 CHK(w == 2 && L == 102); 133 htrdr_planets_source_spectrum_at(&spectrum, 1, &w, &L); 134 CHK(w == 2.01 && eq_eps(L, 102.01, 1.e-6)); 135 136 range[0] = 5.1; range[1] = 6; 137 CHK(htrdr_planets_source_get_spectrum(source, range, &spectrum) == RES_OK); 138 CHK(spectrum.size == 2); 139 htrdr_planets_source_spectrum_at(&spectrum, 0, &w, &L); 140 CHK(w == 5.1 && eq_eps(L, 105.1, 1.e-6)); 141 htrdr_planets_source_spectrum_at(&spectrum, 1, &w, &L); 142 CHK(w == 6 && L == 106); 143 144 range[0] = 7.5; range[1] = 9; 145 CHK(htrdr_planets_source_get_spectrum(source, range, &spectrum) == RES_OK); 146 CHK(spectrum.size == 3); 147 htrdr_planets_source_spectrum_at(&spectrum, 0, &w, &L); 148 CHK(w == 7.5 && eq_eps(L, 107.5, 1.e-6)); 149 htrdr_planets_source_spectrum_at(&spectrum, 1, &w, &L); 150 CHK(w == 8 && L == 108); 151 htrdr_planets_source_spectrum_at(&spectrum, 2, &w, &L); 152 CHK(w == 9 && L == 109); 153 154 range[0] = 0.9; range[1] = 7.456; 155 CHK(htrdr_planets_source_get_spectrum(source, range, &spectrum) == RES_OK); 156 CHK(spectrum.size == 9); 157 htrdr_planets_source_spectrum_at(&spectrum, 0, &w, &L); 158 CHK(w == 0.9 && eq_eps(L, 100.9, 1.e-6)); 159 htrdr_planets_source_spectrum_at(&spectrum, 1, &w, &L); 160 CHK(w == 1 && eq_eps(L, 101, 1.e-6)); 161 htrdr_planets_source_spectrum_at(&spectrum, 2, &w, &L); 162 CHK(w == 2 && eq_eps(L, 102, 1.e-6)); 163 htrdr_planets_source_spectrum_at(&spectrum, 3, &w, &L); 164 CHK(w == 3 && eq_eps(L, 103, 1.e-6)); 165 htrdr_planets_source_spectrum_at(&spectrum, 4, &w, &L); 166 CHK(w == 4 && eq_eps(L, 104, 1.e-6)); 167 htrdr_planets_source_spectrum_at(&spectrum, 5, &w, &L); 168 CHK(w == 5 && eq_eps(L, 105, 1.e-6)); 169 htrdr_planets_source_spectrum_at(&spectrum, 6, &w, &L); 170 CHK(w == 6 && eq_eps(L, 106, 1.e-6)); 171 htrdr_planets_source_spectrum_at(&spectrum, 7, &w, &L); 172 CHK(w == 7 && eq_eps(L, 107, 1.e-6)); 173 htrdr_planets_source_spectrum_at(&spectrum, 8, &w, &L); 174 CHK(w == 7.456 && eq_eps(L, 107.456, 1.e-6)); 175 176 htrdr_planets_source_ref_put(source); 177 } 178 179 static void 180 test_spectrum_1_entry(struct htrdr* htrdr) 181 { 182 struct htrdr_ran_wlen_discrete_create_args distrib_args = 183 HTRDR_RAN_WLEN_DISCRETE_CREATE_ARGS_NULL; 184 struct htrdr_ran_wlen_discrete* distrib = NULL; 185 186 struct htrdr_planets_source_args source_args = HTRDR_PLANETS_SOURCE_ARGS_NULL; 187 struct htrdr_planets_source_spectrum spectrum; 188 struct htrdr_planets_source* source = NULL; 189 190 FILE* fp = NULL; 191 char rnrl_filename[] = "rnrl.bin"; 192 double range[2] = {0,0}; 193 double lambda = 0; 194 double pdf = 0; 195 196 CHK(fp = fopen(rnrl_filename, "w")); 197 write_per_wlen_radiances(fp, 4096, 1, 16, 16); 198 CHK(fclose(fp) == 0); 199 200 source_args.rnrl_filename = rnrl_filename; 201 source_args.longitude = 0; 202 source_args.latitude = 0; 203 source_args.distance = 0; 204 source_args.radius = 1; 205 source_args.temperature = -1; 206 CHK(htrdr_planets_source_create(htrdr, &source_args, &source) == RES_OK); 207 208 range[0] = range[1] = 0; 209 CHK(htrdr_planets_source_get_spectrum(source, range, &spectrum) == RES_OK); 210 211 distrib_args.get = htrdr_planets_source_spectrum_at; 212 distrib_args.nwavelengths = spectrum.size; 213 distrib_args.context = &spectrum; 214 CHK(htrdr_ran_wlen_discrete_create(htrdr, &distrib_args, &distrib) == RES_OK); 215 216 lambda = htrdr_ran_wlen_discrete_sample(distrib, 0.3, 0.5, &pdf); 217 CHK(lambda == 0); 218 CHK(pdf == 1); 219 220 htrdr_planets_source_ref_put(source); 221 htrdr_ran_wlen_discrete_ref_put(distrib); 222 } 223 224 static void 225 test_spectrum_fail(struct htrdr* htrdr) 226 { 227 struct htrdr_planets_source_args source_args = HTRDR_PLANETS_SOURCE_ARGS_NULL; 228 struct htrdr_planets_source* source = NULL; 229 FILE* fp = NULL; 230 char rnrl_filename[] = "rnrl.bin"; 231 double w, L; 232 233 source_args.rnrl_filename = rnrl_filename; 234 source_args.longitude = 0; 235 source_args.latitude = 0; 236 source_args.distance = 0; 237 source_args.radius = 1e8; 238 source_args.temperature = -1; 239 240 /* Wrong item size */ 241 CHK(fp = fopen(rnrl_filename, "w")); 242 write_per_wlen_radiances(fp, 4096, 10, 8, 16); 243 CHK(fclose(fp) == 0); 244 CHK(htrdr_planets_source_create(htrdr, &source_args, &source) == RES_BAD_ARG); 245 246 /* Wrong item alignment */ 247 CHK(fp = fopen(rnrl_filename, "w")); 248 write_per_wlen_radiances(fp, 4096, 10, 16, 32); 249 CHK(fclose(fp) == 0); 250 CHK(htrdr_planets_source_create(htrdr, &source_args, &source) == RES_BAD_ARG); 251 252 CHK(fp = fopen(rnrl_filename, "w")); 253 write_per_wlen_radiances(fp, 4096, 4, 16, 16); 254 255 /* Overwrite sorted items by unsorted items */ 256 CHK(fseek(fp, 4096, SEEK_SET) == 0); 257 w = 10; L = 1; 258 CHK(fwrite(&w, sizeof(w), 1, fp) == 1); 259 CHK(fwrite(&L, sizeof(L), 1, fp) == 1); 260 w = 11; L = 2; 261 CHK(fwrite(&w, sizeof(w), 1, fp) == 1); 262 CHK(fwrite(&L, sizeof(L), 1, fp) == 1); 263 w = 9; L = 3; 264 CHK(fwrite(&w, sizeof(w), 1, fp) == 1); 265 CHK(fwrite(&L, sizeof(L), 1, fp) == 1); 266 w = 12; L = 4; 267 CHK(fwrite(&w, sizeof(w), 1, fp) == 1); 268 CHK(fwrite(&L, sizeof(L), 1, fp) == 1); 269 CHK(fclose(fp) == 0); 270 271 /* Unsorted items */ 272 CHK(htrdr_planets_source_create(htrdr, &source_args, &source) == RES_BAD_ARG); 273 } 274 275 static void 276 test_spectrum_from_files(struct htrdr* htrdr, int argc, char** argv) 277 { 278 struct htrdr_ran_wlen_discrete_create_args distrib_args = 279 HTRDR_RAN_WLEN_DISCRETE_CREATE_ARGS_NULL; 280 struct htrdr_ran_wlen_discrete* distrib = NULL; 281 282 struct htrdr_planets_source_args source_args = HTRDR_PLANETS_SOURCE_ARGS_NULL; 283 struct htrdr_planets_source_spectrum spectrum = HTRDR_PLANETS_SOURCE_SPECTRUM_NULL; 284 struct htrdr_planets_source* source = NULL; 285 int i; 286 287 source_args.longitude = 0; 288 source_args.latitude = 0; 289 source_args.distance = 0; 290 source_args.radius = 1e8; 291 source_args.temperature = -1; 292 293 FOR_EACH(i, 1, argc) { 294 double range[2]; 295 double lambda, pdf; 296 source_args.rnrl_filename = argv[i]; 297 298 CHK(htrdr_planets_source_create(htrdr, &source_args, &source) == RES_OK); 299 CHK(htrdr_planets_source_does_radiance_vary_spectrally(source)); 300 CHK(htrdr_planets_source_get_spectral_range(source, range) == RES_OK); 301 302 range[0] = 250; 303 range[1] = 850; 304 CHK(htrdr_planets_source_get_spectrum(source, range, &spectrum) == RES_OK); 305 306 printf("`%s' stores %lu entries between [%g, %g] nm\n", 307 argv[i], spectrum.size, SPLIT2(range)); 308 309 distrib_args.get = htrdr_planets_source_spectrum_at; 310 distrib_args.nwavelengths = spectrum.size; 311 distrib_args.context = &spectrum; 312 CHK(htrdr_ran_wlen_discrete_create(htrdr, &distrib_args, &distrib) == RES_OK); 313 314 lambda = htrdr_ran_wlen_discrete_sample(distrib, 0.3, 0.5, &pdf); 315 printf("lambda = %g nm; pdf = %f nm⁻¹\n", lambda, pdf); 316 317 htrdr_planets_source_ref_put(source); 318 htrdr_ran_wlen_discrete_ref_put(distrib); 319 } 320 } 321 322 int 323 main(int argc, char** argv) 324 { 325 struct htrdr_args args = HTRDR_ARGS_DEFAULT; 326 struct htrdr* htrdr = NULL; 327 size_t memsz = 0; 328 329 args.verbose = 1; 330 htrdr_mpi_init(argc, argv); 331 CHK(htrdr_create(NULL, &args, &htrdr) == RES_OK); 332 333 memsz = MEM_ALLOCATED_SIZE(htrdr_get_allocator(htrdr)); 334 335 if(argc > 1) { 336 test_spectrum_from_files(htrdr, argc, argv); 337 } else { 338 test_spectrum(htrdr); 339 test_spectrum_1_entry(htrdr); 340 test_spectrum_fail(htrdr); 341 } 342 343 CHK(MEM_ALLOCATED_SIZE(htrdr_get_allocator(htrdr)) == memsz); 344 345 htrdr_ref_put(htrdr); 346 htrdr_mpi_finalize(); 347 return 0; 348 }