ARTS 2.5.4 (git: 31ce4f0e)
mc_interp.cc
Go to the documentation of this file.
1/* Copyright (C) 2005-2012 Cory Davis <cory@met.ed.ac.uk>
2
3 This program is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 USA. */
17
26/*===========================================================================
27 === External declarations
28 ===========================================================================*/
29#include "mc_interp.h"
30#include "arts_constants.h"
31#include "arts_conversions.h"
32#include "logic.h"
33#include "montecarlo.h"
34
35
38inline constexpr Numeric PI=Constant::pi;
39
41 GridPos gp1, gpl, gpr;
42 Vector itw1(2), itwl(2), itwr(2);
43 Numeric yl, yr;
44 //Get interpolation weights for x1 interpolation
45 gridpos(gp1, this->x1a, x1);
46 interpweights(itw1, gp1);
47 //Get y values on bounding x1 grid points for desired x2
48 gridpos(gpl, this->x2a[gp1.idx], x2);
49 interpweights(itwl, gpl);
50 gridpos(gpr, this->x2a[gp1.idx + 1], x2);
51 interpweights(itwr, gpr);
52 yl = interp(itwl, this->ya[gp1.idx], gpl);
53 yr = interp(itwr, this->ya[gp1.idx + 1], gpr);
54 //interpolate these two y values useing x1 interpolation weights
55 return itw1[0] * yl + itw1[1] * yr;
56}
57
58//void SLIData2::check() const
59//{
60// Index nx1=this->x1a.nelem();
61// ARTS_ASSERT(nx1>0);
62//}
63
64ostream& operator<<(ostream& os, const SLIData2& /* sli */) {
65 os << "SLIData2 : Output operator not implemented";
66 return os;
67}
68
71 const ArrayOfMatrix& a,
72 const GridPos& tc) {
73 DEBUG_ONLY(const Numeric sum_check_epsilon = 1e-6);
74
75 ARTS_ASSERT(is_size(itw, 2)); // We need 2 interpolation
76 // weights.
77
78 // Check that interpolation weights are valid. The sum of all
79 // weights (last dimension) must always be approximately one.
80 ARTS_ASSERT(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
81
82 Index anr = a[0].nrows();
83 Index anc = a[0].ncols();
84
85 ARTS_ASSERT(tia.nrows() == anr);
86 ARTS_ASSERT(tia.ncols() == anc);
87
88 for (Index inr = 0; inr < anr; inr++)
89 for (Index inc = 0; inc < anc; inc++) {
90 tia(inr, inc) =
91 a[tc.idx](inr, inc) * itw[0] + a[tc.idx + 1](inr, inc) * itw[1];
92 }
93}
94
97 const ArrayOfVector& a,
98 const GridPos& tc) {
99 DEBUG_ONLY(const Numeric sum_check_epsilon = 1e-6);
100 ARTS_ASSERT(is_size(itw, 2)); // We need 2 interpolation
101 // weights.
102
103 // Check that interpolation weights are valid. The sum of all
104 // weights (last dimension) must always be approximately one.
105 ARTS_ASSERT(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
106
107 Index an = a[0].nelem();
108
109 ARTS_ASSERT(tia.nelem() == an);
110
111 for (Index i = 0; i < an; ++i) {
112 tia[i] = a[tc.idx][i] * itw[0] + a[tc.idx + 1][i] * itw[1];
113 }
114}
115
117 VectorView pha_mat_int,
118 Numeric& theta_rad,
119 //Input:
120 const SingleScatteringData& scat_data_single,
121 const Numeric& za_sca,
122 const Numeric& aa_sca,
123 const Numeric& za_inc,
124 const Numeric& aa_inc,
125 const Numeric& rtp_temperature) {
126 Numeric ANG_TOL = 1e-7;
127
128 //Calculate scattering angle from incident and scattered directions.
129 //The two special cases are implemented here to avoid NaNs that can
130 //sometimes occur in in the acos... formula in forward and backscatter
131 //cases. CPD 5/10/03.
132
133 if (abs(aa_sca - aa_inc) < ANG_TOL) {
134 theta_rad = DEG2RAD * abs(za_sca - za_inc);
135 } else if (abs(abs(aa_sca - aa_inc) - 180) < ANG_TOL) {
136 theta_rad = DEG2RAD * (za_sca + za_inc);
137 if (theta_rad > PI) {
138 theta_rad = 2 * PI - theta_rad;
139 }
140 } else {
141 const Numeric za_sca_rad = za_sca * DEG2RAD;
142 const Numeric za_inc_rad = za_inc * DEG2RAD;
143 const Numeric aa_sca_rad = aa_sca * DEG2RAD;
144 const Numeric aa_inc_rad = aa_inc * DEG2RAD;
145
146 // cout << "Interpolation on scattering angle" << endl;
147 ARTS_ASSERT(scat_data_single.pha_mat_data.ncols() == 6);
148 // Calculation of the scattering angle:
149 theta_rad =
150 acos(cos(za_sca_rad) * cos(za_inc_rad) +
151 sin(za_sca_rad) * sin(za_inc_rad) * cos(aa_sca_rad - aa_inc_rad));
152 }
153 const Numeric theta = RAD2DEG * theta_rad;
154
155 // Interpolation of the data on the scattering angle:
156
157 GridPos thet_gp;
158 gridpos(thet_gp, scat_data_single.za_grid, theta);
159 GridPos t_gp;
160
161 if (scat_data_single.T_grid.nelem() == 1) {
162 Vector itw(2);
163 interpweights(itw, thet_gp);
164
165 for (Index i = 0; i < 6; i++) {
166 pha_mat_int[i] = interp(
167 itw, scat_data_single.pha_mat_data(0, 0, joker, 0, 0, 0, i), thet_gp);
168 }
169 } else {
170 gridpos(t_gp, scat_data_single.T_grid, rtp_temperature);
171
172 Vector itw(4);
173 interpweights(itw, t_gp, thet_gp);
174
175 for (Index i = 0; i < 6; i++) {
176 pha_mat_int[i] =
177 interp(itw,
178 scat_data_single.pha_mat_data(0, joker, joker, 0, 0, 0, i),
179 t_gp,
180 thet_gp);
181 }
182 }
183}
Constants of physical expressions as constexpr.
Common ARTS conversions.
This can be used to make arrays out of anything.
Definition: array.h:48
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
Index ncols() const noexcept
Definition: matpackVII.h:159
A constant view of a Vector.
Definition: matpackI.h:512
Numeric sum() const ARTS_NOEXCEPT
The sum of all elements of a Vector.
Definition: matpackI.cc:48
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The MatrixView class.
Definition: matpackI.h:1164
A 2D sequential linear interpolation (SLI) lookup table This class holds the gridded for 2D SLI as we...
Definition: mc_interp.h:52
ArrayOfVector x2a
Definition: mc_interp.h:57
Numeric interpolate(Numeric x1, Numeric x2) const
Definition: mc_interp.cc:40
ArrayOfVector ya
Definition: mc_interp.h:59
Vector x1a
Definition: mc_interp.h:55
The VectorView class.
Definition: matpackI.h:663
The Vector class.
Definition: matpackI.h:899
#define DEBUG_ONLY(...)
Definition: debug.h:52
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
#define abs(x)
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:82
bool is_same_within_epsilon(const Numeric &a, const Numeric &b, const Numeric &epsilon)
Check, if two numbers agree within a given epsilon.
Definition: logic.cc:354
Header file for logic.cc.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
constexpr Numeric DEG2RAD
Definition: mc_interp.cc:36
void interp(MatrixView tia, ConstVectorView itw, const ArrayOfMatrix &a, const GridPos &tc)
interp.
Definition: mc_interp.cc:69
ostream & operator<<(ostream &os, const SLIData2 &)
Definition: mc_interp.cc:64
constexpr Numeric RAD2DEG
Definition: mc_interp.cc:37
void interp_scat_angle_temperature(VectorView pha_mat_int, Numeric &theta_rad, const SingleScatteringData &scat_data_single, const Numeric &za_sca, const Numeric &aa_sca, const Numeric &za_inc, const Numeric &aa_inc, const Numeric &rtp_temperature)
interp_scat_angle_temperature.
Definition: mc_interp.cc:116
constexpr Numeric PI
Definition: mc_interp.cc:38
Interpolation classes and functions created for use within Monte Carlo scattering simulations.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric e
Elementary charge convenience name [C].
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
constexpr auto rad2deg(auto x) noexcept
Converts radians to degrees.
Structure to store a grid position.
Definition: interpolation.h:73
Index idx
Definition: interpolation.h:74
#define a