htrdr_spectral.h (5432B)
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 #ifndef HTRDR_SPECTRAL_H 25 #define HTRDR_SPECTRAL_H 26 27 #include "core/htrdr.h" 28 29 #include <rsys/rsys.h> 30 #include <rsys/math.h> /* PI support */ 31 32 #define HTRDR_SUN_TEMPERATURE 5778 /* In K */ 33 #define HTRDR_DEFAULT_LW_REF_TEMPERATURE 290 /* Default longwave temperature in K */ 34 35 enum htrdr_spectral_type { 36 HTRDR_SPECTRAL_LW, /* Longwave */ 37 HTRDR_SPECTRAL_SW, /* Shortwave */ 38 HTRDR_SPECTRAL_SW_CIE_XYZ /* Shortwave wrt the CIE XYZ tristimulus */ 39 }; 40 41 /* Forwar declaration */ 42 struct htrdr; 43 44 static FINLINE double /* In nanometer */ 45 htrdr_wavenumber_to_wavelength(const double nu/*In cm⁻¹*/) 46 { 47 return 1.e7 / nu; 48 } 49 50 static FINLINE double /* In cm⁻¹ */ 51 wavelength_to_wavenumber(const double lambda/*In nanometer*/) 52 { 53 return htrdr_wavenumber_to_wavelength(lambda); 54 } 55 56 static INLINE double 57 htrdr_wiebelt(const double v) 58 { 59 int m; 60 double w, v2, v4; 61 /*.153989717364e+00;*/ 62 const double fifteen_over_pi_power_4 = 15.0/(PI*PI*PI*PI); 63 const double z0 = 1.0/3.0; 64 const double z1 = 1.0/8.0; 65 const double z2 = 1.0/60.0; 66 const double z4 = 1.0/5040.0; 67 const double z6 = 1.0/272160.0; 68 const double z8 = 1.0/13305600.0; 69 70 if(v >= 2.) { 71 w = 0; 72 for(m=1; m<6 ;m++) 73 w+=exp(-m*v)/(m*m*m*m) * (((m*v+3)*m*v+6)*m*v+6); 74 w = w * fifteen_over_pi_power_4; 75 } else { 76 v2 = v*v; 77 v4 = v2*v2; 78 w = z0 - z1*v + z2*v2 - z4*v2*v2 + z6*v4*v2 - z8*v4*v4; 79 w = 1. - fifteen_over_pi_power_4*v2*v*w; 80 } 81 ASSERT(w >= 0.0 && w <= 1.0); 82 return w; 83 } 84 85 static INLINE double 86 htrdr_blackbody_fraction 87 (const double lambda0, /* In meters */ 88 const double lambda1, /* In meters */ 89 const double temperature) /* In Kelvin */ 90 { 91 const double C2 = 1.43877735e-2; /* m.K */ 92 double x0 = C2 / lambda0; 93 double x1 = C2 / lambda1; 94 double v0 = x0 / temperature; 95 double v1 = x1 / temperature; 96 double w0 = htrdr_wiebelt(v0); 97 double w1 = htrdr_wiebelt(v1); 98 return w1 - w0; 99 } 100 101 /* Return the Planck value in W/m²/sr/m at a given wavelength */ 102 static INLINE double 103 htrdr_planck_monochromatic 104 (const double lambda, /* In meters */ 105 const double temperature) /* In Kelvin */ 106 { 107 const double c = 299792458; /* m/s */ 108 const double h = 6.62607015e-34; /* J.s */ 109 const double k = 1.380649e-23; /* J/K */ 110 const double lambda2 = lambda*lambda; 111 const double lambda5 = lambda2*lambda2*lambda; 112 const double B = ((2.0 * h * c*c) / lambda5) /* W/m²/sr/m */ 113 / (exp(h*c/(lambda*k*temperature))-1.0); 114 ASSERT(temperature > 0); 115 return B; 116 } 117 118 /* Return the average Planck value in W/m²/sr/m over the 119 * [lambda_min, lambda_max] interval. */ 120 static INLINE double 121 htrdr_planck_interval 122 (const double lambda_min, /* In meters */ 123 const double lambda_max, /* In meters */ 124 const double temperature) /* In Kelvin */ 125 { 126 const double T2 = temperature*temperature; 127 const double T4 = T2*T2; 128 const double BOLTZMANN_CONSTANT = 5.6696e-8; /* W/m²/K⁴ */ 129 ASSERT(lambda_min < lambda_max && temperature > 0); 130 return htrdr_blackbody_fraction(lambda_min, lambda_max, temperature) 131 * BOLTZMANN_CONSTANT * T4 / (PI * (lambda_max-lambda_min)); /* In W/m²/sr/m */ 132 } 133 134 /* Invoke planck_monochromatic or planck_interval whether the submitted 135 * interval is null or not, respectively. The returned value is in W/m²/sr/m */ 136 static FINLINE double 137 htrdr_planck 138 (const double lambda_min, /* In meters */ 139 const double lambda_max, /* In meters */ 140 const double temperature) /* In Kelvin */ 141 { 142 ASSERT(lambda_min <= lambda_max && temperature > 0); 143 if(lambda_min == lambda_max) { 144 return htrdr_planck_monochromatic(lambda_min, temperature); 145 } else { 146 return htrdr_planck_interval(lambda_min, lambda_max, temperature); 147 } 148 } 149 150 BEGIN_DECLS 151 152 HTRDR_API res_T 153 htrdr_brightness_temperature 154 (struct htrdr* htrdr, 155 const double lambda_min, /* In meters */ 156 const double lambda_max, /* In meters */ 157 /* Averaged over [lambda_min, lambda_max]. In W/m²/sr/m */ 158 const double radiance, 159 double* temperature); 160 161 HTRDR_API double 162 htrdr_radiance_temperature 163 (struct htrdr* htrdr, 164 const double lambda_min, /* In meters */ 165 const double lambda_max, /* In meters */ 166 const double radiance); /* In W/m²/sr */ 167 168 END_DECLS 169 170 #endif /* HTRDR_SPECTRAL_H */