htrdr_spectral.c (3666B)
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 #include "core/htrdr.h" 25 #include "core/htrdr_log.h" 26 #include "core/htrdr_spectral.h" 27 28 /******************************************************************************* 29 * Exported symbols 30 ******************************************************************************/ 31 res_T 32 htrdr_brightness_temperature 33 (struct htrdr* htrdr, 34 const double lambda_min, 35 const double lambda_max, 36 const double radiance, /* In W/m2/sr/m */ 37 double* temperature) 38 { 39 const size_t MAX_ITER = 100; 40 const double epsilon_T = 1e-4; /* In K */ 41 const double epsilon_B = radiance * 1e-8; 42 double T, T0, T1, T2; 43 double B, B0; 44 size_t i; 45 res_T res = RES_OK; 46 ASSERT(temperature && lambda_min <= lambda_max); 47 48 /* Search for a brightness temperature whose radiance is greater than or 49 * equal to the estimated radiance */ 50 T2 = 200; 51 FOR_EACH(i, 0, MAX_ITER) { 52 const double B2 = htrdr_planck(lambda_min, lambda_max, T2); 53 if(B2 >= radiance) break; 54 T2 *= 2; 55 } 56 if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; } 57 58 B0 = T0 = T1 = 0; 59 FOR_EACH(i, 0, MAX_ITER) { 60 T = (T1+T2)*0.5; 61 B = htrdr_planck(lambda_min, lambda_max, T); 62 63 if(B < radiance) { 64 T1 = T; 65 } else { 66 T2 = T; 67 } 68 69 if(fabs(T-T0) < epsilon_T || fabs(B-B0) < epsilon_B) 70 break; 71 72 T0 = T; 73 B0 = B; 74 } 75 if(i >= MAX_ITER) { res = RES_BAD_OP; goto error; } 76 77 *temperature = T; 78 79 exit: 80 return res; 81 error: 82 htrdr_log_err(htrdr, 83 "Could not compute the brightness temperature for the estimated radiance %g " 84 "averaged over [%g, %g] nanometers.\n", 85 radiance, 86 lambda_min*1e9, 87 lambda_max*1e9); 88 goto exit; 89 } 90 91 double 92 htrdr_radiance_temperature 93 (struct htrdr* htrdr, 94 const double lambda_min, /* In meters */ 95 const double lambda_max, /* In meters */ 96 const double radiance) /* In W/m^2/sr */ 97 { 98 double temperature = 0; 99 double radiance_avg = radiance; 100 res_T res = RES_OK; 101 ASSERT(htrdr && radiance >= 0); 102 103 /* From integrated radiance to average radiance in W/m^2/sr/m */ 104 if(lambda_min != lambda_max) { /* !monochromatic */ 105 radiance_avg /= (lambda_max - lambda_min); 106 } 107 108 res = htrdr_brightness_temperature 109 (htrdr, 110 lambda_min, 111 lambda_max, 112 radiance_avg, 113 &temperature); 114 if(res != RES_OK) { 115 htrdr_log_warn(htrdr, 116 "Could not compute the brightness temperature for the radiance %g.\n", 117 radiance_avg); 118 temperature = 0; 119 } 120 return temperature; 121 } 122