Line data Source code
1 : /**
2 : \file radiometermeasure.cpp
3 : Bojan Nikolic <bn204@mrao.cam.ac.uk>, <bojan@bnikolic.co.uk>
4 :
5 : Initial version February 2008
6 : Maintained by ESO since 2013.
7 : Renamed radiometermeasure.cc 2023
8 :
9 : */
10 :
11 : #include <functional>
12 : #include <vector>
13 :
14 : #include "radiometermeasure.h"
15 : #include "radiometer_utils.h"
16 :
17 : namespace LibAIR2 {
18 :
19 0 : Radiometer::Radiometer( const std::vector<double> & FGrid,
20 0 : const std::vector<double> & coeffs):
21 0 : FGrid(FGrid),
22 0 : coeffs_v(1)
23 : {
24 0 : coeffs_v[0]=coeffs;
25 0 : }
26 :
27 0 : Radiometer::Radiometer( const std::vector<double> & FGrid,
28 0 : const std::vector< std::vector<double> > & coeffs_v):
29 0 : FGrid(FGrid),
30 0 : coeffs_v(coeffs_v)
31 : {
32 0 : }
33 :
34 :
35 0 : DSBRadio::DSBRadio( double f_0,
36 0 : double f_if):
37 0 : f_0(f_0),
38 0 : f_if(f_if)
39 : {
40 0 : std::vector<double> FGrid(2);
41 0 : FGrid[0] = f_0- f_if;
42 0 : FGrid[1] = f_0+ f_if;
43 0 : std::vector<double> coeffs(2, 0.5);
44 :
45 0 : r.reset( new Radiometer(FGrid, coeffs));
46 0 : }
47 :
48 0 : DSBRadio::DSBRadio( double f_0,
49 : double f_if,
50 0 : Radiometer * rr):
51 0 : r(rr),
52 0 : f_0(f_0),
53 0 : f_if(f_if)
54 : {
55 :
56 0 : }
57 :
58 0 : const Radiometer & DSBRadio::getRadiometer(void) const
59 : {
60 0 : return *r;
61 : }
62 :
63 :
64 0 : double DSBRadio::eval( const std::vector<double> & skyTb) const
65 : {
66 0 : return r->eval(skyTb);
67 : }
68 :
69 0 : DSBBWRadio::DSBBWRadio( double f_0,
70 : double f_if,
71 0 : double f_bw):
72 : DSBRadio(f_0, f_if,
73 : MkRadio( f_0, f_if, f_bw) ),
74 0 : f_bw(f_bw)
75 : {
76 0 : }
77 :
78 0 : Radiometer * DSBBWRadio::MkRadio( double f_0,
79 : double f_if,
80 : double f_bw)
81 : {
82 0 : const size_t nsamples=30;
83 :
84 0 : std::vector<double> FGrid(nsamples *2 );
85 0 : std::vector<double> coeffs(nsamples *2, 1.0/(nsamples *2));
86 :
87 0 : const double delta = f_bw / (nsamples -1);
88 0 : const double flow = f_if - f_bw *0.5;
89 :
90 0 : for (size_t i =0 ; i < nsamples ; ++i)
91 : {
92 0 : const double fif_c = flow + delta * i;
93 0 : FGrid[nsamples-i-1] = f_0- fif_c;
94 0 : FGrid[nsamples+i] = f_0+ fif_c;
95 :
96 : }
97 0 : return new Radiometer(FGrid, coeffs);
98 :
99 0 : }
100 :
101 0 : DSBBW_QuadRadio::DSBBW_QuadRadio( double f_0,
102 : double f_if,
103 0 : double f_bw):
104 : DSBRadio(f_0, f_if,
105 0 : MkRadio( f_0, f_if, f_bw) )
106 : {
107 0 : }
108 :
109 0 : Radiometer * DSBBW_QuadRadio::MkRadio( double f_0,
110 : double f_if,
111 : double f_bw)
112 : {
113 : // 5 point Gauss-Legendre Rule.
114 : // See http://www.sitmo.com/eq/423
115 0 : const double x[] = {-9.06179845938663992811e-01,
116 : -5.38469310105683091018e-01,
117 : 0.00000000000000000000e+00,
118 : 5.38469310105683091018e-01,
119 : 9.06179845938663992811e-01};
120 :
121 0 : const double w[] = {2.36926885056189087515e-01,
122 : 4.78628670499366468030e-01,
123 : 5.68888888888888888883e-01,
124 : 4.78628670499366468030e-01,
125 : 2.36926885056189087515e-01};
126 :
127 0 : const size_t nsamples = sizeof(x)/sizeof(double);
128 :
129 0 : std::vector<double> FGrid(nsamples*2);
130 0 : std::vector<double> coeffs(nsamples*2);
131 :
132 : /// Note we want the average signal over the bandwidth. Therefore
133 : /// integrathe the function over bandwidth then divide by
134 : /// bandwidth to get average value. As ressult coefficients below
135 : /// do not f_bw in their expressions.
136 :
137 : /// Note that there is an extra factor of * 0.5 becuase we have
138 : /// two bands
139 :
140 : /// Lower band
141 0 : for (size_t i =0 ; i < nsamples ; ++i)
142 : {
143 0 : const double fc = f_0 - f_if;
144 0 : FGrid[i] = fc + 0.5* f_bw * x[i];
145 0 : coeffs[i]= 0.5 * 0.5 * w[i];
146 : }
147 :
148 : /// Upper band
149 0 : for (size_t i =0 ; i < nsamples ; ++i)
150 : {
151 0 : const double fc = f_0 + f_if;
152 0 : FGrid [i+nsamples] = fc + 0.5* f_bw * x[i];
153 0 : coeffs[i+nsamples]= 0.5 * 0.5 * w[i];
154 : }
155 0 : return new Radiometer(FGrid, coeffs);
156 0 : }
157 :
158 0 : InvalidWVRChannel::InvalidWVRChannel(int chlow,
159 : int chhigh,
160 0 : int ch) :
161 0 : chlow(chlow),
162 0 : chhigh(chhigh),
163 0 : ch(ch)
164 : {
165 0 : }
166 :
167 0 : DSBRadio * MkALMARadiometer(int ch,
168 : double cf_off,
169 : double bw_off) noexcept(false)
170 : {
171 0 : const double filter_c[] = { 1.25, 3.25, 5.5, 7.25};
172 0 : const double filter_bw[]= { 1.5 , 2.5 , 2.0, 1.5};
173 :
174 0 : if ( ch<1 or ch > 4)
175 : {
176 0 : throw InvalidWVRChannel(1 , 4, ch);
177 : }
178 :
179 : return new DSBBW_QuadRadio( 183.31,
180 0 : filter_c[ch-1] +cf_off,
181 0 : filter_bw[ch-1]+bw_off);
182 :
183 :
184 : }
185 :
186 0 : DSBRadio * MkALMADickeProto(int ch) noexcept(false)
187 : {
188 :
189 0 : const double filter_c[] = { 0.88 , 1.94 , 3.175, 5.2 };
190 0 : const double filter_bw[]= { 0.16 , 0.75 , 1.25 , 2.5 };
191 :
192 0 : if ( ch<1 or ch > 4)
193 : {
194 0 : throw InvalidWVRChannel(1 , 4, ch);
195 : }
196 :
197 : return new DSBBW_QuadRadio( 183.31,
198 0 : filter_c[ch-1],
199 0 : filter_bw[ch-1]);
200 :
201 :
202 : }
203 :
204 : /***
205 : Merge the four channels into single radiomer
206 : */
207 0 : Radiometer * MkFullWVR( const std::function< DSBRadio* (int ch) > & mkchn )
208 : {
209 0 : std::vector< std::shared_ptr<DSBRadio> > dsbv;
210 0 : std::vector<const Radiometer *> rv;
211 :
212 0 : for (size_t j = 1 ; j < 5; ++j )
213 : {
214 0 : dsbv.push_back( std::shared_ptr<DSBRadio> (mkchn(j)) );
215 0 : rv.push_back( & dsbv[ j-1]->getRadiometer());
216 : }
217 :
218 0 : return MergeRadiometers( rv).release();
219 0 : }
220 :
221 0 : Radiometer * MkFullALMAWVR(void)
222 : {
223 0 : return MkFullWVR(std::bind(MkALMARadiometer,std::placeholders::_1,0.0,0.0));
224 : }
225 :
226 0 : Radiometer * MkALMAWVR_offset(double cf,
227 : double bw)
228 : {
229 0 : return MkFullWVR(std::bind(MkALMARadiometer,
230 : std::placeholders::_1,
231 : cf,
232 0 : bw));
233 : }
234 :
235 0 : ALMAWVRCharacter::ALMAWVRCharacter(void)
236 : {
237 0 : cf1=1.25;
238 0 : cf2=3.25;
239 0 : cf3=5.5;
240 0 : cf4=7.25;
241 :
242 0 : bw1=1.5;
243 0 : bw2=2.5;
244 0 : bw3=2.0;
245 0 : bw4=1.5;
246 :
247 0 : }
248 :
249 0 : Radiometer *MkALMAWVR(const ALMAWVRCharacter &c)
250 : {
251 0 : std::vector< std::shared_ptr<DSBRadio> > dsbv;
252 0 : dsbv.push_back(std::shared_ptr<DSBRadio>(new DSBBW_QuadRadio(183.31,
253 0 : c.cf1,
254 0 : c.bw1))
255 : );
256 0 : dsbv.push_back(std::shared_ptr<DSBRadio>(new DSBBW_QuadRadio(183.31,
257 0 : c.cf2,
258 0 : c.bw2))
259 : );
260 0 : dsbv.push_back(std::shared_ptr<DSBRadio>(new DSBBW_QuadRadio(183.31,
261 0 : c.cf3,
262 0 : c.bw3))
263 : );
264 0 : dsbv.push_back(std::shared_ptr<DSBRadio>(new DSBBW_QuadRadio(183.31,
265 0 : c.cf4,
266 0 : c.bw4))
267 : );
268 :
269 0 : std::vector<const Radiometer *> rv;
270 0 : for(size_t i=0; i<dsbv.size(); ++i)
271 : {
272 0 : rv.push_back(&dsbv[i]->getRadiometer());
273 : }
274 0 : return MergeRadiometers( rv).release();
275 0 : }
276 :
277 :
278 0 : Radiometer * MkFullDickeProtoWVR(void)
279 : {
280 0 : return MkFullWVR( MkALMADickeProto);
281 : }
282 :
283 0 : Radiometer * MkSSBRadio(double f_c,
284 : double f_bw)
285 : {
286 : // 5 point Gauss-Legendre Rule. See http://www.sitmo.com/eq/423
287 0 : const double x[] = {-9.06179845938663992811e-01,
288 : -5.38469310105683091018e-01,
289 : 0.00000000000000000000e+00,
290 : 5.38469310105683091018e-01,
291 : 9.06179845938663992811e-01};
292 :
293 0 : const double w[] = {2.36926885056189087515e-01,
294 : 4.78628670499366468030e-01,
295 : 5.68888888888888888883e-01,
296 : 4.78628670499366468030e-01,
297 : 2.36926885056189087515e-01};
298 :
299 0 : const size_t nsamples = sizeof(x)/sizeof(double);
300 :
301 0 : std::vector<double> FGrid(nsamples);
302 0 : std::vector<double> coeffs(nsamples);
303 :
304 0 : for (size_t i =0 ; i < nsamples ; ++i)
305 : {
306 0 : FGrid[i] = f_c + 0.5* f_bw * x[i];
307 0 : coeffs[i]= 0.5 * w[i];
308 : }
309 :
310 0 : return new Radiometer(FGrid, coeffs);
311 :
312 0 : }
313 :
314 0 : Radiometer * MkIRAM22(void)
315 : {
316 0 : const double filter_c[] = { 19.2, 22, 25.2};
317 0 : const double filter_bw[] = {1.15, 1.15, 1.15};
318 :
319 0 : std::vector< std::shared_ptr<Radiometer> > rv;
320 0 : std::vector< const Radiometer * > rv_obs;
321 :
322 0 : for (size_t i =0 ; i < 3; ++i)
323 : {
324 0 : rv.push_back(std::shared_ptr<Radiometer>(MkSSBRadio(filter_c[i],
325 0 : filter_bw[i])));
326 :
327 0 : rv_obs.push_back(rv[i].get());
328 : }
329 :
330 0 : return MergeRadiometers(rv_obs).release();
331 :
332 0 : }
333 :
334 : }
335 :
336 :
337 :
|
