rnatm_radcoef.c (8470B)
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 #define _POSIX_C_SOURCE 200112L /* nextafter */ 24 25 #include "rnatm_c.h" 26 27 #include <star/sars.h> 28 #include <star/sck.h> 29 30 #include <rsys/float4.h> 31 32 /******************************************************************************* 33 * Helper functions 34 ******************************************************************************/ 35 static int 36 tetra_get_radcoef_gas 37 (const struct rnatm* atm, 38 const struct suvm_primitive* tetra, 39 const size_t iband, 40 const size_t iquad_pt, 41 const enum rnatm_radcoef radcoef, 42 float k[4]) 43 { 44 struct sck_band band; 45 float ka[4]; 46 float ks[4]; 47 ASSERT(atm && tetra && k); 48 ASSERT(iband < sck_get_bands_count(atm->gas.ck)); 49 ASSERT(!SUVM_PRIMITIVE_NONE(tetra)); 50 ASSERT(tetra->nvertices == 4); 51 ASSERT(tetra->indices[0] < sck_get_nodes_count(atm->gas.ck)); 52 ASSERT(tetra->indices[1] < sck_get_nodes_count(atm->gas.ck)); 53 ASSERT(tetra->indices[2] < sck_get_nodes_count(atm->gas.ck)); 54 ASSERT(tetra->indices[3] < sck_get_nodes_count(atm->gas.ck)); 55 56 /* Retrieve the band */ 57 SCK(get_band(atm->gas.ck, iband, &band)); 58 ASSERT(iquad_pt < band.quad_pts_count); 59 60 /* Get the absorption coefficients */ 61 if(radcoef == RNATM_RADCOEF_Ka || radcoef == RNATM_RADCOEF_Kext) { 62 /* Retrieve the quadrature point */ 63 struct sck_quad_pt quad; 64 SCK(band_get_quad_pt(&band, iquad_pt, &quad)); 65 ka[0] = quad.ka_list[tetra->indices[0]]; 66 ka[1] = quad.ka_list[tetra->indices[1]]; 67 ka[2] = quad.ka_list[tetra->indices[2]]; 68 ka[3] = quad.ka_list[tetra->indices[3]]; 69 } 70 71 /* Get the scattering coefficients */ 72 if(radcoef == RNATM_RADCOEF_Ks || radcoef == RNATM_RADCOEF_Kext) { 73 ks[0] = band.ks_list[tetra->indices[0]]; 74 ks[1] = band.ks_list[tetra->indices[1]]; 75 ks[2] = band.ks_list[tetra->indices[2]]; 76 ks[3] = band.ks_list[tetra->indices[3]]; 77 } 78 79 switch(radcoef) { 80 case RNATM_RADCOEF_Ka: f4_set(k, ka); break; 81 case RNATM_RADCOEF_Ks: f4_set(k, ks); break; 82 case RNATM_RADCOEF_Kext: f4_add(k, ka, ks); break; 83 default: FATAL("Unreachable code\n"); 84 } 85 86 return 1; 87 } 88 89 static int 90 tetra_get_radcoef_aerosol 91 (const struct rnatm* atm, 92 const struct suvm_primitive* tetra, 93 const size_t iband, 94 const struct aerosol* aerosol, 95 const enum rnatm_radcoef radcoef, 96 float k[4]) 97 { 98 struct sck_band gas_band; 99 struct sars_band ars_band; 100 double gas_spectral_range[2]; /* In nm */ 101 size_t ars_ibands[2]; 102 float ka[4], ks[4]; 103 ASSERT(atm && aerosol && tetra && k); 104 105 ASSERT(!SUVM_PRIMITIVE_NONE(tetra)); 106 ASSERT(tetra->nvertices == 4); 107 ASSERT(tetra->indices[0] < sars_get_nodes_count(aerosol->sars)); 108 ASSERT(tetra->indices[1] < sars_get_nodes_count(aerosol->sars)); 109 ASSERT(tetra->indices[2] < sars_get_nodes_count(aerosol->sars)); 110 ASSERT(tetra->indices[3] < sars_get_nodes_count(aerosol->sars)); 111 112 /* Look for the aerosol bands covered by the gas band */ 113 SCK(get_band(atm->gas.ck, iband, &gas_band)); 114 gas_spectral_range[0] = gas_band.lower; 115 gas_spectral_range[1] = nextafter(gas_band.upper, 0); /* inclusive */ 116 SARS(find_bands(aerosol->sars, gas_spectral_range, ars_ibands)); 117 118 /* No aerosol band is overlaid by the gas band */ 119 if(ars_ibands[0] > ars_ibands[1]) { 120 k[0] = 0; 121 k[1] = 0; 122 k[2] = 0; 123 k[3] = 0; 124 return 0; 125 } 126 127 SARS(get_band(aerosol->sars, ars_ibands[0], &ars_band)); 128 129 /* The gas band is included in an aerosol band: use the aerosol radiative 130 * coefficients */ 131 if(ars_ibands[0] == ars_ibands[1] 132 && ars_band.lower <= gas_band.lower 133 && ars_band.upper >= gas_band.upper) { 134 135 /* Get the absorption coefficient */ 136 if(radcoef == RNATM_RADCOEF_Ka || radcoef == RNATM_RADCOEF_Kext) { 137 ka[0] = sars_band_get_ka(&ars_band, tetra->indices[0]); 138 ka[1] = sars_band_get_ka(&ars_band, tetra->indices[1]); 139 ka[2] = sars_band_get_ka(&ars_band, tetra->indices[2]); 140 ka[3] = sars_band_get_ka(&ars_band, tetra->indices[3]); 141 } 142 143 /* Get the scattering coefficient */ 144 if(radcoef == RNATM_RADCOEF_Ks || radcoef == RNATM_RADCOEF_Kext) { 145 ks[0] = sars_band_get_ks(&ars_band, tetra->indices[0]); 146 ks[1] = sars_band_get_ks(&ars_band, tetra->indices[1]); 147 ks[2] = sars_band_get_ks(&ars_band, tetra->indices[2]); 148 ks[3] = sars_band_get_ks(&ars_band, tetra->indices[3]); 149 } 150 151 /* The gas band overlaid N aerosol bands (N >= 1) */ 152 } else { 153 float tau_ka[4] = {0, 0, 0, 0}; 154 float tau_ks[4] = {0, 0, 0, 0}; 155 double lambda_min; 156 double lambda_max; 157 float rcp_gas_band_len; 158 size_t iars_band; 159 160 FOR_EACH(iars_band, ars_ibands[0], ars_ibands[1]+1) { 161 float lambda_len; 162 163 SARS(get_band(aerosol->sars, iars_band, &ars_band)); 164 lambda_min = MMAX(gas_band.lower, ars_band.lower); 165 lambda_max = MMIN(gas_band.upper, ars_band.upper); /* exclusive */ 166 lambda_max = nextafter(lambda_max, 0); /* inclusive */ 167 lambda_len = (float)(lambda_max - lambda_min); 168 ASSERT(lambda_len > 0); 169 170 /* Get the absorption coefficient */ 171 if(radcoef == RNATM_RADCOEF_Ka || radcoef == RNATM_RADCOEF_Kext) { 172 tau_ka[0] += sars_band_get_ka(&ars_band, tetra->indices[0]) * lambda_len; 173 tau_ka[1] += sars_band_get_ka(&ars_band, tetra->indices[1]) * lambda_len; 174 tau_ka[2] += sars_band_get_ka(&ars_band, tetra->indices[2]) * lambda_len; 175 tau_ka[3] += sars_band_get_ka(&ars_band, tetra->indices[3]) * lambda_len; 176 } 177 178 /* Get the scattering coefficient */ 179 if(radcoef == RNATM_RADCOEF_Ks || radcoef == RNATM_RADCOEF_Kext) { 180 tau_ks[0] += sars_band_get_ks(&ars_band, tetra->indices[0]) * lambda_len; 181 tau_ks[1] += sars_band_get_ks(&ars_band, tetra->indices[1]) * lambda_len; 182 tau_ks[2] += sars_band_get_ks(&ars_band, tetra->indices[2]) * lambda_len; 183 tau_ks[3] += sars_band_get_ks(&ars_band, tetra->indices[3]) * lambda_len; 184 } 185 } 186 187 /* Compute the radiative coefficients of the tetrahedron */ 188 lambda_min = gas_band.lower; 189 lambda_max = nextafter(gas_band.upper, 0); /* inclusive */ 190 rcp_gas_band_len = 1.f/(float)(lambda_max - lambda_min); 191 f4_mulf(ks, tau_ks, rcp_gas_band_len); 192 f4_mulf(ka, tau_ka, rcp_gas_band_len); 193 } 194 195 switch(radcoef) { 196 case RNATM_RADCOEF_Ka: f4_set(k, ka); break; 197 case RNATM_RADCOEF_Ks: f4_set(k, ks); break; 198 case RNATM_RADCOEF_Kext: f4_add(k, ka, ks); break; 199 default: FATAL("Unreachable code\n"); 200 } 201 202 return 1; 203 } 204 205 /******************************************************************************* 206 * Local functions 207 ******************************************************************************/ 208 int 209 tetra_get_radcoef 210 (const struct rnatm* atm, 211 const struct suvm_primitive* tetra, 212 const size_t component, 213 const size_t iband, 214 const size_t iquad_pt, 215 const enum rnatm_radcoef radcoef, 216 float k[4]) 217 { 218 ASSERT(atm); 219 220 if(SUVM_PRIMITIVE_NONE(tetra)) { 221 k[0] = 0; 222 k[1] = 0; 223 k[2] = 0; 224 k[3] = 0; 225 return 0; 226 } 227 228 if(component == RNATM_GAS) { 229 return tetra_get_radcoef_gas(atm, tetra, iband, iquad_pt, radcoef, k); 230 } else { 231 const struct aerosol* aerosol = NULL; 232 ASSERT(component < darray_aerosol_size_get(&atm->aerosols)); 233 aerosol = darray_aerosol_cdata_get(&atm->aerosols) + component; 234 return tetra_get_radcoef_aerosol(atm, tetra, iband, aerosol, radcoef, k); 235 } 236 }