ARTS 2.5.9 (git: 825fa5f2)
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
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#include "arts_constants.h"
34#include "arts_conversions.h"
35
36inline constexpr Numeric PI=Constant::pi;
39
40Numeric ran_gaussian(Rng& rng, const Numeric sigma) {
41 Numeric x, y, r2;
42
43 do {
44 /* choose x,y in uniform square (-1,-1) to (+1,+1) */
45
46 x = -1 + 2 * rng.draw();
47 y = -1 + 2 * rng.draw();
48
49 /* see if it is in the unit circle */
50 r2 = x * x + y * y;
51 } while (r2 > 1.0 || r2 == 0);
52
53 /* Box-Muller transform */
54 return sigma * y * sqrt(-2.0 * log(r2) / r2);
55}
56
58 Numeric phi;
59
60 phi = 360.0 * rng.draw() - 180.0;
61
62 return phi;
63}
64
65void rotmat_enu(MatrixView R_ant2enu, ConstVectorView prop_los)
66
67{
68 Numeric cza, sza, caa, saa;
69
70 cza = cos(prop_los[0] * DEG2RAD);
71 sza = sin(prop_los[0] * DEG2RAD);
72 caa = cos(prop_los[1] * DEG2RAD);
73 saa = sin(prop_los[1] * DEG2RAD);
74
75 R_ant2enu(0, 0) = -cza * saa;
76 R_ant2enu(0, 1) = caa;
77 R_ant2enu(0, 2) = sza * saa;
78
79 R_ant2enu(1, 0) = -cza * caa;
80 R_ant2enu(1, 1) = -saa;
81 R_ant2enu(1, 2) = sza * caa;
82
83 R_ant2enu(2, 0) = sza;
84 R_ant2enu(2, 1) = 0.0;
85 R_ant2enu(2, 2) = cza;
86}
87
89 const Index& stokes_dim,
90 const Numeric& f1_dir,
91 const Numeric& f2_dir,
92 ConstMatrixView R_f1,
93 ConstMatrixView R_f2)
94
95{
96 const Numeric flip = f1_dir * f2_dir;
97 Numeric cos_pra1, sin_pra1, cos_pra2, sin_pra2;
98
99 cos_pra1 = R_f1(joker, 0) * R_f2(joker, 0);
100 sin_pra1 = f2_dir * (R_f1(joker, 0) * R_f2(joker, 1));
101 sin_pra2 = f1_dir * (R_f1(joker, 1) * R_f2(joker, 0));
102 cos_pra2 = f1_dir * f2_dir * (R_f1(joker, 1) * R_f2(joker, 1));
103
104 R_pra = 0.0;
105 R_pra(0, 0) = 1.0;
106 if (stokes_dim > 1) {
107 R_pra(1, 1) = 2 * cos_pra1 * cos_pra1 - 1.0;
108 if (stokes_dim > 2) {
109 R_pra(1, 2) = flip * 2 * cos_pra1 * sin_pra1;
110 R_pra(2, 1) = 2 * cos_pra2 * sin_pra2;
111 R_pra(2, 2) = flip * (2 * cos_pra2 * cos_pra2 - 1.0);
112 if (stokes_dim > 3) {
113 R_pra(3, 3) = flip * 1.0;
114 }
115 }
116 }
117}
118
120
121void MCAntenna::set_gaussian(const Numeric& za_sigma, const Numeric& aa_sigma) {
123 sigma_za = za_sigma;
124 sigma_aa = aa_sigma;
125}
126
128 const Numeric& aa_fwhm) {
130 sigma_za = za_fwhm / 2.3548;
131 sigma_aa = aa_fwhm / 2.3548;
132}
133
135 ConstVectorView aa_grid_,
136 ConstMatrixView G_lookup_) {
138 za_grid = za_grid_;
139 aa_grid = aa_grid_;
140 G_lookup = G_lookup_;
141}
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 ARTS_USER_ERROR ("invalid Antenna type.")
175 }
176}
177
178void MCAntenna::draw_los(VectorView sampled_rte_los,
179 MatrixView R_los,
180 Rng& rng,
181 ConstMatrixView R_ant2enu,
182 ConstVectorView bore_sight_los) const {
183 Numeric ant_el, ant_az, ant_r;
184 Numeric tel, taz;
185 Vector k_vhk(3);
186
187 switch (atype) {
189 sampled_rte_los = bore_sight_los;
190 R_los = R_ant2enu;
191 break;
192
194
195 ant_el = 91;
196 ant_az = 91;
197
198 // Assume Gaussian is narrow enough that response is 0 beyond 90 degrees
200 while (ant_el >= 90) {
201 ant_el = ran_gaussian(rng, sigma_za);
202 }
203 while (ant_az >= 90) {
204 ant_az = ran_gaussian(rng, sigma_aa);
205 }
206
207 // Propagation direction
208 tel = tan(DEG2RAD * ant_el);
209 taz = tan(DEG2RAD * ant_az);
210 ant_r = sqrt(1 + tel * tel + taz * taz);
211 k_vhk[0] = tel / ant_r;
212 k_vhk[1] = taz / ant_r;
213 k_vhk[2] = (Numeric)1.0 / ant_r;
214 mult(R_los(joker, 2), R_ant2enu, k_vhk);
215
216 sampled_rte_los[0] = acos(R_los(2, 2)) * RAD2DEG;
217
218 // Horizontal polarization basis
219 // If drawn los is at zenith or nadir, assume same azimuth as boresight
220 if (((Numeric)1.0 - abs(R_los(2, 2))) < DBL_EPSILON) {
221 // H is aligned with H of bs, use row not column because tranpose
222 R_los(joker, 1) = R_ant2enu(1, joker);
223 sampled_rte_los[1] = bore_sight_los[1];
224 } else {
225 const Vector uhat{0.0, 0.0, 1.0};
226 Numeric magh;
227 sampled_rte_los[1] = atan2(R_los(0, 2), R_los(1, 2)) * RAD2DEG;
228 cross3(R_los(joker, 1), R_los(joker, 2), uhat);
229 magh = sqrt(R_los(joker, 1) * R_los(joker, 1));
230 R_los(joker, 1) /= magh;
231 }
232
233 // Vertical polarization basis
234 cross3(R_los(joker, 0), R_los(joker, 1), R_los(joker, 2));
235
236 break;
237
238 default:
239 ARTS_USER_ERROR ("invalid Antenna type.")
240 }
241}
242
243ostream& operator<<(ostream& os, const MCAntenna&) {
244 os << "MCAntenna: Output operator not implemented";
245 return os;
246}
Constants of physical expressions as constexpr.
Common ARTS conversions.
A constant view of a Matrix.
Definition: matpackI.h:1065
A constant view of a Vector.
Definition: matpackI.h:521
The MatrixView class.
Definition: matpackI.h:1188
Definition: rng.h:554
double draw()
Draws a double from the uniform distribution [0,1)
Definition: rng.cc:97
The VectorView class.
Definition: matpackI.h:674
The Vector class.
Definition: matpackI.h:910
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
Definition: disort.cc:57
Definition: doit.cc:60
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
cross3
Definition: matpackI.cc:1350
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
Definition: mc_antenna.cc:37
Numeric ran_gaussian(Rng &rng, const Numeric sigma)
ran_gaussian.
Definition: mc_antenna.cc:40
Definition: mc_antenna.cc:38
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:88
ostream & operator<<(ostream &os, const MCAntenna &)
Definition: mc_antenna.cc:243
Numeric ran_uniform(Rng &rng)
ran_uniform.
Definition: mc_antenna.cc:57
constexpr Numeric PI
Definition: mc_antenna.cc:36
void rotmat_enu(MatrixView R_ant2enu, ConstVectorView prop_los)
rotmat_enu.
Definition: mc_antenna.cc:65
Workspace functions for the solution of cloud-box radiative transfer by Monte Carlo methods....
Numeric ran_gaussian(Rng &rng, const Numeric sigma)
ran_gaussian.
Definition: mc_antenna.cc:40
@ ANTENNA_TYPE_PENCIL_BEAM
Definition: mc_antenna.h:41
@ ANTENNA_TYPE_GAUSSIAN
Definition: mc_antenna.h:42
@ ANTENNA_TYPE_LOOKUP
Definition: mc_antenna.h:43
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:705
An Antenna object used by MCGeneral.
Definition: mc_antenna.h:50
Matrix G_lookup
Definition: mc_antenna.h:54
Numeric sigma_aa
Definition: mc_antenna.h:52
Numeric sigma_za
Definition: mc_antenna.h:52
AntennaType atype
Definition: mc_antenna.h:51
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:178
void return_los(Numeric &wgt, ConstMatrixView R_return, ConstMatrixView R_enu2ant) const
return_los
Definition: mc_antenna.cc:143
void set_lookup(ConstVectorView za_grid, ConstVectorView aa_grid, ConstMatrixView G_lookup)
set_lookup.
Definition: mc_antenna.cc:134
void set_gaussian(const Numeric &za_sigma, const Numeric &aa_sigma)
set_gaussian.
Definition: mc_antenna.cc:121
Vector aa_grid
Definition: mc_antenna.h:53
Vector za_grid
Definition: mc_antenna.h:53
void set_gaussian_fwhm(const Numeric &za_fwhm, const Numeric &aa_fwhm)
set_gaussian_fwhm.
Definition: mc_antenna.cc:127
void set_pencil_beam()
set_pencil_beam
Definition: mc_antenna.cc:119