source: trunk/src/STAttr.cpp@ 1544

Last change on this file since 1544 was 1391, checked in by Malte Marquarding, 17 years ago

merge from alma branch to get alma/GBT support. Commented out fluxUnit changes as they are using a chnaged interface to PKSreader/writer. Also commented out progress meter related code.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.4 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{
56 initData(); // state just private 'static' data
57}
58
59STAttr& STAttr::operator=(const STAttr& other)
60{
61 if (this != &other) {
62 ; // state just private 'static' data
63 }
64 return *this;
65}
66
67STAttr::~STAttr()
68{;}
69
70
71Float STAttr::diameter(Instrument inst) const
72{
73 Float D = 1.0;
74 switch (inst) {
75 case ALMA:
76 {
77 D = 12.0;
78 }
79 break;
80
81 case ATMOPRA:
82 {
83 D = 22.0;
84 }
85 break;
86 case ATPKSMB:
87 case ATPKSHOH:
88 {
89 D = 64.0;
90 }
91 break;
92 case TIDBINBILLA:
93 {
94 D = 70.0;
95 }
96 break;
97 case CEDUNA:
98 {
99 D = 30.0;
100 }
101 break;
102 case GBT:
103 {
104 D = 104.9;
105 }
106 break;
107 case HOBART:
108 {
109 D = 26.0;
110 }
111 break;
112 default:
113 {
114 throw(AipsError("Unknown instrument"));
115 }
116 }
117 return D;
118}
119
120Vector<Float> STAttr::beamEfficiency(Instrument inst, const MEpoch& dateObs,
121 const Vector<Float>& freqs) const
122{
123
124 // Look at date where appropriate
125 MVTime t(dateObs.getValue());
126 uInt year = t.year();
127
128 Vector<Float> facs(freqs.nelements(),1.0);
129 switch (inst) {
130 case ATMOPRA:
131 {
132 if (year<2003) {
133 pushLog("There is no beam efficiency data from before 2003"
134 " - using 2003 data");
135 facs = interp(freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2003Y_);
136 } else if (year==2003) {
137 pushLog("Using beam efficiency data from 2003");
138 facs = interp(freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2003Y_);
139 } else {
140 pushLog("Using beam efficiency data from 2004");
141 facs = interp(freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2004Y_);
142 }
143 }
144 break;
145 default:
146 {
147 pushLog("No beam efficiency data for this instrument - assuming unity");
148 }
149 }
150 return facs;
151}
152
153Vector<Float> STAttr::apertureEfficiency(Instrument inst,
154 const MEpoch& dateObs,
155 const Vector<Float>& freqs) const
156{
157
158 // Look at date where appropriate
159 MVTime t(dateObs.getValue());
160 uInt year = t.year();
161
162 Vector<Float> facs(freqs.nelements(),1.0);
163 switch (inst) {
164 case ATMOPRA:
165 {
166 if (year<2004) {
167 pushLog("There is no aperture efficiency data from before 2004"
168 " - using 2004 data");
169 facs = interp(freqs/1.0e9f, MopEtaApX_, MopEtaAp2004Y_);
170 } else {
171 pushLog("Using aperture efficiency data from 2004");
172 facs = interp(freqs/1.0e9f, MopEtaApX_, MopEtaAp2004Y_);
173 }
174 }
175 break;
176 case TIDBINBILLA:
177 {
178 facs = interp(freqs/1.0e9f, TidEtaApX_, TidEtaApY_);
179 }
180 break;
181 default:
182 {
183 pushLog("No aperture efficiency data for this instrument"
184 " - assuming unity");
185 }
186 }
187 return facs;
188}
189
190Vector<Float> STAttr::JyPerK(Instrument inst, const MEpoch& dateObs,
191 const Vector<Float>& freqs) const
192{
193
194 // Find what we need
195 Vector<Float> etaAp = apertureEfficiency(inst, dateObs, freqs);
196 Float D = diameter(inst);
197 // Compute it
198 Vector<Float> facs(freqs.nelements(),1.0);
199 for (uInt i=0; i<freqs.nelements(); i++) {
200 facs(i) = STAttr::findJyPerK(etaAp(i), D);
201 }
202 return facs;
203}
204
205
206Float STAttr::findJyPerK(Float etaAp, Float D)
207 //
208 // Converts K -> Jy
209 // D in m
210 //
211{
212 Double kb = QC::k.getValue(Unit(String("erg/K")));
213 Float gA = C::pi * D * D / 4.0;
214 return (2.0 * 1.0e19 * kb / etaAp / gA);
215}
216
217
218Vector<Float> STAttr::gainElevationPoly(Instrument inst) const
219{
220
221 // Look at date where appropriate
222 switch (inst) {
223 case TIDBINBILLA:
224 {
225 return TidGainElPoly_.copy();
226 }
227 break;
228 // assume HOH for K-band
229 case ATPKSHOH:
230 {
231 return ParkesGainElPoly_.copy();
232 }
233 break;
234 default:
235 {
236 Vector<Float> t;
237 return t.copy();
238 }
239 }
240}
241
242std::string STAttr::feedPolType(Instrument inst) const
243{
244 std::string type;
245 switch (inst) {
246 case ATMOPRA:
247 case ATPKSMB:
248 case ATPKSHOH:
249 {
250 type = "linear";
251 }
252 break;
253 case TIDBINBILLA:
254 {
255 type = "circular";
256 }
257 break;
258 default:
259 {
260 type = "linear";
261 }
262 }
263 return type;
264}
265
266
267// Private
268Vector<Float> STAttr::interp(const Vector<Float>& xOut,
269 const Vector<Float>& xIn,
270 const Vector<Float>& yIn) const
271{
272 Int method = 1; // Linear
273 Vector<Float> yOut;
274 Vector<Bool> mOut;
275
276 Vector<Bool> mIn(xIn.nelements(),True);
277
278 InterpolateArray1D<Float,Float>::interpolate(yOut, mOut, xOut,
279 xIn, yIn, mIn,
280 method, True, True);
281
282 return yOut;
283}
284
285void STAttr::initData()
286 //
287 // Mopra data from Mopra web page
288 // Tid data from Tid web page
289 // X in GHz
290 //
291{
292
293 // Beam efficiency
294 MopEtaBeamX_.resize(3);
295 MopEtaBeamX_(0) = 86.0;
296 MopEtaBeamX_(1) = 100.0;
297 MopEtaBeamX_(2) = 115.0;
298
299 MopEtaBeam2003Y_.resize(3);
300 MopEtaBeam2003Y_(0) = 0.39;
301 MopEtaBeam2003Y_(1) = 0.37;
302 MopEtaBeam2003Y_(2) = 0.37; // replicated from (1)
303
304 MopEtaBeam2004Y_.resize(3);
305 MopEtaBeam2004Y_(0) = 0.49;
306 MopEtaBeam2004Y_(1) = 0.44;
307 MopEtaBeam2004Y_(2) = 0.42;
308
309 // Aperture efficiency
310 MopEtaApX_.resize(2);
311 MopEtaApX_(0) = 86.0;
312 MopEtaApX_(1) = 115.0;
313
314 MopEtaAp2004Y_.resize(2);
315 MopEtaAp2004Y_(0) = 0.33;
316 MopEtaAp2004Y_(1) = 0.24;
317
318 TidEtaApX_.resize(2);
319 TidEtaApY_.resize(2);
320 TidEtaApX_(0) = 18.0;
321 TidEtaApX_(1) = 26.5;
322 TidEtaApY_(0) = 0.4848;
323 TidEtaApY_(1) = 0.4848;
324
325 // Gain elevation correction polynomial coefficients (for elevation
326 // in degrees)
327
328 TidGainElPoly_.resize(3);
329 TidGainElPoly_(0) = 3.58788e-1;
330 TidGainElPoly_(1) = 2.87243e-2;
331 TidGainElPoly_(2) = -3.219093e-4;
332
333 ParkesGainElPoly_.resize(3);
334 ParkesGainElPoly_(0) = 0.296759e-1;
335 ParkesGainElPoly_(1) = -0.293124e-3;
336 ParkesGainElPoly_(2) = 0.264295e-6;
337}
338
339
340Instrument STAttr::convertInstrument(const String& instrument,
341 Bool throwIt)
342{
343 String t(instrument);
344 t.upcase();
345
346 // The strings are what STReader returns, after cunning
347 // interrogation of station names... :-(
348 Instrument inst = asap::UNKNOWNINST;
349 if (t==String("DSS-43")) {
350 inst = TIDBINBILLA;
351 } else if (t==String("ALMA")) {
352 inst = ALMA;
353 } else if (t==String("ATPKSMB")) {
354 inst = ATPKSMB;
355 } else if (t==String("ATPKSHOH")) {
356 inst = ATPKSHOH;
357 } else if (t==String("ATMOPRA")) {
358 inst = ATMOPRA;
359 } else if (t==String("CEDUNA")) {
360 inst = CEDUNA;
361 } else if (t==String("GBT")) {
362 inst = GBT;
363 } else if (t==String("HOBART")) {
364 inst = HOBART;
365 } else {
366 if (throwIt) {
367 throw AipsError("Unrecognized instrument"
368 " - use function scan.set_instrument to set");
369 }
370 }
371 return inst;
372}
Note: See TracBrowser for help on using the repository browser.