source: branches/hpc34/src/STAttr.cpp@ 2623

Last change on this file since 2623 was 2163, checked in by Malte Marquarding, 14 years ago

Remove various compiler warnings
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.5 KB
Line 
1//#---------------------------------------------------------------------------
2//# STAttr.cc: A collection of attributes for different telescopes
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id:
30//#---------------------------------------------------------------------------
31
32#include <casa/aips.h>
33#include <casa/Arrays/Vector.h>
34#include <casa/Arrays/ArrayMath.h>
35#include <casa/Exceptions.h>
36#include <casa/Quanta/QC.h>
37#include <casa/Quanta/Quantum.h>
38#include <casa/Quanta/MVTime.h>
39
40#include <measures/Measures/MEpoch.h>
41
42#include <scimath/Mathematics/InterpolateArray1D.h>
43
44#include "STAttr.h"
45
46using namespace casa;
47using namespace asap;
48
49STAttr::STAttr()
50{
51 initData();
52}
53
54STAttr::STAttr(const STAttr& other):
55 Logger()
56{
57 (void) other; //suppress unused warning
58 initData(); // state just private 'static' data
59}
60
61STAttr& STAttr::operator=(const STAttr& other)
62{
63 if (this != &other) {
64 ; // state just private 'static' data
65 }
66 return *this;
67}
68
69STAttr::~STAttr()
70{;}
71
72
73Float STAttr::diameter(Instrument inst) const
74{
75 Float D = 1.0;
76 switch (inst) {
77 case ALMA:
78 {
79 D = 12.0;
80 }
81 break;
82
83 case ATMOPRA:
84 {
85 D = 22.0;
86 }
87 break;
88 case ATPKSMB:
89 case ATPKSHOH:
90 {
91 D = 64.0;
92 }
93 break;
94 case TIDBINBILLA:
95 {
96 D = 70.0;
97 }
98 break;
99 case CEDUNA:
100 {
101 D = 30.0;
102 }
103 break;
104 case GBT:
105 {
106 D = 104.9;
107 }
108 break;
109 case HOBART:
110 {
111 D = 26.0;
112 }
113 break;
114 default:
115 {
116 throw(AipsError("Unknown instrument"));
117 }
118 }
119 return D;
120}
121
122Vector<Float> STAttr::beamEfficiency(Instrument inst, const MEpoch& dateObs,
123 const Vector<Float>& freqs) const
124{
125
126 // Look at date where appropriate
127 MVTime t(dateObs.getValue());
128 uInt year = t.year();
129
130 Vector<Float> facs(freqs.nelements(),1.0);
131 switch (inst) {
132 case ATMOPRA:
133 {
134 if (year<2003) {
135 pushLog("There is no beam efficiency data from before 2003"
136 " - using 2003 data");
137 facs = interp(freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2003Y_);
138 } else if (year==2003) {
139 pushLog("Using beam efficiency data from 2003");
140 facs = interp(freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2003Y_);
141 } else {
142 pushLog("Using beam efficiency data from 2004");
143 facs = interp(freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2004Y_);
144 }
145 }
146 break;
147 default:
148 {
149 pushLog("No beam efficiency data for this instrument - assuming unity");
150 }
151 }
152 return facs;
153}
154
155Vector<Float> STAttr::apertureEfficiency(Instrument inst,
156 const MEpoch& dateObs,
157 const Vector<Float>& freqs) const
158{
159
160 // Look at date where appropriate
161 MVTime t(dateObs.getValue());
162 uInt year = t.year();
163
164 Vector<Float> facs(freqs.nelements(),1.0);
165 switch (inst) {
166 case ATMOPRA:
167 {
168 if (year<2004) {
169 pushLog("There is no aperture efficiency data from before 2004"
170 " - using 2004 data");
171 facs = interp(freqs/1.0e9f, MopEtaApX_, MopEtaAp2004Y_);
172 } else {
173 pushLog("Using aperture efficiency data from 2004");
174 facs = interp(freqs/1.0e9f, MopEtaApX_, MopEtaAp2004Y_);
175 }
176 }
177 break;
178 case TIDBINBILLA:
179 {
180 facs = interp(freqs/1.0e9f, TidEtaApX_, TidEtaApY_);
181 }
182 break;
183 default:
184 {
185 pushLog("No aperture efficiency data for this instrument"
186 " - assuming unity");
187 }
188 }
189 return facs;
190}
191
192Vector<Float> STAttr::JyPerK(Instrument inst, const MEpoch& dateObs,
193 const Vector<Float>& freqs) const
194{
195
196 // Find what we need
197 Vector<Float> etaAp = apertureEfficiency(inst, dateObs, freqs);
198 Float D = diameter(inst);
199 // Compute it
200 Vector<Float> facs(freqs.nelements(),1.0);
201 for (uInt i=0; i<freqs.nelements(); i++) {
202 facs(i) = STAttr::findJyPerK(etaAp(i), D);
203 }
204 return facs;
205}
206
207
208Float STAttr::findJyPerK(Float etaAp, Float D)
209 //
210 // Converts K -> Jy
211 // D in m
212 //
213{
214 Double kb = QC::k.getValue(Unit(String("erg/K")));
215 Float gA = C::pi * D * D / 4.0;
216 return (2.0 * 1.0e19 * kb / etaAp / gA);
217}
218
219
220Vector<Float> STAttr::gainElevationPoly(Instrument inst) const
221{
222
223 // Look at date where appropriate
224 switch (inst) {
225 case TIDBINBILLA:
226 {
227 return TidGainElPoly_.copy();
228 }
229 break;
230 // assume HOH for K-band
231 case ATPKSHOH:
232 {
233 return ParkesGainElPoly_.copy();
234 }
235 break;
236 default:
237 {
238 Vector<Float> t;
239 return t.copy();
240 }
241 }
242}
243
244std::string STAttr::feedPolType(Instrument inst) const
245{
246 std::string type;
247 switch (inst) {
248 case ATMOPRA:
249 case ATPKSMB:
250 case ATPKSHOH:
251 {
252 type = "linear";
253 }
254 break;
255 case TIDBINBILLA:
256 {
257 type = "circular";
258 }
259 break;
260 default:
261 {
262 type = "linear";
263 }
264 }
265 return type;
266}
267
268
269// Private
270Vector<Float> STAttr::interp(const Vector<Float>& xOut,
271 const Vector<Float>& xIn,
272 const Vector<Float>& yIn) const
273{
274 Int method = 1; // Linear
275 Vector<Float> yOut;
276 Vector<Bool> mOut;
277
278 Vector<Bool> mIn(xIn.nelements(),True);
279
280 InterpolateArray1D<Float,Float>::interpolate(yOut, mOut, xOut,
281 xIn, yIn, mIn,
282 method, True, True);
283
284 return yOut;
285}
286
287void STAttr::initData()
288 //
289 // Mopra data from Mopra web page
290 // Tid data from Tid web page
291 // X in GHz
292 //
293{
294
295 // Beam efficiency
296 MopEtaBeamX_.resize(3);
297 MopEtaBeamX_(0) = 86.0;
298 MopEtaBeamX_(1) = 100.0;
299 MopEtaBeamX_(2) = 115.0;
300
301 MopEtaBeam2003Y_.resize(3);
302 MopEtaBeam2003Y_(0) = 0.39;
303 MopEtaBeam2003Y_(1) = 0.37;
304 MopEtaBeam2003Y_(2) = 0.37; // replicated from (1)
305
306 MopEtaBeam2004Y_.resize(3);
307 MopEtaBeam2004Y_(0) = 0.49;
308 MopEtaBeam2004Y_(1) = 0.44;
309 MopEtaBeam2004Y_(2) = 0.42;
310
311 // Aperture efficiency
312 MopEtaApX_.resize(2);
313 MopEtaApX_(0) = 86.0;
314 MopEtaApX_(1) = 115.0;
315
316 MopEtaAp2004Y_.resize(2);
317 MopEtaAp2004Y_(0) = 0.33;
318 MopEtaAp2004Y_(1) = 0.24;
319
320 TidEtaApX_.resize(2);
321 TidEtaApY_.resize(2);
322 TidEtaApX_(0) = 18.0;
323 TidEtaApX_(1) = 26.5;
324 TidEtaApY_(0) = 0.4848;
325 TidEtaApY_(1) = 0.4848;
326
327 // Gain elevation correction polynomial coefficients (for elevation
328 // in degrees)
329
330 TidGainElPoly_.resize(3);
331 TidGainElPoly_(0) = 3.58788e-1;
332 TidGainElPoly_(1) = 2.87243e-2;
333 TidGainElPoly_(2) = -3.219093e-4;
334
335 // 2009-09-15 - 13mm (22.2GHz) receiver
336 ParkesGainElPoly_.resize(3);
337 ParkesGainElPoly_(0) = -0.194031;
338 ParkesGainElPoly_(1) = 0.457724e-1;
339 ParkesGainElPoly_(2) = -0.438659e-3;
340}
341
342
343Instrument STAttr::convertInstrument(const String& instrument,
344 Bool throwIt)
345{
346 String t(instrument);
347 t.upcase();
348
349 // The strings are what STReader returns, after cunning
350 // interrogation of station names... :-(
351 Instrument inst = asap::UNKNOWNINST;
352 if (t==String("DSS-43")) {
353 inst = TIDBINBILLA;
354 } else if (t==String("ALMA")) {
355 inst = ALMA;
356 } else if (t==String("ATPKSMB")) {
357 inst = ATPKSMB;
358 } else if (t==String("ATPKSHOH")) {
359 inst = ATPKSHOH;
360 } else if (t==String("ATMOPRA")) {
361 inst = ATMOPRA;
362 } else if (t==String("CEDUNA")) {
363 inst = CEDUNA;
364 } else if (t==String("GBT")) {
365 inst = GBT;
366 } else if (t==String("HOBART")) {
367 inst = HOBART;
368 } else {
369 if (throwIt) {
370 throw AipsError("Unrecognized instrument"
371 " - use function scan.set_instrument to set");
372 }
373 }
374 return inst;
375}
Note: See TracBrowser for help on using the repository browser.