source: trunk/src/SDAttr.cc @ 507

Last change on this file since 507 was 507, checked in by kil064, 19 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.