source: tags/asap-4.1.0/src/STAttr.cpp@ 2747

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

merge from trunk into release candidate

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