ARTS  2.4.0(git:4fb77825)
mc_antenna.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 
30 #include "mc_antenna.h"
31 #include <cfloat>
32 #include <sstream>
33 
34 extern const Numeric PI;
35 extern const Numeric DEG2RAD;
36 extern const Numeric RAD2DEG;
37 
38 Numeric ran_gaussian(Rng& rng, const Numeric sigma) {
39  Numeric x, y, r2;
40 
41  do {
42  /* choose x,y in uniform square (-1,-1) to (+1,+1) */
43 
44  x = -1 + 2 * rng.draw();
45  y = -1 + 2 * rng.draw();
46 
47  /* see if it is in the unit circle */
48  r2 = x * x + y * y;
49  } while (r2 > 1.0 || r2 == 0);
50 
51  /* Box-Muller transform */
52  return sigma * y * sqrt(-2.0 * log(r2) / r2);
53 }
54 
56  Numeric phi;
57 
58  phi = 360.0 * rng.draw() - 180.0;
59 
60  return phi;
61 }
62 
63 void rotmat_enu(MatrixView R_ant2enu, ConstVectorView prop_los)
64 
65 {
66  Numeric cza, sza, caa, saa;
67 
68  cza = cos(prop_los[0] * DEG2RAD);
69  sza = sin(prop_los[0] * DEG2RAD);
70  caa = cos(prop_los[1] * DEG2RAD);
71  saa = sin(prop_los[1] * DEG2RAD);
72 
73  R_ant2enu(0, 0) = -cza * saa;
74  R_ant2enu(0, 1) = caa;
75  R_ant2enu(0, 2) = sza * saa;
76 
77  R_ant2enu(1, 0) = -cza * caa;
78  R_ant2enu(1, 1) = -saa;
79  R_ant2enu(1, 2) = sza * caa;
80 
81  R_ant2enu(2, 0) = sza;
82  R_ant2enu(2, 1) = 0.0;
83  R_ant2enu(2, 2) = cza;
84 }
85 
87  const Index& stokes_dim,
88  const Numeric& f1_dir,
89  const Numeric& f2_dir,
90  ConstMatrixView R_f1,
91  ConstMatrixView R_f2)
92 
93 {
94  const Numeric flip = f1_dir * f2_dir;
95  Numeric cos_pra1, sin_pra1, cos_pra2, sin_pra2;
96 
97  cos_pra1 = R_f1(joker, 0) * R_f2(joker, 0);
98  sin_pra1 = f2_dir * (R_f1(joker, 0) * R_f2(joker, 1));
99  sin_pra2 = f1_dir * (R_f1(joker, 1) * R_f2(joker, 0));
100  cos_pra2 = f1_dir * f2_dir * (R_f1(joker, 1) * R_f2(joker, 1));
101 
102  R_pra = 0.0;
103  R_pra(0, 0) = 1.0;
104  if (stokes_dim > 1) {
105  R_pra(1, 1) = 2 * cos_pra1 * cos_pra1 - 1.0;
106  if (stokes_dim > 2) {
107  R_pra(1, 2) = flip * 2 * cos_pra1 * sin_pra1;
108  R_pra(2, 1) = 2 * cos_pra2 * sin_pra2;
109  R_pra(2, 2) = flip * (2 * cos_pra2 * cos_pra2 - 1.0);
110  if (stokes_dim > 3) {
111  R_pra(3, 3) = flip * 1.0;
112  }
113  }
114  }
115 }
116 
118 
119 void MCAntenna::set_gaussian(const Numeric& za_sigma, const Numeric& aa_sigma) {
121  sigma_za = za_sigma;
122  sigma_aa = aa_sigma;
123 }
124 
126  const Numeric& aa_fwhm) {
128  sigma_za = za_fwhm / 2.3548;
129  sigma_aa = aa_fwhm / 2.3548;
130 }
131 
133  ConstVectorView aa_grid_,
134  ConstMatrixView G_lookup_) {
136  za_grid = za_grid_;
137  aa_grid = aa_grid_;
138  G_lookup = G_lookup_;
139 }
140 
142 
144  ConstMatrixView R_return,
145  ConstMatrixView R_enu2ant) const {
146  Numeric z, term_el, term_az;
147  Numeric ant_el, ant_az;
148  Vector k_vhk(3);
149 
150  switch (atype) {
152  wgt = 1.0;
153  break;
154 
156 
157  mult(k_vhk, R_enu2ant, R_return(joker, 2));
158 
159  // Assume Gaussian is narrow enough that response is 0 beyond 90 degrees
160  // Same assumption is made for drawing samples (draw_los)
161  if (k_vhk[2] > 0) {
162  ant_el = atan(k_vhk[0] / k_vhk[2]) * RAD2DEG;
163  ant_az = atan(k_vhk[1] / k_vhk[2]) * RAD2DEG;
164  term_el = ant_el / sigma_za;
165  term_az = ant_az / sigma_aa;
166  z = term_el * term_el + term_az * term_az;
167  wgt = exp(-0.5 * z);
168  } else {
169  wgt = 0.0;
170  }
171  break;
172 
173  default:
174  ostringstream os;
175  os << "invalid Antenna type.";
176  throw runtime_error(os.str());
177  }
178 }
179 
180 void MCAntenna::draw_los(VectorView sampled_rte_los,
181  MatrixView R_los,
182  Rng& rng,
183  ConstMatrixView R_ant2enu,
184  ConstVectorView bore_sight_los) const {
185  Numeric ant_el, ant_az, ant_r;
186  Numeric tel, taz;
187  Vector k_vhk(3);
188 
189  switch (atype) {
191  sampled_rte_los = bore_sight_los;
192  R_los = R_ant2enu;
193  break;
194 
196 
197  ant_el = 91;
198  ant_az = 91;
199 
200  // Assume Gaussian is narrow enough that response is 0 beyond 90 degrees
201  // Same assumption is made for radar return samples (return_los)
202  while (ant_el >= 90) {
203  ant_el = ran_gaussian(rng, sigma_za);
204  }
205  while (ant_az >= 90) {
206  ant_az = ran_gaussian(rng, sigma_aa);
207  }
208 
209  // Propagation direction
210  tel = tan(DEG2RAD * ant_el);
211  taz = tan(DEG2RAD * ant_az);
212  ant_r = sqrt(1 + tel * tel + taz * taz);
213  k_vhk[0] = tel / ant_r;
214  k_vhk[1] = taz / ant_r;
215  k_vhk[2] = (Numeric)1.0 / ant_r;
216  mult(R_los(joker, 2), R_ant2enu, k_vhk);
217 
218  sampled_rte_los[0] = acos(R_los(2, 2)) * RAD2DEG;
219 
220  // Horizontal polarization basis
221  // If drawn los is at zenith or nadir, assume same azimuth as boresight
222  if (((Numeric)1.0 - abs(R_los(2, 2))) < DBL_EPSILON) {
223  // H is aligned with H of bs, use row not column because tranpose
224  R_los(joker, 1) = R_ant2enu(1, joker);
225  sampled_rte_los[1] = bore_sight_los[1];
226  } else {
227  const Vector uhat{0.0, 0.0, 1.0};
228  Numeric magh;
229  sampled_rte_los[1] = atan2(R_los(0, 2), R_los(1, 2)) * RAD2DEG;
230  cross3(R_los(joker, 1), R_los(joker, 2), uhat);
231  magh = sqrt(R_los(joker, 1) * R_los(joker, 1));
232  R_los(joker, 1) /= magh;
233  }
234 
235  // Vertical polarization basis
236  cross3(R_los(joker, 0), R_los(joker, 1), R_los(joker, 2));
237 
238  break;
239 
240  default:
241  ostringstream os;
242  os << "invalid Antenna type.";
243  throw runtime_error(os.str());
244  }
245 }
246 
247 ostream& operator<<(ostream& os, const MCAntenna&) {
248  os << "MCAntenna: Output operator not implemented";
249  return os;
250 }
MCAntenna::get_type
AntennaType get_type() const
AntennaType get_type.
Definition: mc_antenna.cc:141
ran_gaussian
Numeric ran_gaussian(Rng &rng, const Numeric sigma)
ran_gaussian.
Definition: mc_antenna.cc:38
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
MCAntenna::aa_grid
Vector aa_grid
Definition: mc_antenna.h:54
ANTENNA_TYPE_LOOKUP
@ ANTENNA_TYPE_LOOKUP
Definition: mc_antenna.h:44
joker
const Joker joker
abs
#define abs(x)
Definition: legacy_continua.cc:20626
ANTENNA_TYPE_PENCIL_BEAM
@ ANTENNA_TYPE_PENCIL_BEAM
Definition: mc_antenna.h:42
mult
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
MCAntenna::set_pencil_beam
void set_pencil_beam()
set_pencil_beam
Definition: mc_antenna.cc:117
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
ARTS::Var::stokes_dim
Index stokes_dim(Workspace &ws) noexcept
Definition: autoarts.h:6650
MCAntenna
An Antenna object used by MCGeneral.
Definition: mc_antenna.h:51
MCAntenna::set_gaussian_fwhm
void set_gaussian_fwhm(const Numeric &za_fwhm, const Numeric &aa_fwhm)
set_gaussian_fwhm.
Definition: mc_antenna.cc:125
MCAntenna::G_lookup
Matrix G_lookup
Definition: mc_antenna.h:55
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
cross3
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b)
cross3
Definition: matpackI.cc:1393
DEG2RAD
const Numeric DEG2RAD
Rng::draw
double draw()
Draws a double from the uniform distribution [0,1)
Definition: rng.cc:97
MCAntenna::draw_los
void draw_los(VectorView sampled_rte_los, MatrixView R_los, Rng &rng, ConstMatrixView R_ant2enu, ConstVectorView bore_sight_los) const
draw_los.
Definition: mc_antenna.cc:180
MCAntenna::atype
AntennaType atype
Definition: mc_antenna.h:52
MCAntenna::sigma_aa
Numeric sigma_aa
Definition: mc_antenna.h:53
rotmat_stokes
void rotmat_stokes(MatrixView R_pra, const Index &stokes_dim, const Numeric &f1_dir, const Numeric &f2_dir, ConstMatrixView R_f1, ConstMatrixView R_f2)
rotmat_stokes.
Definition: mc_antenna.cc:86
operator<<
ostream & operator<<(ostream &os, const MCAntenna &)
Definition: mc_antenna.cc:247
VectorView
The VectorView class.
Definition: matpackI.h:610
MCAntenna::sigma_za
Numeric sigma_za
Definition: mc_antenna.h:53
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
mc_antenna.h
Workspace functions for the solution of cloud-box radiative transfer by Monte Carlo methods....
RAD2DEG
const Numeric RAD2DEG
MCAntenna::za_grid
Vector za_grid
Definition: mc_antenna.h:54
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
Rng
Definition: rng.h:554
AntennaType
AntennaType
Definition: mc_antenna.h:41
MCAntenna::set_lookup
void set_lookup(ConstVectorView za_grid, ConstVectorView aa_grid, ConstMatrixView G_lookup)
set_lookup.
Definition: mc_antenna.cc:132
ran_uniform
Numeric ran_uniform(Rng &rng)
ran_uniform.
Definition: mc_antenna.cc:55
rotmat_enu
void rotmat_enu(MatrixView R_ant2enu, ConstVectorView prop_los)
rotmat_enu.
Definition: mc_antenna.cc:63
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
PI
const Numeric PI
Vector
The Vector class.
Definition: matpackI.h:860
MCAntenna::return_los
void return_los(Numeric &wgt, ConstMatrixView R_return, ConstMatrixView R_enu2ant) const
return_los
Definition: mc_antenna.cc:143
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
MCAntenna::set_gaussian
void set_gaussian(const Numeric &za_sigma, const Numeric &aa_sigma)
set_gaussian.
Definition: mc_antenna.cc:119
ANTENNA_TYPE_GAUSSIAN
@ ANTENNA_TYPE_GAUSSIAN
Definition: mc_antenna.h:43