test_solstice_simulation_time.c (13791B)
1 /* Copyright (C) 2018-2026 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2016-2018 CNRS 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 _XOPEN_SOURCE /* strptime support */ 18 #define _POSIX_C_SOURCE 200809L /* mkstemp support */ 19 20 #include <star/scem.h> 21 22 #include <rsys/rsys.h> 23 #include <rsys/math.h> 24 #include <rsys/double2.h> 25 26 #include <time.h> 27 #include <stdio.h> 28 #include <stdlib.h> 29 #include <string.h> 30 31 enum side { 32 FRONT, 33 BACK 34 }; 35 36 enum global_result_type { 37 GLOBAL_POTENTIAL, 38 GLOBAL_ABSORBED, 39 GLOBAL_COS, 40 GLOBAL_SHADOW, 41 GLOBAL_MISSING, 42 GLOBAL_ATMOSPHERE, 43 GLOBAL_REFLECTIVITY, 44 GLOBAL_RESULTS_COUNT__ 45 }; 46 47 enum receiver_result_type { 48 FIRST_RECEIVER_RESULT, 49 FRONT_ABSORBED_FLUX = FIRST_RECEIVER_RESULT, 50 FRONT_INCOMING_FLUX, 51 FRONT_ABSORBED_FIELD_GAIN, 52 FRONT_ABSORBED_ATM_GAIN, 53 FRONT_EFFICIENCY, 54 BACK_ABSORBED_FLUX, 55 BACK_INCOMING_FLUX, 56 BACK_ABSORBED_FIELD_GAIN, 57 BACK_ABSORBED_ATM_GAIN, 58 BACK_EFFICIENCY, 59 RECEIVER_RESULTS_COUNT__ 60 }; 61 62 enum primary_result_type { 63 FIRST_PRIMARY_RESULT, 64 PRIMARY_COS = FIRST_PRIMARY_RESULT, 65 PRIMARY_SHADOW, 66 PRIMARY_RESULTS_COUNT__ 67 }; 68 69 struct counts { 70 unsigned long global, receiver, primary, realisation, failed; 71 }; 72 73 static int 74 counts_ok(const struct counts* ref, const struct counts* c) 75 { 76 CHK(ref->global == GLOBAL_RESULTS_COUNT__); 77 CHK(c->global == GLOBAL_RESULTS_COUNT__); 78 return ref->receiver == c->receiver 79 && ref->primary == c->primary 80 && ref->failed >= c->failed; 81 } 82 83 #define MAX_LINE_LEN 2048 84 85 static const char 86 sundir_header [] = "#--- Sun direction:"; 87 static const char 88 suntime_header [] = "#--- Sun location and time:"; 89 90 #define IS_NEW_BLOCK(Line, Header) (!strncmp((Line), (Header), strlen(Header))) 91 92 static int 93 read_line(char* line, size_t max_line_len, FILE* stream) 94 { 95 ASSERT(stream && line && max_line_len); 96 line = fgets(line, (int)max_line_len, stream); 97 if(!line) return 0; 98 CHK(strlen(line) + 1 < max_line_len); 99 return 1; 100 } 101 102 static int 103 get_angles_and_counts 104 (FILE* file, 105 double angles[2], 106 struct counts* counts) 107 { 108 char line[MAX_LINE_LEN]; 109 int r; 110 111 CHK(file != NULL); 112 CHK(angles != NULL); 113 CHK(counts != NULL); 114 115 /* Get sun dir */ 116 r = read_line(line, sizeof(line), file); 117 if (!r) { 118 CHK(feof(file) == 1); 119 return 0; 120 } 121 CHK(IS_NEW_BLOCK(line, sundir_header) == 1); 122 CHK(sscanf(line+strlen(sundir_header), "%lg%lg", &angles[0], &angles[1]) == 2); 123 124 /* Get counts */ 125 CHK(read_line(line, sizeof(line), file) == 1); 126 CHK( 127 sscanf(line, 128 "%lu %lu %lu %lu %lu", 129 &counts->global, &counts->receiver, &counts->primary, 130 &counts->realisation, &counts->failed) == 5); 131 return 1; 132 } 133 134 static int 135 get_dir_and_counts 136 (FILE* file, 137 double dir[3], 138 struct counts* counts) 139 { 140 char line[MAX_LINE_LEN]; 141 int r, n; 142 143 CHK(file != NULL); 144 CHK(dir != NULL); 145 CHK(counts != NULL); 146 147 /* Get sun dir */ 148 r = read_line(line, sizeof(line), file); 149 if (!r) { 150 CHK(feof(file) == 1); 151 return 0; 152 } 153 CHK(IS_NEW_BLOCK(line, suntime_header) == 1); 154 n = sscanf(line+strlen(suntime_header), 155 "%*g %*g %*d-%*d-%*dT%*d:%*d:%*d (%lg %lg %lg)", 156 dir, dir+1, dir+2); 157 CHK(n == 3); 158 159 /* Get counts */ 160 CHK(read_line(line, sizeof(line), file) == 1); 161 CHK( 162 sscanf(line, 163 "%lu %lu %lu %lu %lu", 164 &counts->global, &counts->receiver, &counts->primary, 165 &counts->realisation, &counts->failed) == 5); 166 return 1; 167 } 168 169 static void 170 read_global(FILE* file, double* E, double* SE) 171 { 172 char line[MAX_LINE_LEN]; 173 CHK(read_line(line, sizeof(line), file) == 1); 174 CHK(sscanf(line, "%lg %lg", E, SE) == 2); 175 } 176 177 static void 178 read_recv(FILE* file, char name[], double E[], double SE[]) 179 { 180 char line[MAX_LINE_LEN]; 181 182 CHK(file != NULL); 183 CHK(name != NULL); 184 CHK(E != NULL); 185 CHK(SE != NULL); 186 187 CHK(read_line(line, sizeof(line), file) == 1); 188 CHK( 189 sscanf(line, 190 "%s %*u %*g " 191 "%lg %lg %lg %lg %lg %lg %lg %lg %lg %lg " 192 "%lg %lg %lg %lg %lg %lg %lg %lg %lg %lg", 193 name, /* ID, area */ 194 &E[FRONT_ABSORBED_FLUX], &SE[FRONT_ABSORBED_FLUX], 195 &E[FRONT_INCOMING_FLUX], &SE[FRONT_INCOMING_FLUX], 196 &E[FRONT_ABSORBED_FIELD_GAIN], &SE[FRONT_ABSORBED_FIELD_GAIN], 197 &E[FRONT_ABSORBED_ATM_GAIN], &SE[FRONT_ABSORBED_ATM_GAIN], 198 &E[FRONT_EFFICIENCY], &SE[FRONT_EFFICIENCY], 199 &E[BACK_ABSORBED_FLUX], &SE[BACK_ABSORBED_FLUX], 200 &E[BACK_INCOMING_FLUX], &SE[BACK_INCOMING_FLUX], 201 &E[BACK_ABSORBED_FIELD_GAIN], &SE[BACK_ABSORBED_FIELD_GAIN], 202 &E[BACK_ABSORBED_ATM_GAIN], &SE[BACK_ABSORBED_ATM_GAIN], 203 &E[BACK_EFFICIENCY], &SE[BACK_EFFICIENCY]) == 204 2 * RECEIVER_RESULTS_COUNT__ + 1); 205 } 206 207 static void 208 read_primary 209 (FILE* file, char name[], double* area, double E[], double SE[]) 210 { 211 char line[MAX_LINE_LEN]; 212 213 CHK(file != NULL); 214 CHK(area != NULL); 215 CHK(E != NULL); 216 CHK(SE != NULL); 217 218 CHK(read_line(line, sizeof(line), file) == 1); 219 CHK( 220 sscanf(line, 221 "%s %*u " 222 "%lg %*u " 223 "%lg %lg %lg %lg\n", 224 name, /* ID */ 225 area, /* count, */ 226 &E[PRIMARY_COS], &SE[PRIMARY_COS], 227 &E[PRIMARY_SHADOW], &SE[PRIMARY_SHADOW]) == 228 2 * PRIMARY_RESULTS_COUNT__ + 2); 229 } 230 231 232 static void 233 read_recvXprim 234 (FILE* file, 235 unsigned long* rcv_id, 236 unsigned long* prim_id, 237 double E[], 238 double SE[]) 239 { 240 char line[MAX_LINE_LEN]; 241 242 CHK(file != NULL); 243 CHK(rcv_id != NULL); 244 CHK(prim_id != NULL); 245 CHK(E != NULL); 246 CHK(SE != NULL); 247 248 CHK(read_line(line, sizeof(line), file) == 1); 249 CHK( 250 sscanf(line, 251 "%lu %lu " 252 "%lg %lg %lg %lg %lg %lg %lg %lg " 253 "%lg %lg %lg %lg %lg %lg %lg %lg", 254 rcv_id, prim_id, 255 &E[FRONT_ABSORBED_FLUX], &SE[FRONT_ABSORBED_FLUX], 256 &E[FRONT_INCOMING_FLUX], &SE[FRONT_INCOMING_FLUX], 257 &E[FRONT_ABSORBED_FIELD_GAIN], &SE[FRONT_ABSORBED_FIELD_GAIN], 258 &E[FRONT_ABSORBED_ATM_GAIN], &SE[FRONT_ABSORBED_ATM_GAIN], 259 &E[BACK_ABSORBED_FLUX], &SE[BACK_ABSORBED_FLUX], 260 &E[BACK_INCOMING_FLUX], &SE[BACK_INCOMING_FLUX], 261 &E[BACK_ABSORBED_FIELD_GAIN], &SE[BACK_ABSORBED_FIELD_GAIN], 262 &E[BACK_ABSORBED_ATM_GAIN], &SE[BACK_ABSORBED_ATM_GAIN]) == 263 2 * (RECEIVER_RESULTS_COUNT__ - 2 /* efficiencies not read */) + 2); 264 } 265 266 static void 267 compute_estimate_intersection 268 (double intersection[2], 269 const double scale, 270 const double E0, 271 const double SE0, 272 const double E1, 273 const double SE1) 274 { 275 double interval0[2], interval1[2]; 276 CHK(scale > 0); 277 interval0[0] = E0 - scale*SE0; 278 interval0[1] = E0 + scale*SE0; 279 interval1[0] = E1 - scale*SE1; 280 interval1[1] = E1 + scale*SE1; 281 intersection[0] = MMAX(interval0[0], interval1[0]); 282 intersection[1] = MMIN(interval0[1], interval1[1]); 283 } 284 285 static void 286 check_estimate 287 (double ref_E, 288 double ref_SE, 289 double test_E, 290 double test_SE) 291 { 292 if(ref_E == -1) { 293 CHK(ref_SE == -1); 294 CHK(test_E == -1); 295 CHK(test_SE == -1); 296 } else { 297 double interval[2]; 298 CHK(ref_SE >= 0); 299 CHK(test_E >= 0); 300 CHK(test_SE >= 0); 301 if(ref_SE == 0) ref_SE = ref_E / 1000.0; 302 if(test_SE == 0) test_SE = test_E / 1000.0; 303 compute_estimate_intersection(interval, 2, ref_E, ref_SE, test_E, test_SE); 304 CHK(interval[0] <= interval[1]); 305 } 306 } 307 308 static void 309 check_1_reference 310 (FILE* ref_file, 311 FILE* test_file, 312 const struct counts* counts) 313 { 314 unsigned n; 315 316 CHK(ref_file != NULL); 317 CHK(test_file != NULL); 318 CHK(counts != NULL); 319 320 /* both files' pointer are just past the new bloc header */ 321 322 for(n = 0; n < counts->global; n++) { 323 double reference_E, reference_SE, test_E, test_SE; 324 read_global(ref_file, &reference_E, &reference_SE); 325 read_global(test_file, &test_E, &test_SE); 326 check_estimate(reference_E, reference_SE, test_E, test_SE); 327 } 328 for(n = 0; n < counts->receiver; n++) { 329 char ref_rcv_name[MAX_LINE_LEN], test_rcv_name[MAX_LINE_LEN]; 330 double reference_E[RECEIVER_RESULTS_COUNT__]; 331 double reference_SE[RECEIVER_RESULTS_COUNT__]; 332 double test_E[RECEIVER_RESULTS_COUNT__]; 333 double test_SE[RECEIVER_RESULTS_COUNT__]; 334 enum receiver_result_type r; 335 336 read_recv(ref_file, ref_rcv_name, reference_E, reference_SE); 337 read_recv(test_file, test_rcv_name, test_E, test_SE); 338 CHK(strcmp(ref_rcv_name, test_rcv_name) == 0); 339 FOR_EACH(r, FIRST_RECEIVER_RESULT, RECEIVER_RESULTS_COUNT__) { 340 check_estimate(reference_E[r], reference_SE[r], test_E[r], test_SE[r]); 341 } 342 } 343 for(n = 0; n < counts->primary; n++) { 344 char ref_prim_name[MAX_LINE_LEN], test_prim_name[MAX_LINE_LEN]; 345 double reference_E[PRIMARY_RESULTS_COUNT__]; 346 double reference_SE[PRIMARY_RESULTS_COUNT__]; 347 double test_E[PRIMARY_RESULTS_COUNT__]; 348 double test_SE[PRIMARY_RESULTS_COUNT__]; 349 double ref_area, test_area; 350 enum primary_result_type r; 351 352 read_primary(ref_file, ref_prim_name, &ref_area, reference_E, reference_SE); 353 read_primary(test_file, test_prim_name, &test_area, test_E, test_SE); 354 check_estimate(ref_area, 0, test_area, 0); 355 CHK(strcmp(ref_prim_name, test_prim_name) == 0); 356 FOR_EACH(r, FIRST_PRIMARY_RESULT, PRIMARY_RESULTS_COUNT__) { 357 check_estimate(reference_E[r], reference_SE[r], test_E[r], test_SE[r]); 358 } 359 } 360 for(n = 0; n < counts->receiver * counts->primary; n++) { 361 double reference_E[RECEIVER_RESULTS_COUNT__]; 362 double reference_SE[RECEIVER_RESULTS_COUNT__]; 363 double test_E[RECEIVER_RESULTS_COUNT__]; 364 double test_SE[RECEIVER_RESULTS_COUNT__]; 365 unsigned long ref_rcv_id, ref_prim_id; 366 unsigned long test_rcv_id, test_prim_id; 367 368 enum receiver_result_type r; 369 read_recvXprim(ref_file, &ref_rcv_id, &ref_prim_id, reference_E, reference_SE); 370 read_recvXprim(test_file, &test_rcv_id, &test_prim_id, test_E, test_SE); 371 /* we rely on the order of outputs */ 372 CHK(ref_rcv_id == test_rcv_id); 373 CHK(ref_prim_id == test_prim_id); 374 FOR_EACH(r, FIRST_RECEIVER_RESULT, RECEIVER_RESULTS_COUNT__) { 375 if (r == FRONT_EFFICIENCY || r == BACK_EFFICIENCY) 376 continue; /* not read */ 377 check_estimate(reference_E[r], reference_SE[r], test_E[r], test_SE[r]); 378 } 379 } 380 } 381 382 static FINLINE int 383 create_tmp_file(char* name, const size_t max_sizeof_name) 384 { 385 const char* template = "solstice_tmp_file_XXXXXX"; 386 int fd; 387 CHK(name != NULL); 388 CHK(strlen(template)+1 <= max_sizeof_name-1); 389 strcpy(name, template); 390 fd = mkstemp(name); 391 CHK(fd != -1); 392 return fd; 393 } 394 395 /* Scem an solstice agree on elevation. 396 * On the other hand, scem and solstice convention on azimuth differ: 397 * - scem CW with N=0 398 * - solstice CCW with E=0 */ 399 static FINLINE void 400 scem_pos_to_solstice_pos 401 (const struct scem_sun_pos* scem, double solstice[2]) 402 { 403 ASSERT(scem && solstice); 404 solstice[0] = fmod(360 + 90 - MRAD2DEG(scem->azimuth), 360); 405 solstice[1] = MRAD2DEG(scem->elevation); 406 } 407 408 static void 409 do_check(const char* binary, const char* dir, const char* base_name) 410 { 411 struct counts time_counts, dir_counts; 412 int n, h; 413 int err; 414 double lat, lon; 415 char cmd[512]; 416 char time_file_name[128]; 417 char dir_file_name[128]; 418 char sdate[64]; 419 char* p; 420 struct tm date; 421 struct scem_sun_pos scem_pos; 422 struct scem_location loc; 423 enum scem_sun_algo algorithm; 424 res_T res; 425 double sun_dir[3], dir_angles[2], solstice_pos[2]; 426 FILE* dir_file; 427 FILE* time_file; 428 int dir_fd, time_fd; 429 430 lat = -90 + 180 * ((double)rand() / (double)RAND_MAX); 431 lon = 15 * ((double)rand() / (double)RAND_MAX); /* To allow to sample UTC h */ 432 h = (int)(9 + 10 * ((double)rand() / (double)RAND_MAX)); 433 algorithm = SCEM_SUN_PSA; 434 435 snprintf(sdate, sizeof(sdate), "2020-05-05T%d:00:00", h); 436 p = strptime(sdate, "%Y-%m-%dT%H:%M:%S", &date); 437 CHK(p && *p == '\0'); 438 loc.latitude = lat; loc.longitude = lon; 439 440 time_fd = create_tmp_file(time_file_name, sizeof(time_file_name)); 441 time_file = fdopen(time_fd, "r"); 442 CHK(time_file != NULL); 443 444 n = snprintf(cmd, sizeof(cmd), 445 "%s -o %s -f -L %g,%g -T %s -n %lu -R %s%s_receiver.yaml %s%s.yaml", 446 binary, time_file_name, lat, lon, sdate, 10000L, 447 dir, base_name, dir, base_name); 448 CHK((unsigned)n < sizeof(cmd)); 449 450 err = system(cmd); 451 CHK(err == 0); 452 453 get_dir_and_counts(time_file, sun_dir, &time_counts); 454 455 dir_fd = create_tmp_file(dir_file_name, sizeof(dir_file_name)); 456 dir_file = fdopen(dir_fd, "r"); 457 CHK(dir_file != NULL); 458 459 res = scem_sun_position_from_earth(&date, &loc, algorithm, &scem_pos); 460 CHK(res == RES_OK); 461 scem_pos_to_solstice_pos(&scem_pos, solstice_pos); 462 n = snprintf(cmd, sizeof(cmd), 463 "%s -o %s -f -D %g,%g -n %lu -R %s%s_receiver.yaml %s%s.yaml", 464 binary, dir_file_name, 465 SPLIT2(solstice_pos), 466 10000L, dir, base_name, dir, base_name); 467 CHK((unsigned)n < sizeof(cmd)); 468 469 err = system(cmd); 470 CHK(err == 0); 471 472 get_angles_and_counts(dir_file, dir_angles, &dir_counts); 473 474 CHK(counts_ok(&time_counts, &dir_counts) == 1); 475 check_1_reference(time_file, dir_file, &time_counts); 476 477 fclose(time_file); 478 fclose(dir_file); 479 remove(time_file_name); 480 remove(dir_file_name); 481 } 482 483 int 484 main(int argc, char** argv) 485 { 486 int err = 0; 487 488 if(argc != 4) { 489 printf("Usage: %s <solstice-binary> <file-path> <file-base-name>\n", argv[0]); 490 goto error; 491 } 492 493 do_check(argv[1], argv[2], argv[3]); 494 495 exit: 496 return err; 497 error: 498 err = 1; 499 goto exit; 500 } 501