source: trunk/src/STAttr.cpp@ 3007

Last change on this file since 3007 was 2658, checked in by Malte Marquarding, 12 years ago

Ticket #199: Excised Logger::pushLog; now everything is using LogIO
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.8 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#include <casa/Logging/LogIO.h>
40
41#include <measures/Measures/MEpoch.h>
42
43#include <scimath/Mathematics/InterpolateArray1D.h>
44
45#include "STAttr.h"
46
47using namespace casa;
48using namespace asap;
49
50STAttr::STAttr()
51{
52 initData();
53}
54
55STAttr::STAttr(const STAttr& other)
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 casa::LogIO os( casa::LogOrigin( "STAttr", "beamEfficiency()" ) );
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 os << "There is no beam efficiency data from before 2003"
136 <<" - using 2003 data" << casa::LogIO::POST;
137 facs = interp(freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2003Y_);
138 } else if (year==2003) {
139 os << "Using beam efficiency data from 2003" << casa::LogIO::POST;
140 facs = interp(freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2003Y_);
141 } else {
142 os << "Using beam efficiency data from 2004" << casa::LogIO::POST;
143 facs = interp(freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2004Y_);
144 }
145 }
146 break;
147 default:
148 {
149 os << "No beam efficiency data for this instrument - assuming unity"
150 << casa::LogIO::POST;
151 }
152 }
153 return facs;
154}
155
156Vector<Float> STAttr::apertureEfficiency(Instrument inst,
157 const MEpoch& dateObs,
158 const Vector<Float>& freqs) const
159{
160 casa::LogIO os( casa::LogOrigin( "STAttr", "apertureEfficiency()" ) );
161 // Look at date where appropriate
162 MVTime t(dateObs.getValue());
163 uInt year = t.year();
164
165 Vector<Float> facs(freqs.nelements(),1.0);
166 switch (inst) {
167 case ATMOPRA:
168 {
169 if (year<2004) {
170 os << "There is no aperture efficiency data from before 2004"
171 << " - using 2004 data" << casa::LogIO::POST;
172 facs = interp(freqs/1.0e9f, MopEtaApX_, MopEtaAp2004Y_);
173 } else {
174 os << "Using aperture efficiency data from 2004" << casa::LogIO::POST;
175 facs = interp(freqs/1.0e9f, MopEtaApX_, MopEtaAp2004Y_);
176 }
177 }
178 break;
179 case TIDBINBILLA:
180 {
181 facs = interp(freqs/1.0e9f, TidEtaApX_, TidEtaApY_);
182 }
183 break;
184 default:
185 {
186 os << "No aperture efficiency data for this instrument"
187 << " - assuming unity" << casa::LogIO::POST;
188 }
189 }
190 return facs;
191}
192
193Vector<Float> STAttr::JyPerK(Instrument inst, const MEpoch& dateObs,
194 const Vector<Float>& freqs) const
195{
196
197 // Find what we need
198 Vector<Float> etaAp = apertureEfficiency(inst, dateObs, freqs);
199 Float D = diameter(inst);
200 // Compute it
201 Vector<Float> facs(freqs.nelements(),1.0);
202 for (uInt i=0; i<freqs.nelements(); i++) {
203 facs(i) = STAttr::findJyPerK(etaAp(i), D);
204 }
205 return facs;
206}
207
208
209Float STAttr::findJyPerK(Float etaAp, Float D)
210 //
211 // Converts K -> Jy
212 // D in m
213 //
214{
215 Double kb = QC::k.getValue(Unit(String("erg/K")));
216 Float gA = C::pi * D * D / 4.0;
217 return (2.0 * 1.0e19 * kb / etaAp / gA);
218}
219
220
221Vector<Float> STAttr::gainElevationPoly(Instrument inst) const
222{
223
224 // Look at date where appropriate
225 switch (inst) {
226 case TIDBINBILLA:
227 {
228 return TidGainElPoly_.copy();
229 }
230 break;
231 // assume HOH for K-band
232 case ATPKSHOH:
233 {
234 return ParkesGainElPoly_.copy();
235 }
236 break;
237 default:
238 {
239 Vector<Float> t;
240 return t.copy();
241 }
242 }
243}
244
245std::string STAttr::feedPolType(Instrument inst) const
246{
247 std::string type;
248 switch (inst) {
249 case ATMOPRA:
250 case ATPKSMB:
251 case ATPKSHOH:
252 {
253 type = "linear";
254 }
255 break;
256 case TIDBINBILLA:
257 {
258 type = "circular";
259 }
260 break;
261 default:
262 {
263 type = "linear";
264 }
265 }
266 return type;
267}
268
269
270// Private
271Vector<Float> STAttr::interp(const Vector<Float>& xOut,
272 const Vector<Float>& xIn,
273 const Vector<Float>& yIn) const
274{
275 Int method = 1; // Linear
276 Vector<Float> yOut;
277 Vector<Bool> mOut;
278
279 Vector<Bool> mIn(xIn.nelements(),True);
280
281 InterpolateArray1D<Float,Float>::interpolate(yOut, mOut, xOut,
282 xIn, yIn, mIn,
283 method, True, True);
284
285 return yOut;
286}
287
288void STAttr::initData()
289 //
290 // Mopra data from Mopra web page
291 // Tid data from Tid web page
292 // X in GHz
293 //
294{
295
296 // Beam efficiency
297 MopEtaBeamX_.resize(3);
298 MopEtaBeamX_(0) = 86.0;
299 MopEtaBeamX_(1) = 100.0;
300 MopEtaBeamX_(2) = 115.0;
301
302 MopEtaBeam2003Y_.resize(3);
303 MopEtaBeam2003Y_(0) = 0.39;
304 MopEtaBeam2003Y_(1) = 0.37;
305 MopEtaBeam2003Y_(2) = 0.37; // replicated from (1)
306
307 MopEtaBeam2004Y_.resize(3);
308 MopEtaBeam2004Y_(0) = 0.49;
309 MopEtaBeam2004Y_(1) = 0.44;
310 MopEtaBeam2004Y_(2) = 0.42;
311
312 // Aperture efficiency
313 MopEtaApX_.resize(2);
314 MopEtaApX_(0) = 86.0;
315 MopEtaApX_(1) = 115.0;
316
317 MopEtaAp2004Y_.resize(2);
318 MopEtaAp2004Y_(0) = 0.33;
319 MopEtaAp2004Y_(1) = 0.24;
320
321 TidEtaApX_.resize(2);
322 TidEtaApY_.resize(2);
323 TidEtaApX_(0) = 18.0;
324 TidEtaApX_(1) = 26.5;
325 TidEtaApY_(0) = 0.4848;
326 TidEtaApY_(1) = 0.4848;
327
328 // Gain elevation correction polynomial coefficients (for elevation
329 // in degrees)
330
331 TidGainElPoly_.resize(3);
332 TidGainElPoly_(0) = 3.58788e-1;
333 TidGainElPoly_(1) = 2.87243e-2;
334 TidGainElPoly_(2) = -3.219093e-4;
335
336 // 2009-09-15 - 13mm (22.2GHz) receiver
337 ParkesGainElPoly_.resize(3);
338 ParkesGainElPoly_(0) = -0.194031;
339 ParkesGainElPoly_(1) = 0.457724e-1;
340 ParkesGainElPoly_(2) = -0.438659e-3;
341}
342
343
344Instrument STAttr::convertInstrument(const String& instrument,
345 Bool throwIt)
346{
347 String t(instrument);
348 t.upcase();
349
350 // The strings are what STReader returns, after cunning
351 // interrogation of station names... :-(
352 Instrument inst = asap::UNKNOWNINST;
353 if (t==String("DSS-43")) {
354 inst = TIDBINBILLA;
355 } else if (t==String("ALMA")) {
356 inst = ALMA;
357 } else if (t==String("ATPKSMB")) {
358 inst = ATPKSMB;
359 } else if (t==String("ATPKSHOH")) {
360 inst = ATPKSHOH;
361 } else if (t==String("ATMOPRA")) {
362 inst = ATMOPRA;
363 } else if (t==String("CEDUNA")) {
364 inst = CEDUNA;
365 } else if (t==String("GBT")) {
366 inst = GBT;
367 } else if (t==String("HOBART")) {
368 inst = HOBART;
369 } else {
370 if (throwIt) {
371 throw AipsError("Unrecognized instrument"
372 " - use function scan.set_instrument to set");
373 }
374 }
375 return inst;
376}
Note: See TracBrowser for help on using the repository browser.