source: trunk/src/STAttr.cpp@ 931

Last change on this file since 931 was 878, checked in by mar637, 19 years ago

SDAttr->STAttr name change

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.9 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
[878]54STAttr::STAttr(const STAttr& other)
[353]55{
[717]56 initData(); // state just private 'static' data
[353]57}
58
[878]59STAttr& STAttr::operator=(const STAttr& other)
[353]60{
61 if (this != &other) {
[717]62 ; // state just private 'static' data
[353]63 }
64 return *this;
65}
66
[878]67STAttr::~STAttr()
[353]68{;}
69
70
[878]71Float STAttr::diameter(Instrument inst) const
[353]72{
73 Float D = 1.0;
[362]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 }
[353]105 }
106 return D;
107}
108
[878]109Vector<Float> STAttr::beamEfficiency(Instrument inst, const MEpoch& dateObs,
[717]110 const Vector<Float>& freqs) const
[353]111{
112
[717]113 // Look at date where appropriate
[394]114 MVTime t(dateObs.getValue());
115 uInt year = t.year();
[717]116
[353]117 Vector<Float> facs(freqs.nelements(),1.0);
[362]118 switch (inst) {
[717]119 case ATMOPRA:
[362]120 {
[717]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 }
[362]132 }
133 break;
[717]134 default:
135 {
136 pushLog("No beam efficiency data for this instrument - assuming unity");
137 }
[353]138 }
139 return facs;
140}
141
[878]142Vector<Float> STAttr::apertureEfficiency(Instrument inst,
[717]143 const MEpoch& dateObs,
144 const Vector<Float>& freqs) const
[353]145{
[717]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;
[353]177}
178
[878]179Vector<Float> STAttr::JyPerK(Instrument inst, const MEpoch& dateObs,
[717]180 const Vector<Float>& freqs) const
[353]181{
[717]182
183 // Find what we need
184 Vector<Float> etaAp = apertureEfficiency(inst, dateObs, freqs);
185 Float D = diameter(inst);
186 // Compute it
[353]187 Vector<Float> facs(freqs.nelements(),1.0);
188 for (uInt i=0; i<freqs.nelements(); i++) {
[878]189 facs(i) = STAttr::findJyPerK(etaAp(i), D);
[353]190 }
191 return facs;
192}
193
194
[878]195Float STAttr::findJyPerK(Float etaAp, Float D)
[717]196 //
197 // Converts K -> Jy
198 // D in m
199 //
[353]200{
[717]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);
[353]204}
205
206
[878]207Vector<Float> STAttr::gainElevationPoly(Instrument inst) const
[362]208{
[353]209
[717]210 // Look at date where appropriate
211 switch (inst) {
212 case TIDBINBILLA:
213 {
214 return TidGainElPoly_.copy();
215 }
216 break;
217 default:
218 {
219 Vector<Float> t;
220 return t.copy();
221 }
222 }
[362]223}
224
[878]225FeedPolType STAttr::feedPolType(Instrument inst) const
[507]226{
[717]227 FeedPolType type = UNKNOWNFEED;
228 switch (inst) {
229 case ATMOPRA:
230 case ATPKSMB:
231 case ATPKSHOH:
232 {
233 type = LINEAR;
234 }
235 break;
236 case TIDBINBILLA:
237 {
238 type = CIRCULAR;
239 }
240 break;
241 default:
242 {
243 type = UNKNOWNFEED;
244 }
245 }
246 return type;
[507]247}
[362]248
249
[353]250// Private
[878]251Vector<Float> STAttr::interp(const Vector<Float>& xOut,
[717]252 const Vector<Float>& xIn,
253 const Vector<Float>& yIn) const
254{
255 Int method = 1; // Linear
256 Vector<Float> yOut;
257 Vector<Bool> mOut;
[353]258
[717]259 Vector<Bool> mIn(xIn.nelements(),True);
260
261 InterpolateArray1D<Float,Float>::interpolate(yOut, mOut, xOut,
262 xIn, yIn, mIn,
263 method, True, True);
264
265 return yOut;
[353]266}
267
[878]268void STAttr::initData()
[717]269 //
270 // Mopra data from Mopra web page
271 // Tid data from Tid web page
272 // X in GHz
273 //
[362]274{
275
[717]276 // Beam efficiency
[353]277 MopEtaBeamX_.resize(3);
278 MopEtaBeamX_(0) = 86.0;
279 MopEtaBeamX_(1) = 100.0;
280 MopEtaBeamX_(2) = 115.0;
[717]281
[353]282 MopEtaBeam2003Y_.resize(3);
283 MopEtaBeam2003Y_(0) = 0.39;
284 MopEtaBeam2003Y_(1) = 0.37;
[717]285 MopEtaBeam2003Y_(2) = 0.37; // replicated from (1)
286
[353]287 MopEtaBeam2004Y_.resize(3);
288 MopEtaBeam2004Y_(0) = 0.49;
289 MopEtaBeam2004Y_(1) = 0.44;
[361]290 MopEtaBeam2004Y_(2) = 0.42;
[362]291
[717]292 // Aperture efficiency
[353]293 MopEtaApX_.resize(2);
294 MopEtaApX_(0) = 86.0;
295 MopEtaApX_(1) = 115.0;
[717]296
[353]297 MopEtaAp2004Y_.resize(2);
298 MopEtaAp2004Y_(0) = 0.33;
299 MopEtaAp2004Y_(1) = 0.24;
[717]300
[467]301 TidEtaApX_.resize(2);
302 TidEtaApY_.resize(2);
303 TidEtaApX_(0) = 18.0;
304 TidEtaApX_(1) = 26.5;
305 TidEtaApY_(0) = 0.4848;
306 TidEtaApY_(1) = 0.4848;
[362]307
[717]308 // Gain elevation correction polynomial coefficients (for elevation
309 // in degrees)
[362]310
311 TidGainElPoly_.resize(3);
312 TidGainElPoly_(0) = 3.58788e-1;
313 TidGainElPoly_(1) = 2.87243e-2;
314 TidGainElPoly_(2) = -3.219093e-4;
[353]315}
[475]316
317
[878]318Instrument STAttr::convertInstrument(const String& instrument,
[475]319 Bool throwIt)
320{
[717]321 String t(instrument);
322 t.upcase();
323
[878]324 // The strings are what STReader returns, after cunning
[717]325 // interrogation of station names... :-(
326 Instrument inst = asap::UNKNOWNINST;
327 if (t==String("DSS-43")) {
328 inst = TIDBINBILLA;
329 } else if (t==String("ATPKSMB")) {
330 inst = ATPKSMB;
331 } else if (t==String("ATPKSHOH")) {
332 inst = ATPKSHOH;
333 } else if (t==String("ATMOPRA")) {
334 inst = ATMOPRA;
335 } else if (t==String("CEDUNA")) {
336 inst = CEDUNA;
337 } else if (t==String("HOBART")) {
338 inst = HOBART;
339 } else {
340 if (throwIt) {
341 throw AipsError("Unrecognized instrument"
342 " - use function scan.set_instrument to set");
343 }
344 }
345 return inst;
[475]346}
Note: See TracBrowser for help on using the repository browser.