rnatm

Load and structure data describing an atmosphere
git clone git://git.meso-star.fr/rnatm.git
Log | Files | Refs | README | LICENSE

rnatm.c (21977B)


      1 /* Copyright (C) 2022, 2023, 2025 Centre National de la Recherche Scientifique
      2  * Copyright (C) 2022, 2023, 2025 Institut Pierre-Simon Laplace
      3  * Copyright (C) 2022, 2023, 2025 Institut de Physique du Globe de Paris
      4  * Copyright (C) 2022, 2023, 2025 |Méso|Star> (contact@meso-star.com)
      5  * Copyright (C) 2022, 2023, 2025 Observatoire de Paris
      6  * Copyright (C) 2022, 2023, 2025 Université de Reims Champagne-Ardenne
      7  * Copyright (C) 2022, 2023, 2025 Université de Versaille Saint-Quentin
      8  * Copyright (C) 2022, 2023, 2025 Université Paul Sabatier
      9  *
     10  * This program is free software: you can redistribute it and/or modify
     11  * it under the terms of the GNU General Public License as published by
     12  * the Free Software Foundation, either version 3 of the License, or
     13  * (at your option) any later version.
     14  *
     15  * This program is distributed in the hope that it will be useful,
     16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     18  * GNU General Public License for more details.
     19  *
     20  * You should have received a copy of the GNU General Public License
     21  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     22 
     23 #include "rnatm.h"
     24 #include "rnatm_c.h"
     25 #include "rnatm_log.h"
     26 #include "rnatm_voxel.h"
     27 
     28 #include <rad-net/rnsf.h>
     29 
     30 #include <star/sars.h>
     31 #include <star/sbuf.h>
     32 #include <star/sck.h>
     33 #include <star/ssf.h>
     34 #include <star/suvm.h>
     35 #include <star/svx.h>
     36 
     37 #include <rsys/cstr.h>
     38 #include <rsys/mem_allocator.h>
     39 #include <rsys/mutex.h>
     40 
     41 #include <omp.h>
     42 
     43 /*******************************************************************************
     44  * Helper functions
     45  ******************************************************************************/
     46 static INLINE res_T
     47 check_rnatm_gas_args(const struct rnatm_gas_args* args)
     48 {
     49   if(!args) return RES_BAD_ARG;
     50 
     51   /* Filenames cannot be NULL */
     52   if(!args->smsh_filename
     53   || !args->sck_filename
     54   || !args->temperatures_filename)
     55     return RES_BAD_ARG;
     56 
     57   return RES_OK;
     58 }
     59 
     60 static INLINE res_T
     61 check_rnatm_aerosol_args(const struct rnatm_aerosol_args* args)
     62 {
     63   if(!args) return RES_BAD_ARG;
     64 
     65   /* Filenames cannot be NULL */
     66   if(!args->smsh_filename
     67   || !args->sars_filename
     68   || !args->phase_fn_ids_filename
     69   || !args->phase_fn_lst_filename)
     70     return RES_BAD_ARG;
     71 
     72   return RES_OK;
     73 }
     74 
     75 static res_T
     76 check_rnatm_create_args(const struct rnatm_create_args* args)
     77 {
     78   size_t i;
     79   res_T res = RES_OK;
     80 
     81   /* Invalid args */
     82   if(!args) return RES_BAD_ARG;
     83 
     84   /* Invalid total number of components */
     85   if(args->naerosols + 1/*gas*/ >= RNATM_MAX_COMPONENTS_COUNT)
     86     return RES_BAD_ARG;
     87 
     88   /* Invalid gas */
     89   res = check_rnatm_gas_args(&args->gas);
     90   if(res != RES_OK) return res;
     91 
     92   /* Invalid aerosols */
     93   FOR_EACH(i, 0, args->naerosols) {
     94     res = check_rnatm_aerosol_args(args->aerosols+i);
     95     if(res != RES_OK) return res;
     96   }
     97 
     98   /* Invalid spectral range */
     99   if(args->spectral_range[0] > args->spectral_range[1])
    100     return RES_BAD_ARG;
    101 
    102   /* Invalid requirements for loading octrees */
    103   if(!args->octrees_storage && args->load_octrees_from_storage)
    104     return RES_BAD_ARG;
    105 
    106   /* Check miscalleneous arguments */
    107   if(!args->name
    108   || args->optical_thickness < 0
    109   || !args->grid_definition_hint
    110   || !args->nthreads)
    111     return RES_BAD_ARG;
    112 
    113   return RES_OK;
    114 }
    115 
    116 static res_T
    117 setup_aerosol_name
    118   (struct rnatm* atm,
    119    const struct rnatm_create_args* args,
    120    const size_t iaerosol)
    121 {
    122   struct aerosol* aerosol = NULL;
    123   res_T res = RES_OK;
    124   ASSERT(atm && args);
    125   ASSERT(iaerosol < args->naerosols);
    126   ASSERT(args->naerosols == darray_aerosol_size_get(&atm->aerosols));
    127 
    128   aerosol = darray_aerosol_data_get(&atm->aerosols)+iaerosol;
    129 
    130   /* Use user-defined name */
    131   if(args->aerosols[iaerosol].name) {
    132     res = str_set(&aerosol->name, args->aerosols[iaerosol].name);
    133     if(res != RES_OK) {
    134       log_err(atm, "could not set the name of the aerosol %lu to `%s' -- %s\n",
    135         (unsigned long)iaerosol,
    136         args->aerosols[iaerosol].name,
    137         res_to_cstr(res));
    138       goto error;
    139     }
    140 
    141   /* Use default name */
    142   } else {
    143     res = str_printf(&aerosol->name, "aerosol%lu", (unsigned long)iaerosol);
    144     if(res != RES_OK) {
    145       log_err(atm,
    146         "could not set the name of the aerosol %lu to `aerosol%lu' -- %s\n",
    147         (unsigned long)iaerosol,
    148         (unsigned long)iaerosol,
    149         res_to_cstr(res));
    150       goto error;
    151     }
    152   }
    153 
    154 exit:
    155   return res;
    156 error:
    157   goto exit;
    158 }
    159 
    160 static res_T
    161 create_rnatm
    162   (const struct rnatm_create_args* args,
    163    struct rnatm** out_atm)
    164 {
    165   struct rnatm* atm = NULL;
    166   struct mem_allocator* allocator = NULL;
    167   size_t i = 0;
    168   int nthreads_max = 0;
    169   res_T res = RES_OK;
    170 
    171   if(!out_atm) { res = RES_BAD_ARG; goto error;}
    172   res = check_rnatm_create_args(args);
    173   if(res != RES_OK) goto error;
    174 
    175   allocator = args->allocator ? args->allocator : &mem_default_allocator;
    176   atm = MEM_CALLOC(allocator, 1, sizeof(*atm));
    177   if(!atm) {
    178     if(args->verbose) {
    179       #define ERR_STR \
    180         "could not allocate the device of the Rad-Net ATMosphere library"
    181       if(args->logger) {
    182         logger_print(args->logger, LOG_ERROR, ERR_STR);
    183       } else {
    184         fprintf(stderr, MSG_ERROR_PREFIX ERR_STR);
    185       }
    186       #undef ERR_STR
    187     }
    188     res = RES_MEM_ERR;
    189     goto error;
    190   }
    191   ref_init(&atm->ref);
    192   atm->allocator = allocator;
    193   atm->verbose = args->verbose;
    194   str_init(atm->allocator, &atm->name);
    195   gas_init(atm->allocator, &atm->gas);
    196   darray_aerosol_init(atm->allocator, &atm->aerosols);
    197   darray_accel_struct_init(atm->allocator, &atm->accel_structs);
    198   darray_band_init(atm->allocator, &atm->bands);
    199 
    200   /* Setup the number of threads */
    201   nthreads_max = MMAX(omp_get_max_threads(), omp_get_num_procs());
    202   atm->nthreads = MMIN((unsigned)nthreads_max, args->nthreads);
    203 
    204   if(args->logger) {
    205     atm->logger = args->logger;
    206   } else {
    207     res = setup_log_default(atm);
    208     if(res != RES_OK) {
    209       if(args->verbose) {
    210         fprintf(stderr, MSG_ERROR_PREFIX
    211           "could not setup the default logger of the "
    212           "Rad-Net ATMopshere library\n");
    213       }
    214       goto error;
    215     }
    216   }
    217 
    218   res = darray_aerosol_resize(&atm->aerosols, args->naerosols);
    219   if(res != RES_OK) {
    220     log_err(atm, "could not allocate aerosol list -- %s\n", res_to_cstr(res));
    221     goto error;
    222   }
    223 
    224   res = str_set(&atm->name, args->name);
    225   if(res != RES_OK) {
    226     log_err(atm, "could not setup the atmosphere name to `%s' -- %s\n",
    227       args->name, res_to_cstr(res));
    228     goto error;
    229   }
    230 
    231   FOR_EACH(i, 0, args->naerosols) {
    232     res = setup_aerosol_name(atm, args, i);
    233     if(res != RES_OK) goto error;
    234   }
    235 
    236   res = mem_init_regular_allocator(&atm->svx_allocator);
    237   if(res != RES_OK) {
    238     log_err(atm,
    239       "unable to initialize the allocator used to manage the memory of the "
    240       "Star-VoXel library -- %s\n", res_to_cstr(res));
    241     goto error;
    242   }
    243   atm->svx_allocator_is_init = 1;
    244 
    245   atm->mutex = mutex_create();
    246   if(!atm->mutex) {
    247     log_err(atm, "unable to create Rad-Net ATMopshere library mutex\n");
    248     res = RES_MEM_ERR;
    249     goto error;
    250   }
    251 
    252 exit:
    253   if(out_atm) *out_atm = atm;
    254   return res;
    255 error:
    256   if(atm) { RNATM(ref_put(atm)); atm = NULL; }
    257   goto exit;
    258 }
    259 
    260 static void
    261 release_rnatm(ref_T* ref)
    262 {
    263   struct rnatm* atm = CONTAINER_OF(ref, struct rnatm, ref);
    264   ASSERT(ref);
    265   if(atm->logger == &atm->logger__) logger_release(&atm->logger__);
    266   if(atm->svx) SVX(device_ref_put(atm->svx));
    267   if(atm->mutex) mutex_destroy(atm->mutex);
    268   darray_aerosol_release(&atm->aerosols);
    269   darray_accel_struct_release(&atm->accel_structs);
    270   darray_band_release(&atm->bands);
    271   if(atm->svx_allocator_is_init) {
    272     ASSERT(MEM_ALLOCATED_SIZE(&atm->svx_allocator) == 0);
    273     mem_shutdown_regular_allocator(&atm->svx_allocator);
    274   }
    275   gas_release(&atm->gas);
    276   str_release(&atm->name);
    277   MEM_RM(atm->allocator, atm);
    278 }
    279 
    280 /*******************************************************************************
    281  * Exported symbols
    282  ******************************************************************************/
    283 res_T
    284 rnatm_create
    285   (const struct rnatm_create_args* args,
    286    struct rnatm** out_atm)
    287 {
    288   struct rnatm* atm = NULL;
    289   res_T res = RES_OK;
    290 
    291   res = create_rnatm(args, &atm);
    292   if(res != RES_OK) goto error;
    293 
    294   res = setup_meshes(atm, args);
    295   if(res != RES_OK) goto error;
    296   res = setup_properties(atm, args);
    297   if(res != RES_OK) goto error;
    298   res = setup_octrees(atm, args);
    299   if(res != RES_OK) goto error;
    300 
    301 exit:
    302   if(out_atm) *out_atm = atm;
    303   return res;
    304 error:
    305   if(atm) { RNATM(ref_put(atm)); atm = NULL; }
    306   goto exit;
    307 }
    308 
    309 res_T
    310 rnatm_ref_get(struct rnatm* atm)
    311 {
    312   if(!atm) return RES_BAD_ARG;
    313   ref_get(&atm->ref);
    314   return RES_OK;
    315 }
    316 
    317 res_T
    318 rnatm_ref_put(struct rnatm* atm)
    319 {
    320   if(!atm) return RES_BAD_ARG;
    321   ref_put(&atm->ref, release_rnatm);
    322   return RES_OK;
    323 }
    324 
    325 res_T
    326 rnatm_validate(const struct rnatm* atm)
    327 {
    328   res_T res = RES_OK;
    329 
    330   if(!atm) { res = RES_BAD_ARG; goto error; }
    331 
    332   res = check_properties(atm);
    333   if(res != RES_OK) goto error;
    334 
    335 exit:
    336   return res;
    337 error:
    338   goto exit;
    339 }
    340 
    341 double
    342 rnatm_get_k_svx_voxel
    343   (const struct rnatm* atm,
    344    const struct svx_voxel* voxel,
    345    const enum rnatm_radcoef radcoef,
    346    const enum rnatm_svx_op op)
    347 {
    348   const float* vx;
    349   ASSERT(atm && voxel);
    350   ASSERT((unsigned)radcoef < RNATM_RADCOEFS_COUNT__);
    351   ASSERT((unsigned)op < RNATM_SVX_OPS_COUNT__);
    352   (void)atm;
    353 
    354   vx = voxel->data;
    355 
    356   /* Absorption/Scattering coefficient */
    357   if(radcoef != RNATM_RADCOEF_Kext) {
    358     const float k_min = vx[voxel_idata(radcoef, RNATM_SVX_OP_MIN)];
    359     const float k_max = vx[voxel_idata(radcoef, RNATM_SVX_OP_MAX)];
    360     if(k_min > k_max) return 0; /* Empty voxel => null radiative coefficient */
    361     switch(op) {
    362       case RNATM_SVX_OP_MIN: return k_min;
    363       case RNATM_SVX_OP_MAX: return k_max;
    364       default: FATAL("Unreachable code\n"); break;
    365     }
    366 
    367   /* Extinction coefficient */
    368   } else {
    369     const float ka_min = vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MIN)];
    370     const float ka_max = vx[voxel_idata(RNATM_RADCOEF_Ka, RNATM_SVX_OP_MAX)];
    371     const float ks_min = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MIN)];
    372     const float ks_max = vx[voxel_idata(RNATM_RADCOEF_Ks, RNATM_SVX_OP_MAX)];
    373 
    374     /* Unlike absorption and diffusion coefficients, the range of the
    375      * extinction coefficient is not stored by voxel but is recalculated from
    376      * ka and ks by voxel; k_ext = ka + ks. However, the extinction coefficient
    377      * of a tetrahedron is calculated by first adding ka and ks per component
    378      * before adding them to obtain the k_ext of the mixture. Although the
    379      * results should be identical, the reorganization of sums produces a
    380      * numerical inconsistency. Therefore, the k_ext tetrahedron might not be
    381      * strictly included in the range of k_ext computed for the voxel. To deal
    382      * with this numerical inaccuracy, we therefore use the following epsilon
    383      * to slightly increase the returned range of k_ext */
    384     const float epsilon = 1e-6f;
    385 
    386     /* Empty voxel => null radiative coefficient */
    387     if(ka_min > ka_max) { ASSERT(ks_min > ks_max); return 0; }
    388 
    389     switch(op) {
    390       case RNATM_SVX_OP_MIN: return (ka_min + ks_min)*(1.f-epsilon);
    391       case RNATM_SVX_OP_MAX: return (ka_max + ks_max)*(1.f+epsilon);
    392       default: FATAL("Unreachable code\n"); break;
    393     }
    394   }
    395 }
    396 
    397 size_t
    398 rnatm_get_aerosols_count(const struct rnatm* atm)
    399 {
    400   ASSERT(atm);
    401   return darray_aerosol_size_get(&atm->aerosols);
    402 }
    403 
    404 size_t
    405 rnatm_get_spectral_items_count(const struct rnatm* atm)
    406 {
    407   ASSERT(atm);
    408   return darray_accel_struct_size_get(&atm->accel_structs);
    409 }
    410 
    411 res_T
    412 rnatm_find_bands
    413   (const struct rnatm* rnatm,
    414    const double range[2], /* In nm. Limits are inclusive */
    415    size_t ibands[2]) /* Range of overlaped bands. Limits are inclusive */
    416 {
    417   double range_adjusted[2];
    418 
    419   if(!rnatm
    420   || !range
    421   || range[0] > rnatm->spectral_range[1]
    422   || range[1] < rnatm->spectral_range[0])
    423     return RES_BAD_ARG;
    424 
    425   /* Clamp the submitted range to the spectral domain of the atmosphere */
    426   range_adjusted[0] = MMAX(range[0], rnatm->spectral_range[0]);
    427   range_adjusted[1] = MMIN(range[1], rnatm->spectral_range[1]);
    428 
    429   return sck_find_bands(rnatm->gas.ck, range_adjusted, ibands);
    430 }
    431 
    432 res_T
    433 rnatm_band_sample_quad_pt
    434   (const struct rnatm* rnatm,
    435    const double r, /* Canonical random number in [0, 1[ */
    436    const size_t iband,
    437    size_t* iquad)
    438 {
    439   struct sck_band band;
    440   res_T res = RES_OK;
    441 
    442   if(!rnatm || !iquad) {
    443     res = RES_BAD_ARG;
    444     goto error;
    445   }
    446 
    447   res = sck_get_band(rnatm->gas.ck, iband, &band);
    448   if(res != RES_OK) goto error;
    449 
    450   /* Reject the band if it does not overlap the atmosphere spectral domain */
    451   if(band.lower/*Inclusive*/  > rnatm->spectral_range[1]
    452   || band.upper/*Exclusive*/ <= rnatm->spectral_range[0]) {
    453     return RES_BAD_ARG;
    454   }
    455 
    456   res = sck_band_sample_quad_pt(&band, r, iquad);
    457   if(res != RES_OK) goto error;
    458 
    459 exit:
    460   return res;
    461 error:
    462   goto exit;
    463 }
    464 
    465 res_T
    466 rnatm_band_get_desc
    467   (const struct rnatm* rnatm,
    468    const size_t iband,
    469    struct rnatm_band_desc* band_desc)
    470 {
    471   struct sck_band band;
    472   res_T res = RES_OK;
    473 
    474   if(!rnatm || !band_desc) {
    475     res = RES_BAD_ARG;
    476     goto error;
    477   }
    478 
    479   res = sck_get_band(rnatm->gas.ck, iband, &band);
    480   if(res != RES_OK) goto error;
    481 
    482   /* Reject the band if it does not overlap the atmosphere spectral domain */
    483   if(band.lower/*Inclusive*/  > rnatm->spectral_range[1]
    484   || band.upper/*Exclusive*/ <= rnatm->spectral_range[0]) {
    485     return RES_BAD_ARG;
    486   }
    487 
    488   /* TODO Should the range of the band be clamped to the atmospheric spectral
    489    * domain? */
    490   band_desc->lower = band.lower;
    491   band_desc->upper = band.upper;
    492   band_desc->quad_pts_count = band.quad_pts_count;
    493 
    494 exit:
    495   return res;
    496 error:
    497   goto exit;
    498 }
    499 
    500 /*******************************************************************************
    501  * Local functions
    502  ******************************************************************************/
    503 res_T
    504 phase_init(struct mem_allocator* allocator, struct ssf_phase** phase)
    505 {
    506   ASSERT(phase);
    507   (void)allocator;
    508   *phase = NULL;
    509   return RES_OK;
    510 }
    511 
    512 void
    513 phase_release(struct ssf_phase** phase)
    514 {
    515   ASSERT(phase);
    516   if(*phase) SSF(phase_ref_put(*phase));
    517 }
    518 
    519 res_T
    520 phase_copy(struct ssf_phase** dst, struct ssf_phase* const* src)
    521 {
    522   ASSERT(dst && src);
    523   if(*dst) SSF(phase_ref_put(*dst));
    524   *dst = *src;
    525   if(*dst) SSF(phase_ref_get(*dst));
    526   return RES_OK;
    527 }
    528 
    529 res_T
    530 phase_copy_and_release(struct ssf_phase** dst, struct ssf_phase** src)
    531 {
    532   ASSERT(dst && src);
    533   if(*dst) SSF(phase_ref_put(*dst));
    534   *dst = *src;
    535   *src = NULL;
    536   return RES_OK;
    537 }
    538 
    539 res_T
    540 phase_fn_init(struct mem_allocator* allocator, struct phase_fn* phase_fn)
    541 {
    542   ASSERT(phase_fn);
    543   phase_fn->rnsf = NULL;
    544   darray_phase_init(allocator, &phase_fn->phase_lst);
    545   return RES_OK;
    546 }
    547 
    548 void
    549 phase_fn_release(struct phase_fn* phase_fn)
    550 {
    551   ASSERT(phase_fn);
    552   if(phase_fn->rnsf) RNSF(ref_put(phase_fn->rnsf));
    553   darray_phase_release(&phase_fn->phase_lst);
    554 }
    555 
    556 res_T
    557 phase_fn_copy(struct phase_fn* dst, const struct phase_fn* src)
    558 {
    559   ASSERT(dst && src);
    560   if(dst->rnsf) RNSF(ref_put(dst->rnsf));
    561   dst->rnsf = src->rnsf;
    562   if(dst->rnsf) RNSF(ref_get(dst->rnsf));
    563   return darray_phase_copy(&dst->phase_lst, &src->phase_lst);
    564 }
    565 
    566 res_T
    567 phase_fn_copy_and_release(struct phase_fn* dst, struct phase_fn* src)
    568 {
    569   ASSERT(dst && src);
    570   if(dst->rnsf) RNSF(ref_put(dst->rnsf));
    571   dst->rnsf = src->rnsf;
    572   src->rnsf = NULL;
    573   return darray_phase_copy_and_release(&dst->phase_lst, &src->phase_lst);
    574 }
    575 
    576 res_T
    577 gas_init(struct mem_allocator* allocator, struct gas* gas)
    578 {
    579   (void)allocator;
    580   ASSERT(gas);
    581   gas->volume = NULL;
    582   gas->rayleigh = NULL;
    583   gas->temperatures = NULL;
    584   gas->ck = NULL;
    585   gas->ntetrahedra = 0;
    586   gas->nvertices = 0;
    587   return RES_OK;
    588 }
    589 
    590 void
    591 gas_release(struct gas* gas)
    592 {
    593   ASSERT(gas);
    594   if(gas->volume) SUVM(volume_ref_put(gas->volume));
    595   if(gas->rayleigh) SSF(phase_ref_put(gas->rayleigh));
    596   if(gas->temperatures) SBUF(ref_put(gas->temperatures));
    597   if(gas->ck) SCK(ref_put(gas->ck));
    598 }
    599 
    600 res_T
    601 gas_copy(struct gas* dst, const struct gas* src)
    602 {
    603   ASSERT(dst && src);
    604   if(dst->volume) SUVM(volume_ref_put(dst->volume));
    605   if(dst->rayleigh) SSF(phase_ref_put(dst->rayleigh));
    606   if(dst->temperatures) SBUF(ref_put(dst->temperatures));
    607   if(dst->ck) SCK(ref_put(dst->ck));
    608   dst->volume = src->volume;
    609   dst->rayleigh = src->rayleigh;
    610   dst->temperatures = src->temperatures;
    611   dst->ck = src->ck;
    612   dst->ntetrahedra = src->ntetrahedra;
    613   dst->nvertices = src->nvertices;
    614   if(dst->volume) SUVM(volume_ref_get(dst->volume));
    615   if(dst->rayleigh) SSF(phase_ref_get(dst->rayleigh));
    616   if(dst->temperatures) SBUF(ref_get(dst->temperatures));
    617   if(dst->ck) SCK(ref_get(dst->ck));
    618   return RES_OK;
    619 }
    620 
    621 res_T
    622 gas_copy_and_release(struct gas* dst, struct gas* src)
    623 {
    624   ASSERT(dst && src);
    625   if(dst->volume) SUVM(volume_ref_put(dst->volume));
    626   if(dst->rayleigh) SSF(phase_ref_put(dst->rayleigh));
    627   if(dst->temperatures) SBUF(ref_put(dst->temperatures));
    628   if(dst->ck) SCK(ref_put(dst->ck));
    629   dst->volume = src->volume;
    630   dst->rayleigh = src->rayleigh;
    631   dst->temperatures = src->temperatures;
    632   dst->ck = src->ck;
    633   dst->ntetrahedra = src->ntetrahedra;
    634   dst->nvertices = src->nvertices;
    635   src->volume = NULL;
    636   src->rayleigh = NULL;
    637   src->temperatures = NULL;
    638   src->ck = NULL;
    639   return RES_OK;
    640 }
    641 
    642 res_T
    643 aerosol_init(struct mem_allocator* allocator, struct aerosol* aerosol)
    644 {
    645   (void)allocator;
    646   ASSERT(aerosol);
    647   darray_phase_fn_init(allocator, &aerosol->phase_fn_lst);
    648   str_init(allocator, &aerosol->name);
    649   aerosol->volume = NULL;
    650   aerosol->phase_fn_ids = NULL;
    651   aerosol->sars = NULL;
    652   aerosol->ntetrahedra = 0;
    653   aerosol->nvertices = 0;
    654   return RES_OK;
    655 }
    656 
    657 void
    658 aerosol_release(struct aerosol* aerosol)
    659 {
    660   ASSERT(aerosol);
    661   darray_phase_fn_release(&aerosol->phase_fn_lst);
    662   str_release(&aerosol->name);
    663   if(aerosol->volume) SUVM(volume_ref_put(aerosol->volume));
    664   if(aerosol->phase_fn_ids) SBUF(ref_put(aerosol->phase_fn_ids));
    665   if(aerosol->sars) SARS(ref_put(aerosol->sars));
    666 }
    667 
    668 res_T
    669 aerosol_copy(struct aerosol* dst, const struct aerosol* src)
    670 {
    671   res_T res = RES_OK;
    672   ASSERT(dst && src);
    673 
    674   if(dst->volume) SUVM(volume_ref_put(dst->volume));
    675   if(dst->phase_fn_ids) SBUF(ref_put(dst->phase_fn_ids));
    676   if(dst->sars) SARS(ref_put(dst->sars));
    677 
    678   dst->volume = src->volume;
    679   dst->phase_fn_ids = src->phase_fn_ids;
    680   dst->sars = src->sars;
    681   dst->ntetrahedra = src->ntetrahedra;
    682   dst->nvertices = src->nvertices;
    683 
    684   if(dst->volume) SUVM(volume_ref_get(dst->volume));
    685   if(dst->phase_fn_ids) SBUF(ref_get(dst->phase_fn_ids));
    686   if(dst->sars) SARS(ref_get(dst->sars));
    687 
    688   res = darray_phase_fn_copy(&dst->phase_fn_lst, &src->phase_fn_lst);
    689   if(res != RES_OK) return res;
    690   res = str_copy(&dst->name, &src->name);
    691   if(res != RES_OK) return res;
    692 
    693   return RES_OK;
    694 }
    695 
    696 res_T
    697 aerosol_copy_and_release(struct aerosol* dst, struct aerosol* src)
    698 {
    699   res_T res = RES_OK;
    700   ASSERT(dst && src);
    701 
    702   if(dst->volume) SUVM(volume_ref_put(dst->volume));
    703   if(dst->phase_fn_ids) SBUF(ref_put(dst->phase_fn_ids));
    704   if(dst->sars) SARS(ref_put(dst->sars));
    705   dst->volume = src->volume;
    706   dst->phase_fn_ids = src->phase_fn_ids;
    707   dst->sars = src->sars;
    708   dst->ntetrahedra = src->ntetrahedra;
    709   dst->nvertices = src->nvertices;
    710   src->volume = NULL;
    711   src->phase_fn_ids = NULL;
    712   src->sars = NULL;
    713 
    714   res = darray_phase_fn_copy_and_release(&dst->phase_fn_lst, &src->phase_fn_lst);
    715   if(res != RES_OK) return res;
    716   res = str_copy_and_release(&dst->name, &src->name);
    717   if(res != RES_OK) return res;
    718 
    719   return RES_OK;
    720 }
    721 
    722 res_T
    723 accel_struct_init
    724   (struct mem_allocator* allocator,
    725    struct accel_struct* accel_struct)
    726 {
    727   (void)allocator;
    728   ASSERT(accel_struct);
    729   accel_struct->octree = NULL;
    730   accel_struct->iband = 0;
    731   accel_struct->iquad_pt = 0;
    732   memset(&accel_struct->fpos, 0, sizeof(accel_struct->fpos));
    733   return RES_OK;
    734 }
    735 
    736 void
    737 accel_struct_release(struct accel_struct* accel_struct)
    738 {
    739   ASSERT(accel_struct);
    740   if(accel_struct->octree) SVX(tree_ref_put(accel_struct->octree));
    741 }
    742 
    743 res_T
    744 accel_struct_copy(struct accel_struct* dst, const struct accel_struct* src)
    745 {
    746   ASSERT(dst && src);
    747   if(dst->octree) SVX(tree_ref_put(dst->octree));
    748   dst->octree = src->octree;
    749   dst->fpos = src->fpos;
    750   dst->iband = src->iband;
    751   dst->iquad_pt = src->iquad_pt;
    752   if(dst->octree) SVX(tree_ref_get(dst->octree));
    753   return RES_OK;
    754 }
    755 
    756 res_T
    757 accel_struct_copy_and_release(struct accel_struct* dst, struct accel_struct* src)
    758 {
    759   ASSERT(dst && src);
    760   if(dst->octree) SVX(tree_ref_put(dst->octree));
    761   dst->octree = src->octree;
    762   dst->fpos = src->fpos;
    763   dst->iband = src->iband;
    764   dst->iquad_pt = src->iquad_pt;
    765   src->octree = NULL;
    766   return RES_OK;
    767 }
    768 
    769 res_T
    770 make_sure_octree_is_loaded(struct rnatm* atm, const size_t istruct)
    771 {
    772   struct accel_struct* accel_struct = NULL;
    773   int err = 0;
    774   res_T res = RES_OK;
    775   ASSERT(atm && istruct < darray_accel_struct_size_get(&atm->accel_structs));
    776 
    777   accel_struct = darray_accel_struct_data_get(&atm->accel_structs) + istruct;
    778 
    779   mutex_lock(atm->mutex);
    780 
    781   if(accel_struct->octree) {
    782     mutex_unlock(atm->mutex);
    783     goto exit; /* The octree is already loaded */
    784   }
    785 
    786   /* Prepare to read the octree's data */
    787   err = fsetpos(atm->octrees_storage, &accel_struct->fpos);
    788   if(err != 0) {
    789     res = RES_IO_ERR;
    790     mutex_unlock(atm->mutex);
    791     goto error;
    792   }
    793 
    794   /* Deserialize the octree */
    795   res = svx_tree_create_from_stream
    796     (atm->svx, atm->octrees_storage, &accel_struct->octree);
    797   if(res != RES_OK) {
    798     mutex_unlock(atm->mutex);
    799     goto error;
    800   }
    801 
    802   mutex_unlock(atm->mutex);
    803 
    804 exit:
    805   return res;
    806 error:
    807   log_err(atm, "error loading octree %lu -- %s",
    808     (unsigned long)istruct, res_to_cstr(res));
    809   goto exit;
    810 }
    811 
    812 void
    813 unload_octree(struct rnatm* atm, const size_t istruct)
    814 {
    815   struct accel_struct* accel_struct = NULL;
    816   ASSERT(atm);
    817   ASSERT(istruct < darray_accel_struct_size_get(&atm->accel_structs));
    818 
    819   accel_struct = darray_accel_struct_data_get(&atm->accel_structs) + istruct;
    820   if(accel_struct->octree) {
    821     SVX(tree_ref_put(accel_struct->octree));
    822     accel_struct->octree = NULL;
    823   }
    824 }