38#include <Faddeeva/Faddeeva.hh>
66 const Numeric Const2 =
b2 + c2 + d2 - u2 - v2 - w2;
69 Const1 =
b2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
70 Const1 += c2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
71 Const1 += d2 * (d2 * 0.5 + u2 - v2 - w2);
72 Const1 += u2 * (u2 * 0.5 + v2 + w2);
73 Const1 += v2 * (v2 * 0.5 + w2);
75 Const1 += 8 * (
b *
d *
u *
w -
b *
c *
v *
w -
c *
d *
u *
v);
79 Const1 =
sqrt(Const1);
94 const Numeric inv_x2y2 = 1.0 / x2y2;
96 std::cout << x <<
" " << y <<
" " << Const1 <<
" " << Const2 <<
"\n";
99 Numeric inv_y = 0.0, inv_x = 0.0;
106 C2 = (1.0 - cos_y) * inv_x2y2;
107 C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
108 }
else if (y == 0.0) {
112 C2 = (cosh_x - 1.0) * inv_x2y2;
113 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
118 C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
119 C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
120 C2 = (cosh_x - cos_y) * inv_x2y2;
121 C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
124 std::cout << C0 <<
" " << C1 <<
" " << C2 <<
" " << C3 <<
"\n";
126 Matrix F(4, 4, 0), A(4, 4, 0);
129 Eigen::Matrix4d eigA;
130 eigA << 0,
b,
c,
d,
b, 0,
u,
v,
c, -
u, 0,
w,
d, -
v, -
w, 0;
132 eigF = C1 * eigA + C2 * eigA * eigA + C3 * eigA * eigA * eigA;
139 std::cout << F <<
"\n";
145 std::cout <<
"Initialized TransmissionMatrix(2, 4):\n" <<
a <<
"\n";
149 A << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16;
150 std::cout <<
"New Matrix:\n" << A <<
"\n\n";
152 std::cout <<
"Updated TransmissionMatrix Position 1 wit New Matrix:\n"
157 "1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 125 26 27 28 29 30 31 32";
158 std::cout <<
"Stream:\n" << S <<
"\n\n";
159 std::istringstream astream(S);
161 std::cout <<
"Streamed into TransmissionMatrix:\n" <<
a <<
"\n";
165 std::cout <<
"Initialized RadiationVector(3, 3)\n" <<
b <<
"\n";
170 std::cout <<
"New Vector:\n" << B <<
"\n\n";
171 b.Vec3(1).noalias() += B;
172 std::cout <<
"Updated RadiationVector Position 1 with New Vector:\n"
176 String T =
"1 2 3 4 5 6 7 8 90";
177 std::cout <<
"Stream:\n" << T <<
"\n\n";
178 std::istringstream bstream(T);
180 std::cout <<
"Streamed into RadiationVector:\n" <<
b <<
"\n";
184 auto f = [](
Numeric x) {
return 0.1 * x; };
185 auto df = [](
Numeric x) {
return 0.1 + 0 * x; };
191 const Numeric r_extra1 = r_normal + f(x1);
192 const Numeric r_extra2 = r_normal + f(x2);
195 a.Kjj() = 10 +
Numeric(rand() % 1000) / 100;
197 a.K12() = 2 +
Numeric(rand() % 1000) / 100;
199 a.K13() = 3 +
Numeric(rand() % 1000) / 100;
200 a.K23() = -1 +
Numeric(rand() % 1000) / 100;
202 a.K14() = 5 +
Numeric(rand() % 1000) / 100;
203 a.K24() = -3 +
Numeric(rand() % 1000) / 100;
204 a.K34() = -2 +
Numeric(rand() % 1000) / 100;
211 b.Kjj() = 5 +
Numeric(rand() % 1000) / 100;
213 b.K12() = -1 +
Numeric(rand() % 1000) / 100;
215 b.K13() = -3 +
Numeric(rand() % 1000) / 100;
216 b.K23() = 4 +
Numeric(rand() % 1000) / 100;
218 b.K14() = 2 +
Numeric(rand() % 1000) / 100;
219 b.K24() = -1 +
Numeric(rand() % 1000) / 100;
220 b.K34() = 3 +
Numeric(rand() % 1000) / 100;
228 Tensor3 T_normal(1, nstokes, nstokes), T_extra(1, nstokes, nstokes);
229 Tensor4 dT1(1, 1, nstokes, nstokes), dT2(1, 1, nstokes, nstokes);
232 T_normal, dT1, dT2, r_normal,
a,
b, da, da, df(x1), df(x2), 0);
234 std::cout <<
"Transmission at r=" << r_normal <<
":\n"
237 std::cout <<
"First derivative:\n"
240 std::cout <<
"Second derivative:\n"
246 std::cout <<
"Transmission at perturbed r1=" << r_extra1 <<
":\n"
251 std::cout <<
"First derivative perturbed:\n"
255 std::cout <<
"First derivative perturbed relative:\n"
261 std::cout <<
"Transmission at perturbed r2=" << r_extra2 <<
":\n"
266 std::cout <<
"Second derivative perturbed:\n"
270 std::cout <<
"Second derivative perturbed relative:\n"
303 std::cout <<
test1 <<
"\n\n"
363 std::cout << ans1 <<
"\n\n"
372 for (
auto& pm : propmats) {
390 std::cout <<
"Propmats:\n" << propmats <<
"\n\n";
394 for (i = 0; i < 4; i++) {
408 std::cout <<
"Layers:\n" << layers <<
"\n\n";
415 std::cout <<
"Forward accumulation:\n" << cumulative_forward <<
"\n\n";
416 std::cout <<
"Reflect accumulation:\n" << cumulative_reflect <<
"\n\n";
424 Vector x(n), sx(n), shx(n), cx(n), chx(n);
427 for (
int i = 0; i < n; i++) {
428 sx[i] = std::sin(x[i]);
429 cx[i] = std::cos(x[i]);
430 shx[i] = std::sinh(x[i]);
431 chx[i] = std::cosh(x[i]);
435 << std::scientific << std::setprecision(15)
436 <<
"x\tabs(sx/x-1)\tabs(shx/x-1)\tabs((chx-cx)/2x^2-1/2)\tabs((shx/x-sx/x)/2x^2-1/6)\n";
438 for (
int i = 0; i < n; i++)
439 std::cout << x[i] <<
'\t' <<
std::abs(sx[i] / x[i] - 1.0) <<
'\t'
440 <<
std::abs(shx[i] / x[i] - 1.0) <<
'\t'
441 <<
std::abs((chx[i] - cx[i]) / (2 * x[i] * x[i]) - 0.5) <<
'\t'
442 <<
std::abs((shx[i] / x[i] - sx[i] / x[i]) / (2 * x[i] * x[i]) -
451 constexpr Index nf = 501;
452 constexpr Numeric fstart = 25e9;
453 constexpr Numeric fend = 165e9;
476 constexpr auto df = 1000;
477 constexpr auto dt = 0.1;
479 Matrix pxsec_dt(nf, 1, 0);
480 Matrix pxsec_df(nf, 1, 0);
483 PWR93O2AbsModel(pxsec, 0, 1, 1, 1,
"user",
"PWR98", f, {p}, {t}, {0.5}, {1},
Verbosity());
484 PWR93O2AbsModel(pxsec_dt, 0, 1, 1, 1,
"user",
"PWR98", f, {p}, {t+dt}, {0.5}, {1},
Verbosity());
485 PWR93O2AbsModel(pxsec_df, 0, 1, 1, 1,
"user",
"PWR98", f_pert, {p}, {t}, {0.5}, {1},
Verbosity());
487 std::cout<<
"xr = np.array([";
488 for (
Index i=0; i<nf; i++)
489 std::cout<<
'['<<xsec(i, 0)<<
','<<
' '<<pxsec(i, 0) <<
']' <<
',' <<
' ';
490 std::cout <<
']' <<
')' <<
'\n';
492 std::cout<<
"dxr_dt = np.array([";
493 for (
Index i=0; i<nf; i++)
494 std::cout<<
'['<<dxsec[0](i, 0)<<
','<<
' '<<(pxsec_dt(i, 0)-pxsec(i, 0))/dt <<
']' <<
',' <<
' ';
495 std::cout <<
']' <<
')' <<
'\n';
497 std::cout<<
"dxr_df = np.array([";
498 for (
Index i=0; i<nf; i++)
499 std::cout<<
'['<<dxsec[1](i, 0)<<
','<<
' '<<(pxsec_df(i, 0)-pxsec(i, 0))/df <<
']' <<
',' <<
' ';
500 std::cout <<
']' <<
')' <<
'\n';
512 const Numeric stotmax = 0.1e-21;
514 const Index nsig =
Index(((sigmax - sigmin) / dsig) + 0.5) + 1;
518 for (
Index isig=0; isig<nsig; isig++) {
520 invcm_grid[isig] = sigc;
526 std::array<std::pair<lm_hitran_2017::ModeOfLineMixing, lm_hitran_2017::calctype>, 6>{
539 for (
Index i=0;i<n; i++) {
543 Vector vmrs = {1-xco2/100-xh2o/100, xh2o/100, xco2/100};
551 for (
Index i=0; i<nsig; i++) {
552 for (
Index j=0; j<n; j++) {
553 std::cout<<absorption[j][i]<<
' ';
568 if (n == 2 and
String(argc[1]) ==
"new") {
569 std::cout<<
"new test\n";
573 std::cout<<
"old test\n";
Declarations required for the calculation of absorption coefficients.
The global header file for ARTS.
This can be used to make arrays out of anything.
void PWR93O2AbsModel(MatrixView pxsec, const Numeric CCin, const Numeric CLin, const Numeric CWin, const Numeric COin, const String &model, const String &version, ConstVectorView f_grid, ConstVectorView abs_p, ConstVectorView abs_t, ConstVectorView vmrh2o, ConstVectorView vmr, const Verbosity &verbosity)
Oxygen complex at 60 GHz plus mm O2 lines plus O2 continuum.
MatrixViewMap MapToEigen(Tensor3View &A, const Index &i)
Namespace and functions to deal with HITRAN linemixing.
Constains various line scaling functions.
Contains the line shape namespace.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
void nlogspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlogspace
Eigen::Map< MatrixType, 0, StrideType > MatrixViewMap
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
std::complex< Numeric > Complex
my_basic_string< char > String
The String type for ARTS.
void compute(PropagationMatrix &propmat_clearsky, ArrayOfPropagationMatrix &dpropmat_clearsky_dx, const SpeciesIsotopeRecord &model, const Vector &f_grid, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const VMRS &vmr, const ArrayOfRetrievalQuantity &jacobian_quantities)
Compute the predefined model.
constexpr auto atm2pa(T x) noexcept -> decltype(x *101 '325.0)
Conversion from Atm to Pa.
constexpr auto kaycm_per_cmsquared2hz_per_msquared(T x) noexcept -> decltype(x *kaycm2freq(1e-4))
Conversion from cm-1 per molecule per cm^2 to Hz per molecule per m^2.
constexpr auto kaycm2freq(T x) noexcept -> decltype(x *(100 *c))
Conversion from Kayser wavenumber to Hz.
SpeciesIsotopologueRatios isotopologue_ratios()
Vector vmrs(const ConstVectorView &atmospheric_vmrs, const ArrayOfArrayOfSpeciesTag &atmospheric_species, const ArrayOfSpecies &lineshape_species) ARTS_NOEXCEPT
Returns a VMR vector for this model's main calculations.
constexpr std::array Isotopologues
A list of all ARTS isotopologues, note how the species enum class input HAS to be sorted.
constexpr Index find_species_index(const Species spec, const std::string_view isot) noexcept
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
constexpr Rational end(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the largest M for a polarization type of this transition.
void read(HitranRelaxationMatrixData &hitran, ArrayOfAbsorptionLines &bands, const SpeciesIsotopologueRatios &isotopologue_ratio, const String &basedir, const Numeric linemixinglimit, const Numeric fmin, const Numeric fmax, const Numeric stot, const ModeOfLineMixing mode)
Read from HITRAN online line mixing file.
Vector compute(const Numeric p, const Numeric t, const Numeric xco2, const Numeric xh2o, const ConstVectorView &invcm_grid, const Numeric stotmax, const calctype type)
void compute_transmission_matrix_and_derivative(Tensor3View T, Tensor4View dT_dx_upper_level, Tensor4View dT_dx_lower_level, const Numeric &r, const PropagationMatrix &upper_level, const PropagationMatrix &lower_level, const ArrayOfPropagationMatrix &dupper_level_dx, const ArrayOfPropagationMatrix &dlower_level_dx, const Numeric &dr_dTu, const Numeric &dr_dTl, const Index it, const Index iz, const Index ia)
void compute_transmission_matrix(Tensor3View T, const Numeric &r, const PropagationMatrix &upper_level, const PropagationMatrix &lower_level, const Index iz, const Index ia)
Compute the matrix exponent as the transmission matrix of this propagation matrix.
Array< PropagationMatrix > ArrayOfPropagationMatrix
Numeric sqrt(const Rational r)
Square root.
Contains known required VMR values.
Holds all information required for individual partial derivatives.
Radiation Vector for Stokes dimension 1-4.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
void test_transmat_to_cumulativetransmat()
void test_r_deriv_propagationmatrix()
void test_sinc_likes_0limit()
int main(int n, char **argc)
void test_hitran2017(bool newtest=true)
void test_transmissionmatrix()
void test_matrix_buildup()
void test_transmat_from_propmat()
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
void stepwise_transmission(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dtemp1, const Numeric &dr_dtemp2, const Index temp_deriv_pos)
Set the stepwise transmission matrix.
Stuff related to the transmission matrix.
Index make_wigner_ready(int largest, int fastest, int size)
Ready Wigner.
Wigner symbol interactions.
Headers and class definition of Zeeman modeling.