source: branches/Release12/src/SDAttr.cc@ 3037

Last change on this file since 3037 was 717, checked in by mar637, 19 years ago

implemented use of SDLog
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.9 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDAttr.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 "SDAttr.h"
45
46using namespace casa;
47using namespace asap;
48
49SDAttr::SDAttr()
50{
51 initData();
52}
53
54SDAttr::SDAttr(const SDAttr& other)
55{
56 initData(); // state just private 'static' data
57}
58
59SDAttr& SDAttr::operator=(const SDAttr& other)
60{
61 if (this != &other) {
62 ; // state just private 'static' data
63 }
64 return *this;
65}
66
67SDAttr::~SDAttr()
68{;}
69
70
71Float SDAttr::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> SDAttr::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> SDAttr::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> SDAttr::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) = SDAttr::findJyPerK(etaAp(i), D);
190 }
191 return facs;
192}
193
194
195Float SDAttr::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> SDAttr::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 default:
218 {
219 Vector<Float> t;
220 return t.copy();
221 }
222 }
223}
224
225FeedPolType SDAttr::feedPolType(Instrument inst) const
226{
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;
247}
248
249
250// Private
251Vector<Float> SDAttr::interp(const Vector<Float>& xOut,
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;
258
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;
266}
267
268void SDAttr::initData()
269 //
270 // Mopra data from Mopra web page
271 // Tid data from Tid web page
272 // X in GHz
273 //
274{
275
276 // Beam efficiency
277 MopEtaBeamX_.resize(3);
278 MopEtaBeamX_(0) = 86.0;
279 MopEtaBeamX_(1) = 100.0;
280 MopEtaBeamX_(2) = 115.0;
281
282 MopEtaBeam2003Y_.resize(3);
283 MopEtaBeam2003Y_(0) = 0.39;
284 MopEtaBeam2003Y_(1) = 0.37;
285 MopEtaBeam2003Y_(2) = 0.37; // replicated from (1)
286
287 MopEtaBeam2004Y_.resize(3);
288 MopEtaBeam2004Y_(0) = 0.49;
289 MopEtaBeam2004Y_(1) = 0.44;
290 MopEtaBeam2004Y_(2) = 0.42;
291
292 // Aperture efficiency
293 MopEtaApX_.resize(2);
294 MopEtaApX_(0) = 86.0;
295 MopEtaApX_(1) = 115.0;
296
297 MopEtaAp2004Y_.resize(2);
298 MopEtaAp2004Y_(0) = 0.33;
299 MopEtaAp2004Y_(1) = 0.24;
300
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;
307
308 // Gain elevation correction polynomial coefficients (for elevation
309 // in degrees)
310
311 TidGainElPoly_.resize(3);
312 TidGainElPoly_(0) = 3.58788e-1;
313 TidGainElPoly_(1) = 2.87243e-2;
314 TidGainElPoly_(2) = -3.219093e-4;
315}
316
317
318Instrument SDAttr::convertInstrument(const String& instrument,
319 Bool throwIt)
320{
321 String t(instrument);
322 t.upcase();
323
324 // The strings are what SDReader returns, after cunning
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;
346}
Note: See TracBrowser for help on using the repository browser.