stardis_smeteo.c (15647B)
1 /* Copyright (C) 2025 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is dismshbuted in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #define _POSIX_C_SOURCE 200112L /* localtime_r & getopt support */ 17 18 #include "smeteo.h" 19 #include "stardis_smeteo.h" 20 #include "stardis_smeteo_library.h" 21 22 #include <star/scem.h> 23 24 #include <rsys/algorithm.h> 25 #include <rsys/cstr.h> 26 #include <rsys/math.h> 27 #include <rsys/mem_allocator.h> 28 29 /* FIXME required by the UINT_MAX constant used in the stardis-prog-properties 30 * header. It should by stardis, not by the caller. */ 31 #include <limits.h> 32 33 #include <time.h> /* localtime_r */ 34 #include <unistd.h> /* getopt */ 35 36 #define EARTH_TO_SUN_DST 149.5978707e9 /* [m] */ 37 38 enum flux_cpnt { 39 FLUX_SW = BIT(0), /* Shortwave flux */ 40 FLUX_LATENT = BIT(1), /* Latent flux */ 41 FLUX_ALL = ~0, /* All components are handled in flux computation */ 42 FLUX_NONE = 0 /* No flux is imposed */ 43 }; 44 45 struct args { 46 double temperature; /* Take precedence on meteorological data [K] */ 47 int convection; /* Enable convection */ 48 int imposed_flux; /* Imposed flux mask. combination of flux_cpnt */ 49 }; 50 #define ARGS_DEFAULT__ {STARDIS_TEMPERATURE_NONE, 0, 0} 51 static const struct args ARGS_DEFAULT = ARGS_DEFAULT__; 52 53 struct boundary_condition { 54 struct args args; 55 struct stardis_smeteo_lib* lib; 56 struct stardis_smeteo_lib_desc lib_desc; 57 }; 58 static const struct boundary_condition BOUNDARY_CONDITION_NULL = { 59 ARGS_DEFAULT__, NULL, STARDIS_SMETEO_LIB_DESC_NULL__ 60 }; 61 62 /******************************************************************************* 63 * Helper functions 64 ******************************************************************************/ 65 static void 66 usage(FILE* stream, const char* name) 67 { 68 ASSERT(stream && name); 69 fprintf(stream, "usage: %s [-hls] [-t temperature]\n", name); 70 } 71 72 static res_T 73 args_init(struct args* args, int argc, char* argv[]) 74 { 75 int opt = 0; 76 res_T res = RES_OK; 77 ASSERT(args && argc && argv); 78 79 *args = ARGS_DEFAULT; 80 81 optind = 1; 82 while((opt=getopt(argc, argv, "hlst:")) != -1) { 83 switch(opt) { 84 case 'h': args->convection = 1; break; 85 case 'l': args->imposed_flux |= FLUX_LATENT; break; 86 case 's': args->imposed_flux |= FLUX_SW; break; 87 case 't': 88 res = cstr_to_double(optarg, &args->temperature); 89 break; 90 default: res = RES_BAD_ARG; break; 91 } 92 if(res != RES_OK) { 93 if(optarg) { 94 fprintf(stderr, "%s: invalid option argument '%s' -- '%c'\n", 95 argv[0], optarg, opt); 96 } 97 goto error; 98 } 99 } 100 exit: 101 return res; 102 error: 103 usage(stderr, argv[0]); 104 goto exit; 105 } 106 107 static int 108 cmp_day_1850(const void* key, const void* val) 109 { 110 const double day_1850 = *(double*)key; 111 const struct smeteo_entry* entry = val; 112 if(day_1850 < entry->day_1850) return -1; 113 else if(day_1850 > entry->day_1850) return +1; 114 else /* day_1850 == entry->day_1850 */ return 0; 115 } 116 117 static size_t 118 get_meteo_entry_id 119 (const struct boundary_condition* bcond, 120 const double time) /* [s] */ 121 { 122 struct smeteo_entry* found = NULL; 123 double time_1850 = 0; 124 ASSERT(bcond); 125 126 /* Define the input time in relation to January 1, 1850, i.e., in terms of the 127 * number of days that have elapsed since 00:00 on January 1, 1850. 128 * Note that this simple conversion ignored leap seconds. */ 129 time_1850 = time/*[s]*/ / (24/*hours*/*3600/*[s]*/); 130 131 /* Limit to 0, i.e. January 1, 1850, i.e. the initial condition */ 132 time_1850 = time_1850 < 0 ? 0 : time_1850; 133 134 /* Search for the time interval whose center date is greater than or equal to 135 * the input time */ 136 found = search_lower_bound 137 (&time_1850, 138 bcond->lib_desc.smeteo_desc.entries, 139 bcond->lib_desc.smeteo_desc.nentries, 140 sizeof(struct smeteo_entry), 141 cmp_day_1850); 142 143 /* There are no intervals whose center time is greater than the input time. 144 * The time is therefore either in the upper half of the interval or beyond 145 * it. In both cases, return the last interval, which, in the case of an input 146 * time beyond the available time intervals, amounts to repeating the last 147 * interval indefinitely. */ 148 if(!found) { 149 return bcond->lib_desc.smeteo_desc.nentries -1; 150 151 } else { 152 double half_interval_duration = 0; 153 154 /* Define whether the input time is included in the interval found or in 155 * the previous one. If the time elapsed between the input time and the 156 * central time of the current interval is less than half the duration of an 157 * interval, then the time belongs to the interval found. Otherwise, it 158 * belongs to the previous one. 159 * 160 * Note that there is no need to check that the interval found is the first 161 * one and that there is therefore no previous interval, because the input 162 * time has already been truncated to 0, i.e. midnight on January 1, 1850 */ 163 half_interval_duration = bcond->lib_desc.smeteo_desc.entries[0].day_1850; 164 if(found->day_1850 - time_1850 > half_interval_duration) --found; 165 166 /* Check that the input time is within the selected time interval */ 167 ASSERT(time_1850 >= found->day_1850 - half_interval_duration); 168 ASSERT(time_1850 <= found->day_1850 + half_interval_duration); 169 170 /* Calculate the index of the selected interval */ 171 ASSERT((intptr_t)found >= (intptr_t)bcond->lib_desc.smeteo_desc.entries); 172 return (size_t)(found - bcond->lib_desc.smeteo_desc.entries); 173 } 174 } 175 176 static void 177 spherical_to_cartesian_dir 178 (const double azimuth, /* [In radians */ 179 const double elevation, /* In radians */ 180 double dir[3]) 181 { 182 double cos_azimuth; 183 double sin_azimuth; 184 double cos_elevation; 185 double sin_elevation; 186 ASSERT(azimuth >= 0 && azimuth < 2*PI); 187 ASSERT(elevation >= -PI/2.0 && elevation <= PI/2.0); 188 ASSERT(dir); 189 190 cos_azimuth = cos(azimuth); 191 sin_azimuth = sin(azimuth); 192 cos_elevation = cos(elevation); 193 sin_elevation = sin(elevation); 194 195 dir[0] = cos_elevation * cos_azimuth; 196 dir[1] = cos_elevation * sin_azimuth; 197 dir[2] = sin_elevation; 198 } 199 200 static void 201 compute_sun_position 202 (const struct boundary_condition* bcond, 203 const double time, /* UTC +00:00 [s] */ 204 struct scem_sun_pos* sun) 205 { 206 /* Star-CElestial-Mechanics */ 207 struct scem_location pos = SCEM_LOCATION_NULL; 208 209 /* Time and date */ 210 const struct tm* date_utc0 = NULL; 211 size_t i = 0; 212 213 res_T res = RES_OK; 214 ASSERT(bcond && sun); 215 216 /* Get the time associated with the interval. Thus, the position of the sun is 217 * assumed to be constant during the interval. This is necessary to ensure 218 * energy conservation given how solar flux is defined, namely as an average 219 * over the time interval of the solar flux incident perpendicular to a 220 * horizontal surface. */ 221 i = get_meteo_entry_id(bcond, time); 222 date_utc0 = &bcond->lib_desc.smeteo_desc.entries[i].time; 223 224 /* Compute the solar position at the specified location and date */ 225 pos.latitude = bcond->lib_desc.smeteo_desc.latitude; 226 pos.longitude = bcond->lib_desc.smeteo_desc.longitude; 227 res = scem_sun_position_from_earth 228 (date_utc0, &pos, bcond->lib_desc.algo, sun); 229 CHK(res == RES_OK); 230 } 231 232 /******************************************************************************* 233 * Exported symbols 234 ******************************************************************************/ 235 void* 236 stardis_create_data 237 (const struct stardis_description_create_context* ctx, 238 void* lib_data, 239 size_t argc, 240 char* argv[]) 241 { 242 struct boundary_condition* bcond = NULL; 243 res_T res = RES_OK; 244 245 if(!ctx || !lib_data || !argc || !argv) goto error; 246 247 if(!(bcond = mem_calloc(1, sizeof(*bcond)))) goto error; 248 *bcond = BOUNDARY_CONDITION_NULL; 249 250 if((res = args_init(&bcond->args, (int)argc, argv)) != RES_OK) goto error; 251 bcond->lib = lib_data; 252 stardis_smeteo_lib_ref_get(bcond->lib); 253 stardis_smeteo_lib_get_desc(bcond->lib, &bcond->lib_desc); 254 255 exit: 256 return bcond; 257 error: 258 if(bcond) { mem_rm(bcond); bcond = NULL; } 259 goto exit; 260 } 261 262 void 263 stardis_release_data(void* data) 264 { 265 struct boundary_condition* bcond = data; 266 ASSERT(data); 267 268 stardis_smeteo_lib_ref_put(bcond->lib); 269 mem_rm(bcond); 270 } 271 272 double /* [W/K/m^2] */ 273 stardis_convection_coefficient 274 (const struct stardis_interface_fragment* frag, 275 void* data) 276 { 277 struct boundary_condition* bcond = data; 278 ASSERT(frag && data); 279 280 if(!bcond->args.convection) { 281 return 0; 282 283 } else { 284 size_t i = 0; /* Index of the meteo entry including fragment time */ 285 286 i = get_meteo_entry_id(bcond, frag->time); 287 return bcond->lib_desc.smeteo_desc.entries[i].H; 288 } 289 } 290 291 double /* [W/K/m^2] */ 292 stardis_max_convection_coefficient(void* data) 293 { 294 struct boundary_condition* bcond = data; 295 ASSERT(data); 296 return bcond->lib_desc.max_convection_coef; 297 } 298 299 double /* [W/m^2] */ 300 stardis_boundary_flux 301 (const struct stardis_interface_fragment* frag, 302 void* data) 303 { 304 struct boundary_condition* bcond = data; 305 double flux = 0; /* [W/m^2] */ 306 size_t i = 0; 307 ASSERT(frag && data); 308 309 if(bcond->args.imposed_flux == FLUX_NONE) { 310 return STARDIS_FLUX_NONE; 311 } 312 313 i = get_meteo_entry_id(bcond, frag->time); 314 315 if(bcond->args.imposed_flux & FLUX_LATENT) { 316 flux += -bcond->lib_desc.smeteo_desc.entries[i].LE; 317 } 318 319 if(bcond->args.imposed_flux & FLUX_SW) { 320 double SWdn = 0; /* shortwave downward flux [W/m^2] */ 321 SWdn = bcond->lib_desc.smeteo_desc.entries[i].SWdn_direct 322 + bcond->lib_desc.smeteo_desc.entries[i].SWdn_diffuse; 323 324 /* Flux on the ground side, i.e. the downward flux - upward flux */ 325 flux += SWdn - bcond->lib_desc.smeteo_desc.entries[i].SWup; 326 } 327 328 return flux; 329 } 330 331 double /* [K] */ 332 stardis_medium_temperature(const struct stardis_vertex* vtx, void* data) 333 { 334 struct boundary_condition* bcond = data; 335 size_t i = 0; 336 ASSERT(vtx && data); 337 338 i = get_meteo_entry_id(bcond, vtx->time); 339 return bcond->lib_desc.smeteo_desc.entries[i].Tatm; 340 } 341 342 /* Ground emissivity */ 343 double 344 stardis_emissivity 345 (const struct stardis_interface_fragment* frag, 346 const unsigned source_id, 347 void* data) 348 { 349 struct boundary_condition* bcond = data; 350 ASSERT(data); 351 (void)frag; /* Avoid "unused variable" warning */ 352 353 if(source_id == STARDIS_INTERN_SOURCE_ID) { 354 return 1; 355 } else { 356 return 1 - bcond->lib_desc.smeteo_desc.albedo; 357 } 358 } 359 360 /* Specular part of the BRDF of the ground */ 361 double 362 stardis_specular_fraction 363 (const struct stardis_interface_fragment* frag, 364 const unsigned source_id, 365 void* data) 366 { 367 (void)data, (void)frag, (void)source_id; /* Avoid "unused variable" warning */ 368 return 0; /* <=> The ground is purely diffuse */ 369 } 370 371 /* The reference ground temperature is set here to the surface temperature 372 * obtained from meteorological data. */ 373 double /* [K] */ 374 stardis_reference_temperature 375 (const struct stardis_interface_fragment* frag, 376 void* data) 377 { 378 return stardis_boundary_temperature(frag, data); 379 } 380 381 double 382 stardis_calorific_capacity(const struct stardis_vertex* vtx, void* data) 383 { 384 (void)vtx, (void)data; 385 return 1; /* Dummy */ 386 } 387 388 double 389 stardis_volumic_mass(const struct stardis_vertex* vtx, void* data) 390 { 391 (void)vtx, (void)data; 392 return 1; /* Dummy */ 393 } 394 395 /* Ground temperature */ 396 double /* [K] */ 397 stardis_boundary_temperature 398 (const struct stardis_interface_fragment* frag, 399 void* data) 400 { 401 struct boundary_condition* bcond = data; 402 ASSERT(frag && data); 403 404 if(STARDIS_TEMPERATURE_IS_KNOWN(bcond->args.temperature)) { 405 return bcond->args.temperature; 406 } else { 407 size_t i = 0; /* Index of the meteo entry including fragment time */ 408 409 i = get_meteo_entry_id(bcond, frag->time); 410 return bcond->lib_desc.smeteo_desc.entries[i].Tsrf; 411 } 412 } 413 414 double /* [K] */ 415 stardis_radiative_env_temperature 416 (const double time, /* [s] */ 417 const double dir[3], 418 void* data) 419 { 420 struct boundary_condition* bcond = data; 421 size_t i = 0; /* Index of the meteo entry including fragment time */ 422 ASSERT(time && dir && data); 423 (void)dir, (void)data; /* Avoid "unused variable" warning */ 424 425 i = get_meteo_entry_id(bcond, time); 426 return bcond->lib_desc.smeteo_desc.entries[i].Trad; 427 } 428 429 double /* [K] */ 430 stardis_radiative_env_reference_temperature 431 (const double time, /* [s] */ 432 const double dir[3], 433 void* data) 434 { 435 return stardis_radiative_env_temperature(time, dir, data); 436 } 437 438 /* Range of reference temperature, i.e. range of ground temperature from the 439 * entire set of surface temperatures */ 440 double* 441 stardis_t_range(void* data, double range[2] /* [K] */) 442 { 443 struct boundary_condition* bcond = data; 444 ASSERT(data && range); 445 446 range[0] = MMIN(bcond->lib_desc.Tsrf_range[0], bcond->lib_desc.Trad_range[0]); 447 range[1] = MMAX(bcond->lib_desc.Tsrf_range[1], bcond->lib_desc.Trad_range[1]); 448 return range; 449 } 450 451 double* 452 stardis_spherical_source_position 453 (const double time, /* [s] */ 454 double sun_pos[3], /* [m] */ 455 void* data) 456 { 457 struct scem_sun_pos sun = SCEM_SUN_POS_NULL; 458 double sun_dir[3] = {0,0,0}; 459 ASSERT(sun_pos); 460 461 compute_sun_position(data, time, &sun); 462 463 /* Convert the sun's position to Cartesian coords as expected by Stardis */ 464 spherical_to_cartesian_dir(sun.azimuth, sun.zenith, sun_dir); 465 sun_pos[0] = sun_dir[0] * EARTH_TO_SUN_DST; 466 sun_pos[1] = sun_dir[1] * EARTH_TO_SUN_DST; 467 sun_pos[2] = sun_dir[2] * EARTH_TO_SUN_DST; 468 return sun_pos; 469 } 470 471 double /* [W] */ 472 stardis_spherical_source_power(const double time /* [s] */, void* data) 473 { 474 struct boundary_condition* bcond = data; 475 struct scem_sun_pos sun = SCEM_SUN_POS_NULL; 476 double power = 0; /* [W] */ 477 size_t i = 0; /* Index of the meteo entry including "time" */ 478 ASSERT(data); 479 480 compute_sun_position(data, time, &sun); 481 482 i = get_meteo_entry_id(bcond, time); 483 484 /* W/m^2 of surface perpendicular to the incident direction */ 485 power = bcond->lib_desc.smeteo_desc.entries[i].SWdn_direct / sin(sun.zenith); 486 /* Limit to 0 to avoid negative power, i.e., when the sun is below the horizon. */ 487 power = MMAX(power, 0); 488 /* Overall solar power */ 489 power = power * 4*PI*EARTH_TO_SUN_DST*EARTH_TO_SUN_DST; 490 491 return power; 492 } 493 494 double /* [W/m^2/sr] */ 495 stardis_spherical_source_diffuse_radiance 496 (const double time /* [s] */, 497 const double dir[3], 498 void* data) 499 { 500 double diffuse_radiance = 0; /* [W/m^2/sr] */ 501 struct boundary_condition* bcond = data; 502 size_t i = 0; /* Index of the meteo entry including "time" */ 503 504 ASSERT(data); 505 (void)dir; /* Avoid "unused variable" warning */ 506 507 i = get_meteo_entry_id(bcond, time); 508 diffuse_radiance = bcond->lib_desc.smeteo_desc.entries[i].SWdn_diffuse / PI; 509 return diffuse_radiance; 510 } 511 512 const char* 513 get_copyright_notice(void* data) 514 { 515 (void)data; /* Avoid "unused variable" warnings */ 516 return "Copyright (C) 2025 |Méso|Star> (contact@meso-star.com)"; 517 } 518 519 const char* 520 get_license_short(void* data) 521 { 522 (void)data; /* Avoid "unused variable" warnings */ 523 return "GNU GPL version 3 or later <http://www.gnu.org/licenses/>"; 524 } 525 526 const char* 527 get_license_text(void* data) 528 { 529 (void)data; /* Avoid "unused variable" warnings */ 530 return 531 "This is free software released under the GPL v3+ license: GNU GPL\n" 532 "version 3 or later. You are welcome to redistribute it under certain\n" 533 "conditions; refer to <http://www.gnu.org/licenses/> for details."; 534 }