ARTS  2.4.0(git:4fb77825)
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 "logic.h"
31 #include "montecarlo.h"
32 
34  GridPos gp1, gpl, gpr;
35  Vector itw1(2), itwl(2), itwr(2);
36  Numeric yl, yr;
37  //Get interpolation weights for x1 interpolation
38  gridpos(gp1, this->x1a, x1);
39  interpweights(itw1, gp1);
40  //Get y values on bounding x1 grid points for desired x2
41  gridpos(gpl, this->x2a[gp1.idx], x2);
42  interpweights(itwl, gpl);
43  gridpos(gpr, this->x2a[gp1.idx + 1], x2);
44  interpweights(itwr, gpr);
45  yl = interp(itwl, this->ya[gp1.idx], gpl);
46  yr = interp(itwr, this->ya[gp1.idx + 1], gpr);
47  //interpolate these two y values useing x1 interpolation weights
48  return itw1[0] * yl + itw1[1] * yr;
49 }
50 
51 //void SLIData2::check() const
52 //{
53 // Index nx1=this->x1a.nelem();
54 // assert(nx1>0);
55 //}
56 
57 ostream& operator<<(ostream& os, const SLIData2& /* sli */) {
58  os << "SLIData2 : Output operator not implemented";
59  return os;
60 }
61 
62 void interp(MatrixView tia,
63  ConstVectorView itw,
64  const ArrayOfMatrix& a,
65  const GridPos& tc) {
66  DEBUG_ONLY(const Numeric sum_check_epsilon = 1e-6);
67 
68  assert(is_size(itw, 2)); // We need 2 interpolation
69  // weights.
70 
71  // Check that interpolation weights are valid. The sum of all
72  // weights (last dimension) must always be approximately one.
73  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
74 
75  Index anr = a[0].nrows();
76  Index anc = a[0].ncols();
77 
78  assert(tia.nrows() == anr);
79  assert(tia.ncols() == anc);
80 
81  for (Index inr = 0; inr < anr; inr++)
82  for (Index inc = 0; inc < anc; inc++) {
83  tia(inr, inc) =
84  a[tc.idx](inr, inc) * itw[0] + a[tc.idx + 1](inr, inc) * itw[1];
85  }
86 }
87 
88 void interp(VectorView tia,
89  ConstVectorView itw,
90  const ArrayOfVector& a,
91  const GridPos& tc) {
92  DEBUG_ONLY(const Numeric sum_check_epsilon = 1e-6);
93  assert(is_size(itw, 2)); // We need 2 interpolation
94  // weights.
95 
96  // Check that interpolation weights are valid. The sum of all
97  // weights (last dimension) must always be approximately one.
98  assert(is_same_within_epsilon(itw.sum(), 1, sum_check_epsilon));
99 
100  Index an = a[0].nelem();
101 
102  assert(tia.nelem() == an);
103 
104  for (Index i = 0; i < an; ++i) {
105  tia[i] = a[tc.idx][i] * itw[0] + a[tc.idx + 1][i] * itw[1];
106  }
107 }
108 
110  VectorView pha_mat_int,
111  Numeric& theta_rad,
112  //Input:
114  const Numeric& za_sca,
115  const Numeric& aa_sca,
116  const Numeric& za_inc,
117  const Numeric& aa_inc,
118  const Numeric& rtp_temperature) {
119  Numeric ANG_TOL = 1e-7;
120 
121  //Calculate scattering angle from incident and scattered directions.
122  //The two special cases are implemented here to avoid NaNs that can
123  //sometimes occur in in the acos... formula in forward and backscatter
124  //cases. CPD 5/10/03.
125 
126  if (abs(aa_sca - aa_inc) < ANG_TOL) {
127  theta_rad = DEG2RAD * abs(za_sca - za_inc);
128  } else if (abs(abs(aa_sca - aa_inc) - 180) < ANG_TOL) {
129  theta_rad = DEG2RAD * (za_sca + za_inc);
130  if (theta_rad > PI) {
131  theta_rad = 2 * PI - theta_rad;
132  }
133  } else {
134  const Numeric za_sca_rad = za_sca * DEG2RAD;
135  const Numeric za_inc_rad = za_inc * DEG2RAD;
136  const Numeric aa_sca_rad = aa_sca * DEG2RAD;
137  const Numeric aa_inc_rad = aa_inc * DEG2RAD;
138 
139  // cout << "Interpolation on scattering angle" << endl;
140  assert(scat_data_single.pha_mat_data.ncols() == 6);
141  // Calculation of the scattering angle:
142  theta_rad =
143  acos(cos(za_sca_rad) * cos(za_inc_rad) +
144  sin(za_sca_rad) * sin(za_inc_rad) * cos(aa_sca_rad - aa_inc_rad));
145  }
146  const Numeric theta = RAD2DEG * theta_rad;
147 
148  // Interpolation of the data on the scattering angle:
149 
150  GridPos thet_gp;
151  gridpos(thet_gp, scat_data_single.za_grid, theta);
152  GridPos t_gp;
153 
154  if (scat_data_single.T_grid.nelem() == 1) {
155  Vector itw(2);
156  interpweights(itw, thet_gp);
157 
158  for (Index i = 0; i < 6; i++) {
159  pha_mat_int[i] = interp(
160  itw, scat_data_single.pha_mat_data(0, 0, joker, 0, 0, 0, i), thet_gp);
161  }
162  } else {
164 
165  Vector itw(4);
166  interpweights(itw, t_gp, thet_gp);
167 
168  for (Index i = 0; i < 6; i++) {
169  pha_mat_int[i] =
170  interp(itw,
171  scat_data_single.pha_mat_data(0, joker, joker, 0, 0, 0, i),
172  t_gp,
173  thet_gp);
174  }
175  }
176 }
ARTS::Var::scat_data_single
SingleScatteringData scat_data_single(Workspace &ws) noexcept
Definition: autoarts.h:5926
SLIData2::x1a
Vector x1a
Definition: mc_interp.h:59
is_same_within_epsilon
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:351
interp
void interp(MatrixView tia, ConstVectorView itw, const ArrayOfMatrix &a, const GridPos &tc)
interp.
Definition: mc_interp.cc:62
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
operator<<
ostream & operator<<(ostream &os, const SLIData2 &)
Definition: mc_interp.cc:57
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:156
mc_interp.h
Interpolation classes and functions created for use within Monte Carlo scattering simulations.
RAD2DEG
const Numeric RAD2DEG
Global constant, conversion from radians to degrees.
joker
const Joker joker
abs
#define abs(x)
Definition: legacy_continua.cc:20626
DEBUG_ONLY
#define DEBUG_ONLY(...)
Definition: debug.h:36
SLIData2::x2a
ArrayOfVector x2a
Definition: mc_interp.h:61
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
interp_scat_angle_temperature
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:109
DEG2RAD
const Numeric DEG2RAD
Global constant, conversion from degrees to radians.
montecarlo.h
is_size
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:79
SingleScatteringData
Definition: optproperties.h:80
Array
This can be used to make arrays out of anything.
Definition: array.h:108
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
VectorView
The VectorView class.
Definition: matpackI.h:610
ConstVectorView::sum
Numeric sum() const
The sum of all elements of a Vector.
Definition: matpackI.cc:53
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
SLIData2::ya
ArrayOfVector ya
Definition: mc_interp.h:63
GridPos
Structure to store a grid position.
Definition: interpolation.h:73
ARTS::Var::rtp_temperature
Numeric rtp_temperature(Workspace &ws) noexcept
Definition: autoarts.h:5807
logic.h
Header file for logic.cc.
x1
#define x1
Definition: lineshapemodel.h:414
SLIData2
A 2D sequential linear interpolation (SLI) lookup table This class holds the gridded for 2D SLI as we...
Definition: mc_interp.h:56
SLIData2::interpolate
Numeric interpolate(Numeric x1, Numeric x2) const
Definition: mc_interp.cc:33
GridPos::idx
Index idx
Definition: interpolation.h:74
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Vector
The Vector class.
Definition: matpackI.h:860
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
x2
#define x2
Definition: lineshapemodel.h:415
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
PI
const Numeric PI
Global constant, pi.