ARTS 2.5.0 (git: 9ee3ac6c)
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
34extern const Numeric PI;
35extern const Numeric DEG2RAD;
36extern const Numeric RAD2DEG;
37
38Numeric 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
63void 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
119void 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 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
199 // Same assumption is made for radar return samples (return_los)
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}
A constant view of a Matrix.
Definition: matpackI.h:1014
A constant view of a Vector.
Definition: matpackI.h:489
An Antenna object used by MCGeneral.
Definition: mc_antenna.h:51
Matrix G_lookup
Definition: mc_antenna.h:55
Numeric sigma_aa
Definition: mc_antenna.h:53
Numeric sigma_za
Definition: mc_antenna.h:53
AntennaType atype
Definition: mc_antenna.h:52
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
AntennaType get_type() const
AntennaType get_type.
Definition: mc_antenna.cc:141
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:132
void set_gaussian(const Numeric &za_sigma, const Numeric &aa_sigma)
set_gaussian.
Definition: mc_antenna.cc:119
Vector aa_grid
Definition: mc_antenna.h:54
Vector za_grid
Definition: mc_antenna.h:54
void set_gaussian_fwhm(const Numeric &za_fwhm, const Numeric &aa_fwhm)
set_gaussian_fwhm.
Definition: mc_antenna.cc:125
void set_pencil_beam()
set_pencil_beam
Definition: mc_antenna.cc:117
The MatrixView class.
Definition: matpackI.h:1125
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:626
The Vector class.
Definition: matpackI.h:876
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define abs(x)
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b) ARTS_NOEXCEPT
cross3
Definition: matpackI.cc:1393
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
const Joker joker
const Numeric RAD2DEG
Numeric ran_gaussian(Rng &rng, const Numeric sigma)
ran_gaussian.
Definition: mc_antenna.cc:38
const Numeric PI
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
ostream & operator<<(ostream &os, const MCAntenna &)
Definition: mc_antenna.cc:243
Numeric ran_uniform(Rng &rng)
ran_uniform.
Definition: mc_antenna.cc:55
const Numeric DEG2RAD
void rotmat_enu(MatrixView R_ant2enu, ConstVectorView prop_los)
rotmat_enu.
Definition: mc_antenna.cc:63
Workspace functions for the solution of cloud-box radiative transfer by Monte Carlo methods....
AntennaType
Definition: mc_antenna.h:41
@ ANTENNA_TYPE_PENCIL_BEAM
Definition: mc_antenna.h:42
@ ANTENNA_TYPE_GAUSSIAN
Definition: mc_antenna.h:43
@ ANTENNA_TYPE_LOOKUP
Definition: mc_antenna.h:44
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739