atrstm_radcoefs.h (6642B)
1 /* Copyright (C) 2022, 2023 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2020, 2021 Centre National de la Recherche Scientifique 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 #ifndef ATRSTM_RADCOEFS_H 18 #define ATRSTM_RADCOEFS_H 19 20 #include "atrstm.h" 21 #include "atrstm_rdgfa.h" 22 23 #include <rsys/math.h> 24 #include <rsys/rsys.h> 25 26 /* Radiative coefficients in m^-1 */ 27 struct radcoefs { 28 double ka; /* Absorption coefficient */ 29 double ks; /* Scattering coefficient */ 30 double kext; /* Extinction coefficient */ 31 }; 32 #define RADCOEFS_NULL__ {0,0,0} 33 static const struct radcoefs RADCOEFS_NULL = RADCOEFS_NULL__; 34 35 struct radcoefs_compute_args { 36 double lambda; /* wavelength [nm] */ 37 double n; /* Refractive index (real component) */ 38 double kappa; /* Refractive index (imaginary component) */ 39 double fractal_prefactor; 40 double fractal_dimension; 41 double soot_volumic_fraction; /* In m^3(soot) / m^3 */ 42 double soot_primary_particles_count; 43 double soot_primary_particles_diameter; /* In nm */ 44 int radcoefs_mask; /* Combination of atrstm_radcoef_flag */ 45 }; 46 #define RADCOEFS_COMPUTE_ARGS_NULL__ {0,0,0,0,0,0,0,0,0} 47 static const struct radcoefs_compute_args RADCOEFS_COMPUTE_ARGS_NULL = 48 RADCOEFS_COMPUTE_ARGS_NULL__; 49 50 /* Forward declarations */ 51 struct atrri_refractive_index; 52 struct atrstm; 53 struct suvm_primitive; 54 55 extern LOCAL_SYM int 56 check_fetch_radcoefs_args 57 (const struct atrstm* atrstm, 58 const struct atrstm_fetch_radcoefs_args* args, 59 const char* func_name); 60 61 extern LOCAL_SYM res_T 62 fetch_radcoefs 63 (const struct atrstm* atrstm, 64 const struct atrstm_fetch_radcoefs_args* args, 65 atrstm_radcoefs_T radcoefs); 66 67 /* Write per node radiative coefficients regarding the submitted refracted 68 * index */ 69 extern LOCAL_SYM res_T 70 dump_radcoefs 71 (const struct atrstm* atrstm, 72 const struct atrri_refractive_index* refract_id, 73 FILE* stream); 74 75 extern LOCAL_SYM res_T 76 primitive_compute_radcoefs 77 (const struct atrstm* atrstm, 78 const struct atrri_refractive_index* refract_id, 79 const struct suvm_primitive* prim, 80 const int radcoefs_mask, 81 struct radcoefs radcoefs[]); /* Per primitive node radiative coefficients */ 82 83 static INLINE res_T 84 primitive_compute_radcoefs_range 85 (const struct atrstm* atrstm, 86 const struct atrri_refractive_index* refract_id, 87 const struct suvm_primitive* prim, 88 struct radcoefs* radcoefs_min, 89 struct radcoefs* radcoefs_max) 90 { 91 struct radcoefs props[4]; 92 res_T res = RES_OK; 93 ASSERT(radcoefs_min && radcoefs_max); 94 95 res = primitive_compute_radcoefs 96 (atrstm, refract_id, prim, ATRSTM_RADCOEFS_MASK_ALL, props); 97 if(res != RES_OK) return res; 98 99 radcoefs_min->ka = MMIN 100 (MMIN(props[0].ka, props[1].ka), 101 MMIN(props[2].ka, props[3].ka)); 102 radcoefs_min->ks = MMIN 103 (MMIN(props[0].ks, props[1].ks), 104 MMIN(props[2].ks, props[3].ks)); 105 radcoefs_min->kext = MMIN 106 (MMIN(props[0].kext, props[1].kext), 107 MMIN(props[2].kext, props[3].kext)); 108 109 radcoefs_max->ka = MMAX 110 (MMAX(props[0].ka, props[1].ka), 111 MMAX(props[2].ka, props[3].ka)); 112 radcoefs_max->ks = MMAX 113 (MMAX(props[0].ks, props[1].ks), 114 MMAX(props[2].ks, props[3].ks)); 115 radcoefs_max->kext = MMAX 116 (MMAX(props[0].kext, props[1].kext), 117 MMAX(props[2].kext, props[3].kext)); 118 119 return RES_OK; 120 } 121 122 /* Compute the radiative coefficients of the semi transparent medium as described 123 * in "Effects of multiple scattering on radiative properties of soot fractal 124 * aggregates" J. Yon et al., Journal of Quantitative Spectroscopy and 125 * Radiative Transfer Vol 133, pp. 374-381, 2014 */ 126 static INLINE void 127 radcoefs_compute 128 (struct radcoefs* radcoefs, 129 const struct radcoefs_compute_args* args) 130 { 131 /* Input variables */ 132 const double lambda = args->lambda; /* [nm] */ 133 const double n = args->n; 134 const double kappa = args->kappa; 135 const double Fv = args->soot_volumic_fraction; /* [nm^3(soot)/nm] */ 136 const double Np = args->soot_primary_particles_count; 137 const double Dp = args->soot_primary_particles_diameter; /* [nm] */ 138 const double kf = args->fractal_prefactor; 139 const double Df = args->fractal_dimension; 140 const int ka = (args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ka) != 0; 141 const int ks = (args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ks) != 0; 142 const int kext = (args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Kext) != 0; 143 144 /* Clean up radiative coefficients */ 145 radcoefs->ka = 0; 146 radcoefs->ks = 0; 147 radcoefs->kext = 0; 148 149 if(Np != 0 && Fv != 0 && (args->radcoefs_mask & ATRSTM_RADCOEFS_MASK_ALL)) { 150 /* Precomputed values */ 151 const double n2 = n*n; 152 const double kappa2 = kappa*kappa; 153 const double four_n2_kappa2 = 4*n2*kappa2; 154 const double xp = (PI*Dp) / lambda; 155 const double xp3 = xp*xp*xp; 156 const double k = (2*PI) / lambda; 157 const double k2 = k*k; 158 159 const double Va = (Np*PI*Dp*Dp*Dp) / 6; /* [nm^3] */ 160 const double rho = Fv / Va; /* [#aggregates / nm^3] */ 161 const double tmp0 = 1.e9 * rho; 162 163 /* E(m), m = n + i*kappa */ 164 const double tmp1 = (n2 - kappa2 + 2); 165 const double denom = tmp1*tmp1 + four_n2_kappa2; 166 const double Em = (6*n*kappa) / denom; 167 168 if(ka || kext) { 169 /* Absorption cross section, [nm^3/aggrefate] */ 170 const double Cabs = Np * (4*PI*xp3)/(k2) * Em; 171 172 /* Absorption coefficient, [m^-1] */ 173 radcoefs->ka = tmp0 * Cabs; 174 } 175 176 if(ks || kext) { 177 /* F(m), m = n + i*kappa */ 178 const double tmp2 = ((n2-kappa2 - 1)*(n2-kappa2+2) + four_n2_kappa2)/denom; 179 const double Fm = tmp2*tmp2 + Em*Em; 180 181 /* Gyration radius */ 182 const double Rg = compute_gyration_radius(kf, Df, Np, Dp); /* [nm] */ 183 const double Rg2 = Rg*Rg; 184 const double g = pow(1 + 4*k2*Rg2/(3*Df), -Df*0.5); 185 186 /* Scattering cross section, [nm^3/aggrefate] */ 187 const double Csca = Np*Np * (8*PI*xp3*xp3) / (3*k2) * Fm * g; 188 189 /* Scattering coefficient, [m^-1] */ 190 radcoefs->ks = tmp0 * Csca; 191 } 192 193 if(kext) { 194 radcoefs->kext = radcoefs->ka + radcoefs->ks; 195 } 196 ASSERT(denom != 0 && Df != 0 && Dp != 0); 197 } 198 } 199 200 #endif /* ATRSTM_RADCOEFS_H */