source: trunk/src/STAttr.cpp@ 2316

Last change on this file since 2316 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
RevLine 
[353]1//#---------------------------------------------------------------------------
[878]2//# STAttr.cc: A collection of attributes for different telescopes
[353]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>
[394]38#include <casa/Quanta/MVTime.h>
[353]39
40#include <measures/Measures/MEpoch.h>
41
42#include <scimath/Mathematics/InterpolateArray1D.h>
43
[878]44#include "STAttr.h"
[717]45
[353]46using namespace casa;
47using namespace asap;
48
[878]49STAttr::STAttr()
[353]50{
[361]51 initData();
[353]52}
53
[2163]54STAttr::STAttr(const STAttr& other):
55 Logger()
[353]56{
[2163]57 (void) other; //suppress unused warning
[717]58 initData(); // state just private 'static' data
[353]59}
60
[1188]61STAttr& STAttr::operator=(const STAttr& other)
[353]62{
63 if (this != &other) {
[717]64 ; // state just private 'static' data
[353]65 }
66 return *this;
67}
68
[878]69STAttr::~STAttr()
[353]70{;}
71
72
[878]73Float STAttr::diameter(Instrument inst) const
[353]74{
75 Float D = 1.0;
[362]76 switch (inst) {
[1391]77 case ALMA:
78 {
79 D = 12.0;
80 }
81 break;
82
[362]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;
[1391]104 case GBT:
105 {
106 D = 104.9;
107 }
108 break;
[362]109 case HOBART:
110 {
111 D = 26.0;
112 }
113 break;
114 default:
115 {
116 throw(AipsError("Unknown instrument"));
117 }
[353]118 }
119 return D;
120}
121
[1188]122Vector<Float> STAttr::beamEfficiency(Instrument inst, const MEpoch& dateObs,
[717]123 const Vector<Float>& freqs) const
[353]124{
125
[717]126 // Look at date where appropriate
[394]127 MVTime t(dateObs.getValue());
128 uInt year = t.year();
[717]129
[353]130 Vector<Float> facs(freqs.nelements(),1.0);
[362]131 switch (inst) {
[717]132 case ATMOPRA:
[362]133 {
[717]134 if (year<2003) {
135 pushLog("There is no beam efficiency data from before 2003"
136 " - using 2003 data");
[1188]137 facs = interp(freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2003Y_);
[717]138 } else if (year==2003) {
139 pushLog("Using beam efficiency data from 2003");
[1188]140 facs = interp(freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2003Y_);
[717]141 } else {
142 pushLog("Using beam efficiency data from 2004");
[1188]143 facs = interp(freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2004Y_);
[717]144 }
[362]145 }
146 break;
[717]147 default:
148 {
149 pushLog("No beam efficiency data for this instrument - assuming unity");
150 }
[353]151 }
152 return facs;
153}
154
[1188]155Vector<Float> STAttr::apertureEfficiency(Instrument inst,
156 const MEpoch& dateObs,
[717]157 const Vector<Float>& freqs) const
[353]158{
[1188]159
[717]160 // Look at date where appropriate
161 MVTime t(dateObs.getValue());
162 uInt year = t.year();
[1188]163
[717]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;
[353]190}
191
[878]192Vector<Float> STAttr::JyPerK(Instrument inst, const MEpoch& dateObs,
[717]193 const Vector<Float>& freqs) const
[353]194{
[1188]195
[717]196 // Find what we need
197 Vector<Float> etaAp = apertureEfficiency(inst, dateObs, freqs);
198 Float D = diameter(inst);
199 // Compute it
[353]200 Vector<Float> facs(freqs.nelements(),1.0);
201 for (uInt i=0; i<freqs.nelements(); i++) {
[878]202 facs(i) = STAttr::findJyPerK(etaAp(i), D);
[353]203 }
204 return facs;
205}
206
207
[878]208Float STAttr::findJyPerK(Float etaAp, Float D)
[717]209 //
210 // Converts K -> Jy
211 // D in m
212 //
[353]213{
[717]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);
[353]217}
218
219
[878]220Vector<Float> STAttr::gainElevationPoly(Instrument inst) const
[362]221{
[353]222
[717]223 // Look at date where appropriate
224 switch (inst) {
225 case TIDBINBILLA:
226 {
227 return TidGainElPoly_.copy();
228 }
229 break;
[1345]230 // assume HOH for K-band
231 case ATPKSHOH:
232 {
233 return ParkesGainElPoly_.copy();
234 }
235 break;
[717]236 default:
237 {
238 Vector<Float> t;
239 return t.copy();
240 }
241 }
[362]242}
243
[1188]244std::string STAttr::feedPolType(Instrument inst) const
[507]245{
[1188]246 std::string type;
[717]247 switch (inst) {
248 case ATMOPRA:
249 case ATPKSMB:
250 case ATPKSHOH:
251 {
[1188]252 type = "linear";
[717]253 }
254 break;
255 case TIDBINBILLA:
256 {
[1188]257 type = "circular";
[717]258 }
259 break;
260 default:
261 {
[1188]262 type = "linear";
[717]263 }
264 }
265 return type;
[507]266}
[362]267
268
[353]269// Private
[878]270Vector<Float> STAttr::interp(const Vector<Float>& xOut,
[1188]271 const Vector<Float>& xIn,
[717]272 const Vector<Float>& yIn) const
273{
274 Int method = 1; // Linear
275 Vector<Float> yOut;
276 Vector<Bool> mOut;
[353]277
[717]278 Vector<Bool> mIn(xIn.nelements(),True);
[1188]279
[717]280 InterpolateArray1D<Float,Float>::interpolate(yOut, mOut, xOut,
281 xIn, yIn, mIn,
282 method, True, True);
[1188]283
[717]284 return yOut;
[353]285}
286
[1188]287void STAttr::initData()
[717]288 //
289 // Mopra data from Mopra web page
290 // Tid data from Tid web page
291 // X in GHz
292 //
[362]293{
294
[717]295 // Beam efficiency
[353]296 MopEtaBeamX_.resize(3);
297 MopEtaBeamX_(0) = 86.0;
298 MopEtaBeamX_(1) = 100.0;
299 MopEtaBeamX_(2) = 115.0;
[717]300
[353]301 MopEtaBeam2003Y_.resize(3);
302 MopEtaBeam2003Y_(0) = 0.39;
303 MopEtaBeam2003Y_(1) = 0.37;
[717]304 MopEtaBeam2003Y_(2) = 0.37; // replicated from (1)
[1188]305
[353]306 MopEtaBeam2004Y_.resize(3);
307 MopEtaBeam2004Y_(0) = 0.49;
308 MopEtaBeam2004Y_(1) = 0.44;
[361]309 MopEtaBeam2004Y_(2) = 0.42;
[362]310
[717]311 // Aperture efficiency
[353]312 MopEtaApX_.resize(2);
313 MopEtaApX_(0) = 86.0;
314 MopEtaApX_(1) = 115.0;
[717]315
[353]316 MopEtaAp2004Y_.resize(2);
317 MopEtaAp2004Y_(0) = 0.33;
318 MopEtaAp2004Y_(1) = 0.24;
[717]319
[467]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;
[362]326
[717]327 // Gain elevation correction polynomial coefficients (for elevation
328 // in degrees)
[362]329
330 TidGainElPoly_.resize(3);
331 TidGainElPoly_(0) = 3.58788e-1;
332 TidGainElPoly_(1) = 2.87243e-2;
333 TidGainElPoly_(2) = -3.219093e-4;
[1345]334
[1618]335 // 2009-09-15 - 13mm (22.2GHz) receiver
[1345]336 ParkesGainElPoly_.resize(3);
[1618]337 ParkesGainElPoly_(0) = -0.194031;
338 ParkesGainElPoly_(1) = 0.457724e-1;
339 ParkesGainElPoly_(2) = -0.438659e-3;
[353]340}
[475]341
342
[878]343Instrument STAttr::convertInstrument(const String& instrument,
[475]344 Bool throwIt)
345{
[717]346 String t(instrument);
347 t.upcase();
[1188]348
[878]349 // The strings are what STReader returns, after cunning
[717]350 // interrogation of station names... :-(
351 Instrument inst = asap::UNKNOWNINST;
352 if (t==String("DSS-43")) {
353 inst = TIDBINBILLA;
[1391]354 } else if (t==String("ALMA")) {
355 inst = ALMA;
[717]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;
[1391]364 } else if (t==String("GBT")) {
365 inst = GBT;
[717]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;
[475]375}
Note: See TracBrowser for help on using the repository browser.