source: trunk/src/SDAttr.cc@ 490

Last change on this file since 490 was 475, checked in by kil064, 20 years ago

move function convertInstrument into here from SDMemTable
fix error in apertureEfficiency function (missing break)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.0 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 "SDAttr.h"
33#include <casa/aips.h>
34#include <casa/Arrays/Vector.h>
35#include <casa/Arrays/ArrayMath.h>
36#include <casa/Exceptions.h>
37#include <casa/Quanta/QC.h>
38#include <casa/Quanta/Quantum.h>
39#include <casa/Quanta/MVTime.h>
40
41#include <measures/Measures/MEpoch.h>
42
43#include <scimath/Mathematics/InterpolateArray1D.h>
44
45using namespace casa;
46using namespace asap;
47
48SDAttr::SDAttr ()
49{
50 initData();
51}
52
53SDAttr::SDAttr(const SDAttr& other)
54{
55 initData(); // state just private 'static' data
56}
57
58SDAttr& SDAttr::operator=(const SDAttr& other)
59{
60 if (this != &other) {
61 ; // state just private 'static' data
62 }
63 return *this;
64}
65
66SDAttr::~SDAttr()
67{;}
68
69
70Float SDAttr::diameter (Instrument inst) const
71{
72 Float D = 1.0;
73 switch (inst) {
74 case ATMOPRA:
75 {
76 D = 22.0;
77 }
78 break;
79 case ATPKSMB:
80 case ATPKSHOH:
81 {
82 D = 64.0;
83 }
84 break;
85 case TIDBINBILLA:
86 {
87 D = 70.0;
88 }
89 break;
90 case CEDUNA:
91 {
92 D = 30.0;
93 }
94 break;
95 case HOBART:
96 {
97 D = 26.0;
98 }
99 break;
100 default:
101 {
102 throw(AipsError("Unknown instrument"));
103 }
104 }
105//
106 return D;
107}
108
109Vector<Float> SDAttr::beamEfficiency (Instrument inst, const MEpoch& dateObs, const Vector<Float>& freqs) const
110{
111
112// Look at date where appropriate
113
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 cerr << "There is no beam efficiency data from before 2003 - using 2003 data" << endl;
123 facs = interp (freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2003Y_);
124 } else if (year==2003) {
125 cerr << "Using beam efficiency data from 2003" << endl;
126 facs = interp (freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2003Y_);
127 } else {
128 cerr << "Using beam efficiency data from 2004" << endl;
129 facs = interp (freqs/1.0e9f, MopEtaBeamX_, MopEtaBeam2004Y_);
130 }
131 }
132 break;
133 default:
134 {
135 cerr << "No beam efficiency data for this instrument - assuming unity" << endl;
136 }
137 }
138//
139 return facs;
140}
141
142Vector<Float> SDAttr::apertureEfficiency (Instrument inst, const MEpoch& dateObs, const Vector<Float>& freqs) const
143{
144
145// Look at date where appropriate
146
147 MVTime t(dateObs.getValue());
148 uInt year = t.year();
149//
150 Vector<Float> facs(freqs.nelements(),1.0);
151 switch (inst) {
152 case ATMOPRA:
153 {
154 if (year<2004) {
155 cerr << "There is no aperture efficiency data from before 2004 - using 2004 data" << endl;
156 facs = interp (freqs/1.0e9f, MopEtaApX_, MopEtaAp2004Y_);
157 } else {
158 cerr << "Using aperture efficiency data from 2004" << endl;
159 facs = interp (freqs/1.0e9f, MopEtaApX_, MopEtaAp2004Y_);
160 }
161 }
162 break;
163 case TIDBINBILLA:
164 {
165 facs = interp (freqs/1.0e9f, TidEtaApX_, TidEtaApY_);
166 }
167 break;
168 default:
169 {
170 cerr << "No aperture efficiency data for this instrument - assuming unity" << endl;
171 }
172 }
173 return facs;
174}
175
176Vector<Float> SDAttr::JyPerK (Instrument inst, const MEpoch& dateObs, const Vector<Float>& freqs) const
177{
178
179// FInd what we need
180
181 Vector<Float> etaAp = apertureEfficiency (inst, dateObs, freqs);
182 Float D = diameter(inst);
183
184// Compute it
185
186 Vector<Float> facs(freqs.nelements(),1.0);
187 for (uInt i=0; i<freqs.nelements(); i++) {
188 facs(i) = SDAttr::findJyPerK (etaAp(i), D);
189 }
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
212 switch (inst) {
213 case TIDBINBILLA:
214 {
215 return TidGainElPoly_.copy();
216 }
217 break;
218 default:
219 {
220 Vector<Float> t;
221 return t.copy();
222 }
223 }
224}
225
226
227
228
229// Private
230
231Vector<Float> SDAttr::interp (const Vector<Float>& xOut, const Vector<Float>& xIn,
232 const Vector<Float>& yIn) const
233{
234 Int method = 1; // Linear
235 Vector<Float> yOut;
236 Vector<Bool> mOut;
237//
238 Vector<Bool> mIn(xIn.nelements(),True);
239//
240 InterpolateArray1D<Float,Float>::interpolate(yOut, mOut, xOut,
241 xIn, yIn, mIn,
242 method, True, True);
243//
244 return yOut;
245}
246
247void SDAttr::initData ()
248//
249// Mopra data from Mopra web page
250// Tid data from Tid web page
251// X in GHz
252//
253{
254
255// Beam efficiency
256
257 MopEtaBeamX_.resize(3);
258 MopEtaBeamX_(0) = 86.0;
259 MopEtaBeamX_(1) = 100.0;
260 MopEtaBeamX_(2) = 115.0;
261//
262 MopEtaBeam2003Y_.resize(3);
263 MopEtaBeam2003Y_(0) = 0.39;
264 MopEtaBeam2003Y_(1) = 0.37;
265 MopEtaBeam2003Y_(2) = 0.37; // replicated from (1)
266//
267 MopEtaBeam2004Y_.resize(3);
268 MopEtaBeam2004Y_(0) = 0.49;
269 MopEtaBeam2004Y_(1) = 0.44;
270 MopEtaBeam2004Y_(2) = 0.42;
271
272// Aperture efficiency
273
274 MopEtaApX_.resize(2);
275 MopEtaApX_(0) = 86.0;
276 MopEtaApX_(1) = 115.0;
277//
278 MopEtaAp2004Y_.resize(2);
279 MopEtaAp2004Y_(0) = 0.33;
280 MopEtaAp2004Y_(1) = 0.24;
281//
282 TidEtaApX_.resize(2);
283 TidEtaApY_.resize(2);
284 TidEtaApX_(0) = 18.0;
285 TidEtaApX_(1) = 26.5;
286 TidEtaApY_(0) = 0.4848;
287 TidEtaApY_(1) = 0.4848;
288
289// Gain elevation correction polynomial coefficients (for elevation in degrees)
290
291 TidGainElPoly_.resize(3);
292 TidGainElPoly_(0) = 3.58788e-1;
293 TidGainElPoly_(1) = 2.87243e-2;
294 TidGainElPoly_(2) = -3.219093e-4;
295}
296
297
298Instrument SDAttr::convertInstrument(const String& instrument,
299 Bool throwIt)
300{
301 String t(instrument);
302 t.upcase();
303
304 // The strings are what SDReader returns, after cunning interrogation
305 // of station names... :-(
306 Instrument inst = asap::UNKNOWN;
307 if (t==String("DSS-43")) {
308 inst = TIDBINBILLA;
309 } else if (t==String("ATPKSMB")) {
310 inst = ATPKSMB;
311 } else if (t==String("ATPKSHOH")) {
312 inst = ATPKSHOH;
313 } else if (t==String("ATMOPRA")) {
314 inst = ATMOPRA;
315 } else if (t==String("CEDUNA")) {
316 inst = CEDUNA;
317 } else if (t==String("HOBART")) {
318 inst = HOBART;
319 } else {
320 if (throwIt) {
321 throw AipsError("Unrecognized instrument - use function scan.set_instrument to set");
322 }
323 }
324 return inst;
325}
326
327
Note: See TracBrowser for help on using the repository browser.