commit 7285d0c9f12e3be839e372cb94f01f4e2fc9bbad
parent a30b085f8224d7e411070cec55323a8baf12f392
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 29 Apr 2022 11:59:05 +0200
Provide Voigt profile function
Diffstat:
| M | cmake/CMakeLists.txt | | | 5 | +++-- |
| M | src/sln.h | | | 79 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- |
| A | src/sln_voigt.c | | | 199 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
3 files changed, 279 insertions(+), 4 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -48,7 +48,8 @@ set(SLN_FILES_SRC
sln_device.c
sln_log.c
sln_tree.c
- sln_tree_build.c)
+ sln_tree_build.c
+ sln_voigt.c)
set(SLN_FILES_INC
sln_device_c.h
sln_log.h)
@@ -65,7 +66,7 @@ rcmake_prepend_path(SLN_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
add_library(sln SHARED
${SLN_FILES_SRC} ${SLN_FILES_INC} ${SLN_FILES_INC_API})
-target_link_libraries(sln RSys StarHITRAN)
+target_link_libraries(sln RSys StarHITRAN m)
set_target_properties(sln PROPERTIES
DEFINE_SYMBOL SLN_SHARED_BUILD
diff --git a/src/sln.h b/src/sln.h
@@ -20,7 +20,9 @@
#define SLN_H
#include <rsys/rsys.h>
+
#include <float.h>
+#include <math.h>
/* Library symbol management */
#if defined(SLN_SHARED_BUILD) /* Build shared library */
@@ -77,7 +79,7 @@ struct sln_molecule {
size_t nisotopes; /* 0 <=> select all the isotopes of the molecules */
double concentration;
- double cutoff; /* in cm^1 */
+ double cutoff; /* in cm^-1 */
int32_t id; /* Identifier of the molecule into the isotope metadata */
};
@@ -112,10 +114,12 @@ static const struct sln_tree_create_args SLN_TREE_CREATE_ARGS_DEFAULT =
struct sln_tree_desc {
double wavenumber_range[2]; /* Spectral range */
+ double pressure; /* In atm */
+ double temperature; /* In K */
size_t nlines; /* Number of partitionned lines */
size_t max_nlines_per_leaf; /* max #lines per leaf */
};
-#define SLN_TREE_DESC_NULL__ {{0,0}, 0, 0}
+#define SLN_TREE_DESC_NULL__ {{0,0},0,0,0,0}
static const struct sln_tree_desc SLN_TREE_DESC_NULL = SLN_TREE_DESC_NULL__;
/* Forward declarations of opaque data structures */
@@ -186,4 +190,75 @@ sln_leaf_get_line
const size_t iline,
const struct shtr_line** line);
+/*******************************************************************************
+ * Helper functions
+ ******************************************************************************/
+/* Purpose: to calculate the Faddeeva function with relative error less than
+ * 10^(-4).
+ *
+ * Inputs: x and y, parameters for the Voigt function :
+ * - x is defined as x=(nu-nu_c)/gamma_D*sqrt(ln(2)) with nu the current
+ * wavenumber, nu_c the wavenumber at line center, gamma_D the Doppler
+ * linewidth.
+ * - y is defined as y=gamma_L/gamma_D*sqrt(ln(2)) with gamma_L the Lorentz
+ * linewith and gamma_D the Doppler linewidth
+ *
+ * Output: k, the Voigt function; it has to be multiplied by
+ * sqrt(ln(2)/pi)*1/gamma_D so that the result may be interpretable in terms of
+ * line profile.
+ *
+ * TODO check the copyright */
+SLN_API double
+sln_faddeeva
+ (const double x,
+ const double y);
+
+static INLINE double
+sln_compute_line_half_width_doppler
+ (const double nu, /* Line center wrt pressure in cm^-1 */ /* TODO check this */
+ const double molar_mass, /* In kg.mol^-1 */
+ const double temperature) /* In K */
+{
+ /* kb = 1.3806e-23
+ * Na = 6.02214179e23 TODO wrong avogadro constant
+ * c = 299792458
+ * sqrt(2*log(2)*kb*Na)/c */
+ const double sqrt_two_ln2_kb_Na_over_c = 1.1324432520993544423e-08;
+ const double gamma_d = nu * sqrt_two_ln2_kb_Na_over_c * sqrt(temperature/molar_mass);
+ ASSERT(nu > 0 && temperature >= 0 && molar_mass > 0);
+ return gamma_d;
+}
+
+static INLINE double
+sln_compute_line_half_width_lorentz
+ (const double gamma_air, /* Air broadening half width [cm^-1.atm^-1] */
+ const double gamma_self, /* Air broadening half width [cm^-1.atm^-1] */
+ const double pressure, /* [atm^-1] */
+ const double concentration)
+{
+ const double Ps = pressure * concentration;
+ const double gamma_l = (pressure - Ps) * gamma_air + Ps * gamma_self;
+ ASSERT(gamma_air > 0 && gamma_self > 0);
+ ASSERT(pressure > 0 && concentration >= 0 && concentration <= 1);
+ return gamma_l;
+}
+
+static INLINE double
+sln_compute_voigt_profile
+ (const double wavenumber, /* In cm^-1 */
+ const double nu, /* Line center in cm^-1 */
+ const double gamma_d, /* Doppler line half width in cm^-1 */
+ const double gamma_l) /* Lorentz line half width in cm^-1 */
+{
+ /* Constants */
+ const double sqrt_ln2 = 0.83255461115769768821; /* sqrt(log(2)) */
+ const double sqrt_ln2_over_pi = 0.46971863934982566180; /* sqrt(log(2)/M_PI) */
+ const double sqrt_ln2_over_gamma_d = sqrt_ln2 / gamma_d;
+
+ const double x = (wavenumber - nu) * sqrt_ln2_over_gamma_d;
+ const double y = gamma_l * sqrt_ln2_over_gamma_d;
+ const double k = sln_faddeeva(x, y);
+ return k*sqrt_ln2_over_pi/gamma_d;
+}
+
#endif /* SLN_H */
diff --git a/src/sln_voigt.c b/src/sln_voigt.c
@@ -0,0 +1,199 @@
+/* Copyright (C) 2022 CNRS - LMD
+ * Copyright (C) 2022 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2022 Université Paul Sabatier - IRIT
+ * Copyright (C) 2022 Université Paul Sabatier - Laplace
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sln.h"
+
+/*******************************************************************************
+ * Exported function
+ ******************************************************************************/
+double
+sln_faddeeva(const double x, const double Y)
+{
+ /* Constants */
+ const double RRTPI = 0.56418958;
+
+ /* For CPF12 algorithm */
+ const double Y0 = 1.5;
+ const double Y0PY0 = Y0+Y0;
+ const double Y0Q = Y0*Y0;
+
+ const double C[6] = {
+ 1.0117281,
+ -0.75197147,
+ 0.012557727,
+ 0.010022008,
+ -0.00024206814,
+ 0.00000050084806
+ };
+ const double S[6] = {
+ 1.393237,
+ 0.23115241,
+ -0.15535147,
+ 0.0062183662,
+ 0.000091908299,
+ -0.00000062752596
+ };
+ const double T[6] = {
+ 0.31424038,
+ 0.94778839,
+ 1.5976826,
+ 2.2795071,
+ 3.0206370,
+ 3.8897249
+ };
+
+ double ABX, XQ, YQ, YRRTPI;
+ double XLIM0, XLIM1, XLIM2, XLIM3, XLIM4;
+ double A0, D0, D2, E0, E2, E4, H0, H2, H4, H6;
+ double P0, P2, P4, P6, P8, Z0, Z2, Z4, Z6, Z8;
+ double XP[6], XM[6], YP[6], YM[6];
+ double MQ[6], PQ[6], MF[6], PF[6];
+ double D, YF, YPY0, YPY0Q;
+
+ /* Output */
+ double k = 0;
+
+ int i;
+ int RG1, RG2, RG3;
+
+ YQ = Y*Y;
+ YRRTPI = Y*RRTPI;
+
+ if(Y >= 70.55) {
+ XQ = x*x;
+ k=YRRTPI/(XQ+YQ);
+ return k;
+ }
+
+ RG1 = RG2 = RG3 = 1;
+
+ XLIM0 = sqrt(15100.0 + Y*(40.0 - Y*3.6));
+ if(Y >= 8.425){
+ XLIM1 = 0.0;
+ } else {
+ XLIM1 = sqrt(164.0 - Y*(4.3 + Y*1.8));
+ }
+
+ XLIM2 = 6.8-Y;
+ XLIM3 = 2.4*Y;
+ XLIM4 = 18.1*Y+1.65;
+
+ if(Y<=1.0e-6) {
+ XLIM1=XLIM0;
+ XLIM2=XLIM0;
+ }
+
+ ABX = sqrt(x*x);
+ XQ = ABX*ABX;
+ if(ABX >= XLIM0) {
+ k = YRRTPI/(XQ + YQ);
+ } else if (ABX >= XLIM1) {
+ if(RG1 != 0) {
+ RG1 = 0;
+ A0 = YQ + 0.5;
+ D0 = A0*A0;
+ D2 = YQ + YQ - 1.0;
+ }
+ D = RRTPI/(D0 + XQ*(D2 + XQ));
+ k = D*Y*(A0 + XQ);
+
+ } else if (ABX > XLIM2) {
+ if(RG2 != 0) {
+ RG2 = 0;
+ H0 = 0.5625 + YQ*(4.5 + YQ*(10.5 + YQ*(6.0 + YQ)));
+ H2 = -4.5 + YQ*(9.0 + YQ*(6.0 + YQ*4.0));
+ H4 = 10.5 - YQ*(6.0 - YQ*6.0);
+ H6 = -6.0 + YQ*4.0;
+ E0 = 1.875 + YQ*(8.25 + YQ*(5.5+YQ));
+ E2 = 5.25 + YQ*(1.0 + YQ*3.0);
+ E4 = 0.75*H6;
+ }
+ D = RRTPI/(H0 + XQ*(H2 + XQ*(H4 + XQ*(H6 + XQ))));
+ k = D*Y*(E0 + XQ*(E2 + XQ*(E4 + XQ)));
+
+ } else if (ABX<XLIM3) {
+ if(RG3 != 0) {
+ RG3 = 0;
+ Z0 = 272.1014
+ + Y*(1280.829 + Y*(2802.870 + Y*(3764.966
+ + Y*(3447.629 + Y*(2256.981 + Y*(1074.409
+ + Y*(369.1989 + Y*(88.26741 + Y*(13.39880 + Y)))))))));
+ Z2 = 211.678
+ + Y*(902.3066 + Y*(1758.336 + Y*(2037.310
+ + Y*(1549.675 + Y*(793.4273 + Y*(266.2987
+ + Y*(53.59518 + Y*5.0)))))));
+ Z4 = 78.86585
+ + Y*(308.1852 + Y*(497.3014 + Y*(479.2576
+ + Y*(269.2916 + Y*(80.39278 + Y*10.0)))));
+ Z6 = 22.03523
+ + Y*(55.02933 + Y*(92.75679 + Y*(53.59518 + Y*10.0)));
+ Z8 = 1.496460
+ + Y*(13.39880 + Y*5.0);
+ P0 = 153.5168
+ + Y*(549.3954 + Y*(919.4955 + Y*(946.8970
+ + Y*(662.8097 + Y*(328.2151 + Y*(115.3772
+ + Y*(27.93941 + Y*(4.264678 + Y*0.3183291))))))));
+ P2 = -34.16955
+ + Y*(-1.322256+ Y*(124.5975 + Y*(189.7730
+ + Y*(139.4665 + Y*(56.81652 + Y*(12.79458 + Y*1.2733163))))));
+ P4 = 2.584042
+ + Y*(10.46332 + Y*(24.01655 + Y*(29.81482
+ + Y*(12.79568 + Y*1.9099744))));
+ P6 = -0.07272979
+ + Y*(0.9377051 + Y*(4.266322 + Y*1.273316));
+ P8 = 0.0005480304
+ + Y*0.3183291;
+ }
+ D = 1.7724538/(Z0 + XQ*(Z2 + XQ*(Z4 + XQ*(Z6 + XQ*(Z8 + XQ)))));
+ k = D*(P0 + XQ*(P2 + XQ*(P4 + XQ*(P6 + XQ*P8))));
+
+ } else {
+ YPY0 = Y+Y0;
+ YPY0Q = YPY0*YPY0;
+ k=0.0;
+
+ FOR_EACH(i, 0, 6) {
+ D = x - T[i];
+ MQ[i] = D*D;
+ MF[i] = 1.0/(MQ[i] + YPY0Q);
+ XM[i] = MF[i]*D;
+ YM[i] = MF[i]*YPY0;
+
+ D = x + T[i];
+ PQ[i] = D*D;
+ PF[i] = 1.0/(PQ[i] + YPY0Q);
+ XP[i] = PF[i]*D;
+ YP[i] = PF[i]*YPY0;
+ }
+
+ if(ABX <= XLIM4) {
+ FOR_EACH(i, 0, 6) {
+ k = k + C[i]*(YM[i] + YP[i]) - S[i]*(XM[i] - XP[i]);
+ }
+ } else {
+ YF = Y+Y0PY0;
+ FOR_EACH(i, 0, 6) {
+ k = k
+ + (C[i]*(MQ[i]*MF[i] - Y0*YM[i]) + S[i]*YF*XM[i])/(MQ[i] + Y0Q)
+ + (C[i]*(PQ[i]*PF[i] - Y0*YP[i]) - S[i]*YF*XP[i])/(PQ[i] + Y0Q);
+ }
+ k = Y*k + exp(-XQ);
+ }
+ }
+ return k;
+}
