ARTS  2.4.0(git:4fb77825)
surface.cc
Go to the documentation of this file.
1 /* Copyright (C) 2012
2  Patrick Eriksson <Patrick.Eriksson@chalmers.se>
3 
4  This program is free software; you can redistribute it and/or modify it
5  under the terms of the GNU General Public License as published by the
6  Free Software Foundation; either version 2, or (at your option) any
7  later version.
8 
9  This program is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  GNU General Public License for more details.
13 
14  You should have received a copy of the GNU General Public License
15  along with this program; if not, write to the Free Software
16  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17  USA. */
18 
19 /*===========================================================================
20  === File description
21  ===========================================================================*/
22 
31 /*===========================================================================
32  === External declarations
33  ===========================================================================*/
34 
35 #include "surface.h"
36 #include <cmath>
37 #include "auto_md.h"
38 #include "check_input.h"
39 #include "complex.h"
40 #include "geodetic.h"
41 #include "math_funcs.h"
42 #include "matpackI.h"
43 #include "physics_funcs.h"
44 #include "workspace_ng.h"
45 
46 /*===========================================================================
47  === The functions (in alphabetical order)
48  ===========================================================================*/
49 
51  return (180 - abs(rte_los[0]) + abs(specular_los[0])) / 2;
52 }
53 
55  ConstVectorView z_profile) {
56  Index ip = 0;
57  while (z_surface >= z_profile[ip+1]) {
58  ip++;
59  }
60  return ip;
61 }
62 
68  // Some sizes
69  const Index nf = I.nrows();
70  const Index stokes_dim = I.ncols();
71  const Index nlos = surface_los.nrows();
72 
74 
75  // Loop *surface_los*-es. If no such LOS, we are ready.
76  if (nlos > 0) {
77  for (Index ilos = 0; ilos < nlos; ilos++) {
78  Vector rtmp(stokes_dim); // Reflected Stokes vector for 1 frequency
79 
80  for (Index iv = 0; iv < nf; iv++) {
81  mult(rtmp, surface_rmatrix(ilos, iv, joker, joker), I(ilos, iv, joker));
82  iy(iv, joker) += rtmp;
83  }
84  }
85  }
86 }
87 
90  const Complex& Rv,
91  const Complex& Rh,
92  const Numeric& f,
93  const Index& stokes_dim,
94  const Numeric& surface_skin_t) {
95  assert(surface_rmatrix.nrows() == stokes_dim);
96  assert(surface_rmatrix.ncols() == stokes_dim);
97  assert(surface_emission.nelem() == stokes_dim);
98 
99  // Expressions are derived in the surface chapter in the user guide
100 
101  surface_rmatrix = 0.0;
102  surface_emission = 0.0;
103 
104  Numeric B = planck(f, surface_skin_t);
105 
106  const Numeric rv = pow(abs(Rv), 2.0);
107  const Numeric rh = pow(abs(Rh), 2.0);
108  const Numeric rmean = (rv + rh) / 2;
109 
110  surface_rmatrix(0, 0) = rmean;
111  surface_emission[0] = B * (1 - rmean);
112 
113  if (stokes_dim > 1) {
114  const Numeric rdiff = (rv - rh) / 2;
115 
116  surface_rmatrix(1, 0) = rdiff;
117  surface_rmatrix(0, 1) = rdiff;
118  surface_rmatrix(1, 1) = rmean;
119  surface_emission[1] = -B * rdiff;
120 
121  if (stokes_dim > 2) {
122  const Complex a = Rh * conj(Rv);
123  const Complex b = Rv * conj(Rh);
124  const Numeric c = real(a + b) / 2.0;
125 
126  surface_rmatrix(2, 2) = c;
127 
128  if (stokes_dim > 3) {
129  const Numeric d = imag(a - b) / 2.0;
130 
131  surface_rmatrix(2, 3) = d;
132  surface_rmatrix(3, 2) = -d;
133  surface_rmatrix(3, 3) = c;
134  }
135  }
136  }
137 }
138 
140  const Vector& lat_grid,
141  const Vector& lon_grid,
144  // Check sizes
145  if (surface_props_data.npages() != surface_props_names.nelem())
146  throw runtime_error(
147  "The number of pages in *surface_props_data* and "
148  "length of *surface_props_names* differ.");
149  // If no surface properties, then we are ready
150  if (surface_props_names.nelem() == 0) {
151  return;
152  }
153  if (surface_props_data.nrows() !=
154  (atmosphere_dim == 1 ? 1 : lat_grid.nelem()))
155  throw runtime_error("Row-size of *surface_props_data* not as expected.");
156  if (surface_props_data.ncols() !=
157  (atmosphere_dim <= 2 ? 1 : lon_grid.nelem()))
158  throw runtime_error("Column-size of *surface_props_data* not as expected.");
159 
160  for (Index i = 0; i < surface_props_names.nelem(); i++) {
161  if (surface_props_names[i].nelem() == 0) {
162  ostringstream os;
163  os << "Element " << i << " (0-based) of *surface_props_names* is empty.";
164  throw runtime_error(os.str());
165  }
166  for (Index j = i + 1; j < surface_props_names.nelem(); j++) {
168  ostringstream os;
169  os << "Two surface properties with same name found!\n"
170  << "This found for these two properties\n"
171  << " index: " << i << endl
172  << " index: " << j << endl
173  << " name: " << surface_props_names[i];
174  throw runtime_error(os.str());
175  }
176  }
177  }
178 }
179 
181  const String& vname,
182  const Index& atmosphere_dim,
183  const ArrayOfGridPos& gp_lat,
184  const ArrayOfGridPos& gp_lon,
185  const Matrix& itw,
188  assert(v.nelem() == 1);
189  assert(surface_props_data.npages() == surface_props_names.nelem());
190 
191  for (Index i = 0; i < surface_props_names.nelem(); i++) {
192  if (surface_props_names[i] == vname) {
196  gp_lat,
197  gp_lon,
198  itw);
199  return;
200  }
201  }
202 
203  ostringstream os;
204  os << "The following property was requested\n"
205  << " " << vname << endl
206  << "but it could not be found in *surface_props_names*.";
207  throw runtime_error(os.str());
208 }
209 
214  const Index nq = dsurface_names.nelem();
215 
216  if (dsurface_rmatrix_dx.nelem() != nq) {
217  throw runtime_error(
218  "The lengths of *dsurface_names* and *dsurface_rmatrix_dx* differ.");
219  }
220  if (dsurface_emission_dx.nelem() != nq) {
221  throw runtime_error(
222  "The lengths of *dsurface_names* and *dsurface_emission_dx* differ.");
223  }
224 
225  for (Index i = 0; i < nq; i++) {
226  bool found = false;
227  for (Index j = 0; j < surface_props_names.nelem() && !found; j++) {
228  if (dsurface_names[i] == surface_props_names[j]) {
229  found = true;
230  }
231  }
232  if (!found) {
233  ostringstream os;
234  os << "String " << i << " (0-based) of *dsurface_names* is \""
235  << dsurface_names[i] << "\"\n"
236  << "but this string could not be found in *surface_props_names*.\n"
237  << "This is likely due to incorrect choice of quantity when\n"
238  << " calling *jacobianAddSurfaceQuantity*.";
239  throw runtime_error(os.str());
240  }
241  }
242 }
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Var::atmosphere_dim
Index atmosphere_dim(Workspace &ws) noexcept
Definition: autoarts.h:2510
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
auto_md.h
conj
constexpr Complex conj(Complex c)
the conjugate of c
Definition: complex.h:65
calc_incang
Numeric calc_incang(ConstVectorView rte_los, ConstVectorView specular_los)
Calculates the incidence angle for a flat surface, based on rte_los and specular_los.
Definition: surface.cc:50
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
surface.h
ARTS::Var::lat_grid
Vector lat_grid(Workspace &ws) noexcept
Definition: autoarts.h:3962
joker
const Joker joker
ARTS::Var::surface_rmatrix
Tensor4 surface_rmatrix(Workspace &ws) noexcept
Definition: autoarts.h:6800
surface_specular_R_and_b
void surface_specular_R_and_b(MatrixView surface_rmatrix, VectorView surface_emission, const Complex &Rv, const Complex &Rh, const Numeric &f, const Index &stokes_dim, const Numeric &surface_skin_t)
Sets up the surface reflection matrix and emission vector for the case of specular reflection.
Definition: surface.cc:88
surface_props_check
void surface_props_check(const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names)
Peforms basic checks of surface_props_data and surface_props_names
Definition: surface.cc:139
abs
#define abs(x)
Definition: legacy_continua.cc:20626
mult
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
surface_calc
void surface_calc(Matrix &iy, ConstTensor3View I, ConstMatrixView surface_los, ConstTensor4View surface_rmatrix, ConstMatrixView surface_emission)
Weights together downwelling radiation and surface emission.
Definition: surface.cc:63
ARTS::Var::specular_los
Vector specular_los(Workspace &ws) noexcept
Definition: autoarts.h:6587
dsurface_check
void dsurface_check(const ArrayOfString &surface_props_names, const ArrayOfString &dsurface_names, const ArrayOfTensor4 dsurface_rmatrix_dx, const ArrayOfMatrix &dsurface_emission_dx)
Peforms basic checks of the dsurface variables.
Definition: surface.cc:210
ARTS::Var::stokes_dim
Index stokes_dim(Workspace &ws) noexcept
Definition: autoarts.h:6650
index_of_zsurface
Index index_of_zsurface(const Numeric &z_surface, ConstVectorView z_profile)
Lccates the surface with respect to pressure levels.
Definition: surface.cc:54
ARTS::Var::iy
Matrix iy(Workspace &ws) noexcept
Definition: autoarts.h:3690
Complex
std::complex< Numeric > Complex
Definition: complex.h:33
ARTS::Var::dsurface_names
ArrayOfString dsurface_names(Workspace &ws) noexcept
Definition: autoarts.h:3344
ConstTensor4View
A constant view of a Tensor4.
Definition: matpackIV.h:133
matpackI.h
Implementation of Matrix, Vector, and such stuff.
complex.h
A class implementing complex numbers for ARTS.
Array< String >
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
Absorption::nelem
Index nelem(const Lines &l)
Number of lines.
Definition: absorptionlines.h:1820
my_basic_string< char >
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
VectorView
The VectorView class.
Definition: matpackI.h:610
ARTS::Var::rte_los
Vector rte_los(Workspace &ws) noexcept
Definition: autoarts.h:5651
physics_funcs.h
This file contains declerations of functions of physical character.
surface_props_interp
void surface_props_interp(Vector &v, const String &vname, const Index &atmosphere_dim, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, const Matrix &itw, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names)
Peforms an interpolation of surface_props_data
Definition: surface.cc:180
ARTS::Var::dsurface_emission_dx
ArrayOfMatrix dsurface_emission_dx(Workspace &ws) noexcept
Definition: autoarts.h:3333
ARTS::Var::surface_props_data
Tensor3 surface_props_data(Workspace &ws) noexcept
Definition: autoarts.h:6745
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
math_funcs.h
geodetic.h
ARTS::Var::surface_emission
Matrix surface_emission(Workspace &ws) noexcept
Definition: autoarts.h:6700
ConstTensor3View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
ARTS::Var::lon_grid
Vector lon_grid(Workspace &ws) noexcept
Definition: autoarts.h:4090
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
planck
Numeric planck(const Numeric &f, const Numeric &t)
planck
Definition: physics_funcs.cc:269
workspace_ng.h
This file contains the Workspace class.
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:132
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
ARTS::Var::surface_skin_t
Numeric surface_skin_t(Workspace &ws) noexcept
Definition: autoarts.h:6877
interp_atmsurface_by_itw
void interp_atmsurface_by_itw(VectorView x, const Index &atmosphere_dim, ConstMatrixView x_surface, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
Interpolates a surface-type variable with pre-calculated weights by interp_atmsurface_gp2itw.
Definition: special_interp.cc:247
check_input.h
ARTS::Var::surface_los
Matrix surface_los(Workspace &ws) noexcept
Definition: autoarts.h:6714
Vector
The Vector class.
Definition: matpackI.h:860
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
ARTS::Var::dsurface_rmatrix_dx
ArrayOfTensor4 dsurface_rmatrix_dx(Workspace &ws) noexcept
Definition: autoarts.h:3356
ARTS::Var::z_surface
Matrix z_surface(Workspace &ws) noexcept
Definition: autoarts.h:7754
ARTS::Var::surface_props_names
ArrayOfString surface_props_names(Workspace &ws) noexcept
Definition: autoarts.h:6759