source: trunk/src/STAttr.cpp@ 1378

Last change on this file since 1378 was 1345, checked in by mar637, 18 years ago

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