source: branches/Release-2-fixes/src/SDAttr.cc@ 3037

Last change on this file since 3037 was 507, checked in by kil064, 20 years ago

add function feedPolType
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.4 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
226FeedPolType SDAttr::feedPolType (Instrument inst) const
227{
228 FeedPolType type = UNKNOWNFEED;
229 switch (inst) {
230 case ATMOPRA:
231 case ATPKSMB:
232 case ATPKSHOH:
233 {
234 type = LINEAR;
235 }
236 break;
237 case TIDBINBILLA:
238 {
239 type = CIRCULAR;
240 }
241 break;
242 default:
243 {
244 type = UNKNOWNFEED;
245 }
246 }
247 return type;
248}
249
250
251
252
253
254// Private
255
256Vector<Float> SDAttr::interp (const Vector<Float>& xOut, const Vector<Float>& xIn,
257 const Vector<Float>& yIn) const
258{
259 Int method = 1; // Linear
260 Vector<Float> yOut;
261 Vector<Bool> mOut;
262//
263 Vector<Bool> mIn(xIn.nelements(),True);
264//
265 InterpolateArray1D<Float,Float>::interpolate(yOut, mOut, xOut,
266 xIn, yIn, mIn,
267 method, True, True);
268//
269 return yOut;
270}
271
272void SDAttr::initData ()
273//
274// Mopra data from Mopra web page
275// Tid data from Tid web page
276// X in GHz
277//
278{
279
280// Beam efficiency
281
282 MopEtaBeamX_.resize(3);
283 MopEtaBeamX_(0) = 86.0;
284 MopEtaBeamX_(1) = 100.0;
285 MopEtaBeamX_(2) = 115.0;
286//
287 MopEtaBeam2003Y_.resize(3);
288 MopEtaBeam2003Y_(0) = 0.39;
289 MopEtaBeam2003Y_(1) = 0.37;
290 MopEtaBeam2003Y_(2) = 0.37; // replicated from (1)
291//
292 MopEtaBeam2004Y_.resize(3);
293 MopEtaBeam2004Y_(0) = 0.49;
294 MopEtaBeam2004Y_(1) = 0.44;
295 MopEtaBeam2004Y_(2) = 0.42;
296
297// Aperture efficiency
298
299 MopEtaApX_.resize(2);
300 MopEtaApX_(0) = 86.0;
301 MopEtaApX_(1) = 115.0;
302//
303 MopEtaAp2004Y_.resize(2);
304 MopEtaAp2004Y_(0) = 0.33;
305 MopEtaAp2004Y_(1) = 0.24;
306//
307 TidEtaApX_.resize(2);
308 TidEtaApY_.resize(2);
309 TidEtaApX_(0) = 18.0;
310 TidEtaApX_(1) = 26.5;
311 TidEtaApY_(0) = 0.4848;
312 TidEtaApY_(1) = 0.4848;
313
314// Gain elevation correction polynomial coefficients (for elevation in degrees)
315
316 TidGainElPoly_.resize(3);
317 TidGainElPoly_(0) = 3.58788e-1;
318 TidGainElPoly_(1) = 2.87243e-2;
319 TidGainElPoly_(2) = -3.219093e-4;
320}
321
322
323Instrument SDAttr::convertInstrument(const String& instrument,
324 Bool throwIt)
325{
326 String t(instrument);
327 t.upcase();
328
329 // The strings are what SDReader returns, after cunning interrogation
330 // of station names... :-(
331 Instrument inst = asap::UNKNOWNINST;
332 if (t==String("DSS-43")) {
333 inst = TIDBINBILLA;
334 } else if (t==String("ATPKSMB")) {
335 inst = ATPKSMB;
336 } else if (t==String("ATPKSHOH")) {
337 inst = ATPKSHOH;
338 } else if (t==String("ATMOPRA")) {
339 inst = ATMOPRA;
340 } else if (t==String("CEDUNA")) {
341 inst = CEDUNA;
342 } else if (t==String("HOBART")) {
343 inst = HOBART;
344 } else {
345 if (throwIt) {
346 throw AipsError("Unrecognized instrument - use function scan.set_instrument to set");
347 }
348 }
349 return inst;
350}
351
352
Note: See TracBrowser for help on using the repository browser.