source: trunk/src/SDAttr.cc @ 475

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