htrdr

Solving radiative transfer in heterogeneous media
git clone git://git.meso-star.fr/htrdr.git
Log | Files | Refs | README | LICENSE

htrdr_planets.c (24333B)


      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 #define _POSIX_C_SOURCE 200112L /* fdopen, nanosleep, nextafter, rint */
     25 
     26 #include "core/htrdr.h"
     27 #include "core/htrdr_ran_wlen_cie_xyz.h"
     28 #include "core/htrdr_ran_wlen_discrete.h"
     29 #include "core/htrdr_ran_wlen_planck.h"
     30 #include "core/htrdr_log.h"
     31 
     32 #include "planets/htrdr_planets.h"
     33 #include "planets/htrdr_planets_args.h"
     34 #include "planets/htrdr_planets_c.h"
     35 #include "planets/htrdr_planets_source.h"
     36 
     37 #include <rad-net/rnatm.h>
     38 #include <rad-net/rngrd.h>
     39 
     40 #include <star/scam.h>
     41 #include <star/smsh.h>
     42 
     43 #include <rsys/cstr.h>
     44 #include <rsys/double3.h>
     45 #include <rsys/mem_allocator.h>
     46 
     47 #include <fcntl.h> /* open */
     48 #include <math.h> /* nextafter, rint */
     49 #include <unistd.h> /* close */
     50 #include <sys/stat.h>
     51 
     52 #include <mpi.h>
     53 #include <time.h>
     54 
     55 /*******************************************************************************
     56  * Helper function
     57  ******************************************************************************/
     58 /* Calculate the number of fixed size spectral intervals to use for the
     59  * cumulative */
     60 static size_t
     61 compute_nintervals_for_spectral_cdf(const struct htrdr_planets* cmd)
     62 {
     63   double range_size;
     64   size_t nintervals;
     65   ASSERT(cmd);
     66 
     67   range_size =
     68     cmd->spectral_domain.wlen_range[1]
     69   - cmd->spectral_domain.wlen_range[0];
     70 
     71   /* Initially assume ~one interval per nanometer */
     72   nintervals = (size_t)rint(range_size);
     73 
     74   return nintervals;
     75 }
     76 
     77 static res_T
     78 setup_octree_storage
     79   (struct htrdr_planets* cmd,
     80    const struct htrdr_planets_args* args,
     81    struct rnatm_create_args* rnatm_args)
     82 {
     83   struct stat file_stat;
     84   int fd = -1;
     85   int err = 0;
     86   res_T res = RES_OK;
     87   ASSERT(cmd && args && rnatm_args);
     88 
     89   rnatm_args->octrees_storage = NULL;
     90   rnatm_args->load_octrees_from_storage = 0;
     91 
     92   if(!args->accel_struct.storage) goto exit;
     93 
     94   fd = open(args->accel_struct.storage, O_CREAT|O_RDWR, S_IRUSR|S_IWUSR);
     95   if(fd < 0) { res = RES_IO_ERR; goto error; }
     96 
     97   rnatm_args->octrees_storage = fdopen(fd, "w+");
     98   if(!rnatm_args->octrees_storage) { res = RES_IO_ERR; goto error; }
     99 
    100   /* From now on, manage the opened file from its pointer and not from its
    101    * descriptor */
    102   fd = -1;
    103 
    104   err = stat(args->accel_struct.storage, &file_stat);
    105   if(err < 0) { res = RES_IO_ERR; goto error; }
    106 
    107   if(file_stat.st_size != 0) {
    108     /* The file is not empty and therefore must contain valid octrees */
    109     rnatm_args->load_octrees_from_storage = 1;
    110   }
    111 
    112 exit:
    113   cmd->octrees_storage = rnatm_args->octrees_storage;
    114   return res;
    115 error:
    116   htrdr_log_err(cmd->htrdr, "error opening the octree storage `%s' -- %s\n",
    117     args->accel_struct.storage, res_to_cstr(res));
    118 
    119   if(fd >= 0) CHK(close(fd) == 0);
    120   if(rnatm_args->octrees_storage) CHK(fclose(rnatm_args->octrees_storage) == 0);
    121   rnatm_args->octrees_storage = NULL;
    122   rnatm_args->load_octrees_from_storage = 1;
    123   goto exit;
    124 }
    125 
    126 static void
    127 mpi_barrier(void)
    128 {
    129   struct timespec t;
    130   int complete = 0;
    131   MPI_Request req;
    132 
    133   t.tv_sec = 0;
    134   t.tv_nsec = 10000000; /* 10ms */
    135 
    136   /* Use an asynchronous barrier to avoid active waiting,
    137    * and therefore wasted resources */
    138   MPI_Ibarrier(MPI_COMM_WORLD, &req);
    139 
    140   do {
    141     nanosleep(&t, NULL);
    142     MPI_Test(&req, &complete, MPI_STATUS_IGNORE);
    143   } while(!complete);
    144 }
    145 
    146 static res_T
    147 setup_atmosphere
    148   (struct htrdr_planets* cmd,
    149    const struct htrdr_planets_args* args)
    150 {
    151   struct rnatm_create_args rnatm = RNATM_CREATE_ARGS_DEFAULT;
    152   res_T res = RES_OK;
    153   ASSERT(cmd && args);
    154 
    155   rnatm.gas = args->gas;
    156   rnatm.aerosols = args->aerosols;
    157   rnatm.naerosols = args->naerosols;
    158   rnatm.name = "atmosphere";
    159   rnatm.spectral_range[0] = args->spectral_domain.wlen_range[0];
    160   rnatm.spectral_range[1] = args->spectral_domain.wlen_range[1];
    161   rnatm.optical_thickness = args->accel_struct.optical_thickness;
    162   rnatm.grid_definition_hint = args->accel_struct.definition_hint;
    163   rnatm.precompute_normals = args->precompute_normals;
    164   rnatm.logger = htrdr_get_logger(cmd->htrdr);
    165   rnatm.allocator = htrdr_get_allocator(cmd->htrdr);
    166   rnatm.nthreads = args->accel_struct.nthreads;
    167   rnatm.verbose = args->verbose;
    168 
    169   if(!args->accel_struct.storage || !args->accel_struct.master_only) {
    170     /* From now on, either the octrees have to be built in memory, or all the
    171      * processes have to build them. In all cases, all processes must create
    172      * the atmopshere data structures. */
    173     if((res = setup_octree_storage(cmd, args, &rnatm)) != RES_OK) goto error;
    174     if((res = rnatm_create(&rnatm, &cmd->atmosphere)) != RES_OK) goto error;
    175 
    176   } else {
    177 
    178     const int is_master_process = htrdr_get_mpi_rank(cmd->htrdr) == 0;
    179     /* A storage space is used and only the master process must fill it with the
    180      * octrees it builds */
    181 
    182     if(is_master_process) {
    183       /* Octrees are built only by the master process and stored on disk.
    184        * Note that in reality, octrees may already be constructed and stored in
    185        * the storage provided. In any case, we can treat this special case as
    186        * the general case, and therefore simply consider that in such situation,
    187        * "construction" by the master process is only faster */
    188       if((res = setup_octree_storage(cmd, args, &rnatm)) != RES_OK) goto error;
    189       if((res = rnatm_create(&rnatm, &cmd->atmosphere)) != RES_OK) goto error;
    190     }
    191 
    192     /* Wait for octrees to be built by the master process */
    193     mpi_barrier();
    194 
    195     if(!is_master_process) {
    196       /* Octrees have been built by the master process.
    197        * Non master processes simply load them. */
    198       if((res = setup_octree_storage(cmd, args, &rnatm)) != RES_OK) goto error;
    199       ASSERT(rnatm.load_octrees_from_storage == 1);
    200       if((res = rnatm_create(&rnatm, &cmd->atmosphere)) != RES_OK) goto error;
    201     }
    202   }
    203 
    204 exit:
    205   return res;
    206 error:
    207   if(cmd->atmosphere) {
    208     RNATM(ref_put(cmd->atmosphere));
    209     cmd->atmosphere = NULL;
    210   }
    211   goto exit;
    212 }
    213 
    214 static res_T
    215 setup_ground
    216   (struct htrdr_planets* cmd,
    217    const struct htrdr_planets_args* args)
    218 {
    219   struct rngrd_create_args rngrd_args = RNGRD_CREATE_ARGS_DEFAULT;
    220   res_T res = RES_OK;
    221   ASSERT(cmd && args);
    222 
    223   if(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_OCTREES)
    224     goto exit;
    225 
    226   rngrd_args.smsh_filename = args->ground.smsh_filename;
    227   rngrd_args.props_filename = args->ground.props_filename;
    228   rngrd_args.mtllst_filename = args->ground.mtllst_filename;
    229   rngrd_args.name = args->ground.name;
    230   rngrd_args.logger = htrdr_get_logger(cmd->htrdr);
    231   rngrd_args.allocator = htrdr_get_allocator(cmd->htrdr);
    232   rngrd_args.verbose = args->verbose;
    233 
    234   res = rngrd_create(&rngrd_args, &cmd->ground);
    235   if(res != RES_OK) goto error;
    236 
    237 exit:
    238   return res;
    239 error:
    240   if(cmd->ground) {
    241     RNGRD(ref_put(cmd->ground));
    242     cmd->ground = NULL;
    243   }
    244   goto exit;
    245 }
    246 
    247 static res_T
    248 setup_spectral_domain_sw
    249   (struct htrdr_planets* cmd,
    250    const struct htrdr_planets_args* args)
    251 {
    252   res_T res = RES_OK;
    253   ASSERT(cmd && args);
    254   ASSERT(cmd->spectral_domain.type == HTRDR_SPECTRAL_SW);
    255 
    256   /* Discrete distribution */
    257   if(args->source.rnrl_filename) {
    258     struct htrdr_planets_source_spectrum spectrum;
    259     struct htrdr_ran_wlen_discrete_create_args discrete_args;
    260 
    261     res = htrdr_planets_source_get_spectrum
    262       (cmd->source, cmd->spectral_domain.wlen_range, &spectrum);
    263     if(res != RES_OK) goto error;
    264 
    265     discrete_args.get = htrdr_planets_source_spectrum_at;
    266     discrete_args.nwavelengths = spectrum.size;
    267     discrete_args.context = &spectrum;
    268     res = htrdr_ran_wlen_discrete_create
    269       (cmd->htrdr, &discrete_args, &cmd->discrete);
    270     if(res != RES_OK) goto error;
    271 
    272   /* Planck distribution */
    273   } else {
    274     const size_t nintervals = compute_nintervals_for_spectral_cdf(cmd);
    275 
    276     /* Use the source temperature as the reference temperature of the Planck
    277      * distribution */
    278     res = htrdr_ran_wlen_planck_create(cmd->htrdr,
    279       cmd->spectral_domain.wlen_range, nintervals, args->source.temperature,
    280       &cmd->planck);
    281     if(res != RES_OK) goto error;
    282   }
    283 
    284 exit:
    285   return res;
    286 error:
    287   goto exit;
    288 }
    289 
    290 static INLINE res_T
    291 setup_spectral_domain
    292   (struct htrdr_planets* cmd,
    293    const struct htrdr_planets_args* args)
    294 {
    295   double ground_T_range[2];
    296   size_t nintervals;
    297   res_T res = RES_OK;
    298   ASSERT(cmd && args);
    299 
    300   /* No spectral distribution required to write octrees */
    301   if(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_OCTREES)
    302     goto exit;
    303 
    304   cmd->spectral_domain = args->spectral_domain;
    305 
    306   /* Configure the spectral distribution */
    307   switch(cmd->spectral_domain.type) {
    308 
    309     case HTRDR_SPECTRAL_LW:
    310       res = rngrd_get_temperature_range(cmd->ground, ground_T_range);
    311       if(res != RES_OK) goto error;
    312 
    313       /* Use as the reference temperature of the Planck distribution the
    314        * maximum scene temperature which, in fact, should be the maximum ground
    315        * temperature */
    316       nintervals = compute_nintervals_for_spectral_cdf(cmd);
    317       res = htrdr_ran_wlen_planck_create(cmd->htrdr,
    318         cmd->spectral_domain.wlen_range, nintervals, ground_T_range[1],
    319         &cmd->planck);
    320       if(res != RES_OK) goto error;
    321       break;
    322 
    323     case HTRDR_SPECTRAL_SW:
    324       res = setup_spectral_domain_sw(cmd, args);
    325       if(res != RES_OK) goto error;
    326       break;
    327 
    328     case HTRDR_SPECTRAL_SW_CIE_XYZ:
    329       /* CIE XYZ distribution */
    330       nintervals = compute_nintervals_for_spectral_cdf(cmd);
    331       res = htrdr_ran_wlen_cie_xyz_create(cmd->htrdr,
    332         cmd->spectral_domain.wlen_range, nintervals, &cmd->cie);
    333       if(res != RES_OK) goto error;
    334       break;
    335 
    336     default: FATAL("Unreachable code\n"); break;
    337   }
    338 
    339 exit:
    340   return res;
    341 error:
    342   goto exit;
    343 }
    344 
    345 static res_T
    346 setup_output
    347   (struct htrdr_planets* cmd,
    348    const struct htrdr_planets_args* args)
    349 {
    350   const char* output_name = NULL;
    351   res_T res = RES_OK;
    352   ASSERT(cmd && args);
    353 
    354   /* No output stream on non master processes */
    355   if(htrdr_get_mpi_rank(cmd->htrdr) != 0) {
    356     cmd->output = NULL;
    357     output_name = "<null>";
    358 
    359   /* Write results on stdout */
    360   } else if(!args->output) {
    361     cmd->output = stdout;
    362     output_name = "<stdout>";
    363 
    364   /* Open the output stream */
    365   } else {
    366     res = htrdr_open_output_stream(cmd->htrdr, args->output, 0/*read*/,
    367       args->force_output_overwrite, &cmd->output);
    368     if(res != RES_OK) goto error;
    369     output_name = args->output;
    370   }
    371 
    372   res = str_set(&cmd->output_name, output_name);
    373   if(res != RES_OK) {
    374     htrdr_log_err(cmd->htrdr, "error storing output stream name `%s' -- %s\n",
    375       output_name, res_to_cstr(res));
    376     goto error;
    377   }
    378 
    379   cmd->output_type = args->output_type;
    380 
    381 exit:
    382   return res;
    383 error:
    384   str_clear(&cmd->output_name);
    385   if(cmd->output && cmd->output != stdout) {
    386     CHK(fclose(cmd->output) == 0);
    387     cmd->output = NULL;
    388   }
    389   goto exit;
    390 }
    391 
    392 static INLINE res_T
    393 setup_source
    394   (struct htrdr_planets* cmd,
    395    const struct htrdr_planets_args* args)
    396 {
    397   res_T res = RES_OK;
    398   ASSERT(cmd && args);
    399 
    400   if(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_OCTREES
    401   /* No source in Longwave.
    402    * Check the spectral domain type on the input arguments, because when the
    403    * source is configured, the spectral dimension may not yet be set */
    404   || args->spectral_domain.type == HTRDR_SPECTRAL_LW)
    405     goto exit;
    406 
    407   res = htrdr_planets_source_create(cmd->htrdr, &args->source, &cmd->source);
    408   if(res != RES_OK) goto error;
    409 
    410 exit:
    411   return res;
    412 error:
    413   goto exit;
    414 }
    415 
    416 static res_T
    417 setup_camera_orthographic
    418   (struct htrdr_planets* cmd,
    419    const struct htrdr_planets_args* args)
    420 {
    421   struct scam_orthographic_args cam_args = SCAM_ORTHOGRAPHIC_ARGS_DEFAULT;
    422   res_T res = RES_OK;
    423 
    424   /* Check pre-conditions */
    425   ASSERT(cmd && args);
    426   ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_IMAGE);
    427   ASSERT(args->cam_type == HTRDR_ARGS_CAMERA_ORTHOGRAPHIC);
    428   ASSERT(htrdr_args_camera_orthographic_check(&args->cam_ortho) == RES_OK);
    429   ASSERT(htrdr_args_image_check(&args->image) == RES_OK);
    430 
    431   d3_set(cam_args.position, args->cam_ortho.position);
    432   d3_set(cam_args.target, args->cam_ortho.target);
    433   d3_set(cam_args.up, args->cam_ortho.up);
    434   cam_args.height = args->cam_ortho.height;
    435   cam_args.aspect_ratio =
    436     (double)args->image.definition[0]
    437   / (double)args->image.definition[1];
    438 
    439   res = scam_create_orthographic
    440     (htrdr_get_logger(cmd->htrdr),
    441      htrdr_get_allocator(cmd->htrdr),
    442      htrdr_get_verbosity_level(cmd->htrdr),
    443      &cam_args,
    444      &cmd->camera);
    445   if(res != RES_OK) goto error;
    446 
    447 exit:
    448   return res;
    449 error:
    450   goto exit;
    451 }
    452 
    453 static res_T
    454 setup_camera_perspective
    455   (struct htrdr_planets* cmd,
    456    const struct htrdr_planets_args* args)
    457 {
    458   struct scam_perspective_args cam_args = SCAM_PERSPECTIVE_ARGS_DEFAULT;
    459   res_T res = RES_OK;
    460 
    461   /* Check pre-conditions */
    462   ASSERT(cmd && args);
    463   ASSERT(args->output_type == HTRDR_PLANETS_ARGS_OUTPUT_IMAGE);
    464   ASSERT(args->cam_type == HTRDR_ARGS_CAMERA_PERSPECTIVE);
    465   ASSERT(htrdr_args_camera_perspective_check(&args->cam_persp) == RES_OK);
    466   ASSERT(htrdr_args_image_check(&args->image) == RES_OK);
    467 
    468   d3_set(cam_args.position, args->cam_persp.position);
    469   d3_set(cam_args.target, args->cam_persp.target);
    470   d3_set(cam_args.up, args->cam_persp.up);
    471   cam_args.field_of_view = MDEG2RAD(args->cam_persp.fov_y);
    472   cam_args.lens_radius = args->cam_persp.lens_radius;
    473   cam_args.focal_distance = args->cam_persp.focal_dst;
    474   cam_args.aspect_ratio =
    475     (double)args->image.definition[0]
    476   / (double)args->image.definition[1];
    477 
    478   res = scam_create_perspective
    479     (htrdr_get_logger(cmd->htrdr),
    480      htrdr_get_allocator(cmd->htrdr),
    481      htrdr_get_verbosity_level(cmd->htrdr),
    482      &cam_args,
    483      &cmd->camera);
    484   if(res != RES_OK) goto error;
    485 
    486 exit:
    487   return res;
    488 error:
    489   goto exit;
    490 }
    491 
    492 static res_T
    493 setup_camera
    494   (struct htrdr_planets* cmd,
    495    const struct htrdr_planets_args* args)
    496 {
    497   res_T res = RES_OK;
    498   ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_IMAGE);
    499 
    500   switch(args->cam_type) {
    501     case HTRDR_ARGS_CAMERA_ORTHOGRAPHIC:
    502       res = setup_camera_orthographic(cmd, args);
    503       break;
    504     case HTRDR_ARGS_CAMERA_PERSPECTIVE:
    505       res = setup_camera_perspective(cmd, args);
    506       break;
    507     default: FATAL("Unreachable code.\n"); break;
    508   }
    509 
    510   return res;
    511 }
    512 
    513 static res_T
    514 setup_buffer_image
    515   (struct htrdr_planets* cmd,
    516    const struct htrdr_planets_args* args)
    517 {
    518   struct htrdr_pixel_format pixfmt = HTRDR_PIXEL_FORMAT_NULL;
    519   res_T res = RES_OK;
    520   ASSERT(cmd && args);
    521   ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_IMAGE);
    522 
    523   planets_get_pixel_format(cmd, &pixfmt);
    524 
    525   /* Setup buffer layout */
    526   cmd->buf_layout.width = args->image.definition[0];
    527   cmd->buf_layout.height = args->image.definition[1];
    528   cmd->buf_layout.pitch = args->image.definition[0] * pixfmt.size;
    529   cmd->buf_layout.elmt_size = pixfmt.size;
    530   cmd->buf_layout.alignment = pixfmt.alignment;
    531 
    532   /* Save the number of samples per pixel */
    533   cmd->spp = args->image.spp;
    534 
    535   /* Create the image buffer only on the master process; Image parts rendered
    536    * by other processes are collected there */
    537   if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit;
    538 
    539   res = htrdr_buffer_create(cmd->htrdr, &cmd->buf_layout, &cmd->buf);
    540   if(res != RES_OK) goto error;
    541 
    542 exit:
    543   return res;
    544 error:
    545   if(cmd->buf) { htrdr_buffer_ref_put(cmd->buf); cmd->buf = NULL; }
    546   goto exit;
    547 }
    548 
    549 static res_T
    550 setup_buffer_raw
    551   (struct htrdr_planets* cmd,
    552    const struct htrdr_planets_args* args)
    553 {
    554   struct smsh_desc desc = SMSH_DESC_NULL;
    555   size_t sz = 0; /* Size of a voxel storing volumic radiative budget */
    556   size_t al = 0; /* Alignment of a voxel storing volumic radiative budget */
    557   res_T res = RES_OK;
    558 
    559   ASSERT(cmd && args);
    560   ASSERT(cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_VOLUMIC_RADIATIVE_BUDGET);
    561   ASSERT(cmd->volrad_mesh != NULL); /* The volurad mesh must be defined */
    562 
    563   res = smsh_get_desc(cmd->volrad_mesh, &desc);
    564   if(res != RES_OK) goto error;
    565 
    566   /* Setup buffer layout for volumic radiative budget calculation */
    567   sz = sizeof(struct planets_voxel_radiative_budget);
    568   al = ALIGNOF(struct planets_voxel_radiative_budget);
    569   cmd->buf_layout.width = desc.ncells;
    570   cmd->buf_layout.height = 1;
    571   cmd->buf_layout.pitch = desc.ncells * sz;
    572   cmd->buf_layout.elmt_size = sz;
    573   cmd->buf_layout.alignment = al;
    574 
    575   /* Save the number of samples per tetrahedron */
    576   cmd->spt = args->volrad_budget.spt;
    577 
    578   /* Create the raw buffer only on master process; buffer parts calculated by
    579    * other processes are collected there */
    580   if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit;
    581 
    582   res = htrdr_buffer_create(cmd->htrdr, &cmd->buf_layout, &cmd->buf);
    583   if(res != RES_OK) goto error;
    584 
    585 exit:
    586   return res;
    587 error:
    588   if(cmd->buf) { htrdr_buffer_ref_put(cmd->buf); cmd->buf = NULL; }
    589   goto exit;
    590 }
    591 
    592 static res_T
    593 setup_buffer
    594   (struct htrdr_planets* cmd,
    595    const struct htrdr_planets_args* args)
    596 {
    597   res_T res = RES_OK;
    598   ASSERT(cmd && args);
    599 
    600   switch(cmd->output_type) {
    601     case HTRDR_PLANETS_ARGS_OUTPUT_IMAGE:
    602       res = setup_buffer_image(cmd, args);
    603       break;
    604     case HTRDR_PLANETS_ARGS_OUTPUT_VOLUMIC_RADIATIVE_BUDGET:
    605       res = setup_buffer_raw(cmd, args);
    606       break;
    607     default: /* Nothing to do */ break;
    608   }
    609   if(res != RES_OK) goto error;
    610 
    611 exit:
    612   return res;
    613 error:
    614   goto exit;
    615 }
    616 
    617 static res_T
    618 setup_volrad_budget_mesh
    619   (struct htrdr_planets* cmd,
    620    const struct htrdr_planets_args* args)
    621 {
    622   struct smsh_create_args create_args = SMSH_CREATE_ARGS_DEFAULT;
    623   struct smsh_load_args load_args = SMSH_LOAD_ARGS_NULL;
    624   struct smsh_desc desc = SMSH_DESC_NULL;
    625   res_T res = RES_OK;
    626   ASSERT(cmd && args);
    627 
    628   if(cmd->output_type != HTRDR_PLANETS_ARGS_OUTPUT_VOLUMIC_RADIATIVE_BUDGET)
    629     goto exit;
    630 
    631   /* Store the number of samples per tetrahedron to be used */
    632   cmd->spt = args->volrad_budget.spt;
    633 
    634   create_args.logger = htrdr_get_logger(cmd->htrdr);
    635   create_args.allocator = htrdr_get_allocator(cmd->htrdr);
    636   create_args.verbose = htrdr_get_verbosity_level(cmd->htrdr);
    637   res = smsh_create(&create_args, &cmd->volrad_mesh);
    638   if(res != RES_OK) goto error;
    639 
    640   load_args.path = args->volrad_budget.smsh_filename;
    641   res = smsh_load(cmd->volrad_mesh, &load_args);
    642   if(res != RES_OK) goto error;
    643 
    644   res = smsh_get_desc(cmd->volrad_mesh, &desc);
    645   if(res != RES_OK) goto error;
    646 
    647   /* Check that the loaded mesh is effectively a volume mesh */
    648   if(desc.dnode != 3 || desc.dcell != 4) {
    649     htrdr_log_err(cmd->htrdr,
    650       "%s: the volumic radiative budget calculation "
    651       "expects a 3D tetrahedral mesh "
    652       "(dimension of mesh: %u; dimension of the vertices: %u)\n",
    653       args->volrad_budget.smsh_filename,
    654       desc.dnode,
    655       desc.dcell);
    656     res = RES_BAD_ARG;
    657     goto error;
    658   }
    659 
    660 exit:
    661   return res;
    662 error:
    663   if(cmd->volrad_mesh) {
    664     SMSH(ref_put(cmd->volrad_mesh));
    665     cmd->volrad_mesh = NULL;
    666   }
    667   goto exit;
    668 }
    669 
    670 static INLINE res_T
    671 write_vtk_octrees(const struct htrdr_planets* cmd)
    672 {
    673   size_t octrees_range[2];
    674   res_T res = RES_OK;
    675   ASSERT(cmd);
    676 
    677   /* Nothing to do on non master process */
    678   if(htrdr_get_mpi_rank(cmd->htrdr) != 0) goto exit;
    679 
    680   octrees_range[0] = 0;
    681   octrees_range[1] = rnatm_get_spectral_items_count(cmd->atmosphere) - 1;
    682 
    683   res = rnatm_write_vtk_octrees(cmd->atmosphere, octrees_range, cmd->output);
    684   if(res != RES_OK) goto error;
    685 
    686 exit:
    687   return res;
    688 error:
    689   goto exit;
    690 }
    691 
    692 static void
    693 planets_release(ref_T* ref)
    694 {
    695   struct htrdr_planets* cmd = CONTAINER_OF(ref, struct htrdr_planets, ref);
    696   struct htrdr* htrdr = NULL;
    697   ASSERT(ref);
    698 
    699   if(cmd->atmosphere) RNATM(ref_put(cmd->atmosphere));
    700   if(cmd->ground) RNGRD(ref_put(cmd->ground));
    701   if(cmd->source) htrdr_planets_source_ref_put(cmd->source);
    702   if(cmd->cie) htrdr_ran_wlen_cie_xyz_ref_put(cmd->cie);
    703   if(cmd->discrete) htrdr_ran_wlen_discrete_ref_put(cmd->discrete);
    704   if(cmd->planck) htrdr_ran_wlen_planck_ref_put(cmd->planck);
    705   if(cmd->octrees_storage) CHK(fclose(cmd->octrees_storage) == 0);
    706   if(cmd->output && cmd->output != stdout) CHK(fclose(cmd->output) == 0);
    707   if(cmd->buf) htrdr_buffer_ref_put(cmd->buf);
    708   if(cmd->camera) SCAM(ref_put(cmd->camera));
    709   if(cmd->volrad_mesh) SMSH(ref_put(cmd->volrad_mesh));
    710   str_release(&cmd->output_name);
    711 
    712   htrdr = cmd->htrdr;
    713   MEM_RM(htrdr_get_allocator(htrdr), cmd);
    714   htrdr_ref_put(htrdr);
    715 }
    716 
    717 /*******************************************************************************
    718  * Exported functions
    719  ******************************************************************************/
    720 res_T
    721 htrdr_planets_create
    722   (struct htrdr* htrdr,
    723    const struct htrdr_planets_args* args,
    724    struct htrdr_planets** out_cmd)
    725 {
    726   struct htrdr_planets* cmd = NULL;
    727   res_T res = RES_OK;
    728   ASSERT(htrdr && out_cmd);
    729 
    730   res = htrdr_planets_args_check(args);
    731   if(res != RES_OK) {
    732     htrdr_log_err(htrdr, "Invalid htrdr_planets arguments -- %s\n",
    733       res_to_cstr(res));
    734     goto error;
    735   }
    736 
    737   cmd = MEM_CALLOC(htrdr_get_allocator(htrdr), 1, sizeof(*cmd));
    738   if(!cmd) {
    739     htrdr_log_err(htrdr, "Error allocating htrdr_planets command\n");
    740     res = RES_MEM_ERR;
    741     goto error;
    742   }
    743   ref_init(&cmd->ref);
    744   htrdr_ref_get(htrdr);
    745   cmd->htrdr = htrdr;
    746   str_init(htrdr_get_allocator(htrdr), &cmd->output_name);
    747 
    748   res = setup_output(cmd, args);
    749   if(res != RES_OK) goto error;
    750   res = setup_source(cmd, args);
    751   if(res != RES_OK) goto error;
    752   res = setup_camera(cmd, args);
    753   if(res != RES_OK) goto error;
    754   res = setup_ground(cmd, args);
    755   if(res != RES_OK) goto error;
    756   res = setup_atmosphere(cmd, args);
    757   if(res != RES_OK) goto error;
    758   res = setup_spectral_domain(cmd, args);
    759   if(res != RES_OK) goto error;
    760   res = setup_volrad_budget_mesh(cmd, args);
    761   if(res != RES_OK) goto error;
    762   res = setup_buffer(cmd, args);
    763   if(res != RES_OK) goto error;
    764 
    765 exit:
    766   *out_cmd = cmd;
    767   return res;
    768 error:
    769   if(cmd) {
    770     htrdr_planets_ref_put(cmd);
    771     cmd = NULL;
    772   }
    773   goto exit;
    774 }
    775 
    776 void
    777 htrdr_planets_ref_get(struct htrdr_planets* cmd)
    778 {
    779   ASSERT(cmd);
    780   ref_get(&cmd->ref);
    781 }
    782 
    783 void
    784 htrdr_planets_ref_put(struct htrdr_planets* cmd)
    785 {
    786   ASSERT(cmd);
    787   ref_put(&cmd->ref, planets_release);
    788 }
    789 
    790 res_T
    791 htrdr_planets_run(struct htrdr_planets* cmd)
    792 {
    793   res_T res = RES_OK;
    794   ASSERT(cmd);
    795 
    796   switch(cmd->output_type) {
    797     case HTRDR_PLANETS_ARGS_OUTPUT_IMAGE:
    798       res = planets_draw_map(cmd);
    799       break;
    800     case HTRDR_PLANETS_ARGS_OUTPUT_OCTREES:
    801       res = write_vtk_octrees(cmd);
    802       break;
    803     case HTRDR_PLANETS_ARGS_OUTPUT_VOLUMIC_RADIATIVE_BUDGET:
    804       res = planets_solve_volrad_budget(cmd);
    805       break;
    806     default: FATAL("Unreachable code\n"); break;
    807   }
    808   if(res != RES_OK) goto error;
    809 
    810 exit:
    811   return res;
    812 error:
    813   goto exit;
    814 }
    815 
    816 /*******************************************************************************
    817  * Local function
    818  ******************************************************************************/
    819 void
    820 planets_get_pixel_format
    821   (const struct htrdr_planets* cmd,
    822    struct htrdr_pixel_format* fmt)
    823 {
    824   ASSERT(cmd && fmt && cmd->output_type == HTRDR_PLANETS_ARGS_OUTPUT_IMAGE);
    825   (void)cmd;
    826 
    827   switch(cmd->spectral_domain.type) {
    828     case HTRDR_SPECTRAL_LW:
    829     case HTRDR_SPECTRAL_SW:
    830       fmt->size = sizeof(struct planets_pixel_xwave);
    831       fmt->alignment = ALIGNOF(struct planets_pixel_xwave);
    832       break;
    833     case HTRDR_SPECTRAL_SW_CIE_XYZ:
    834       fmt->size = sizeof(struct planets_pixel_image);
    835       fmt->alignment = ALIGNOF(struct planets_pixel_image);
    836       break;
    837     default: FATAL("Unreachable code\n"); break;
    838   }
    839 }