Changeset 717 for trunk/src/SDAttr.cc


Ignore:
Timestamp:
11/17/05 14:37:54 (19 years ago)
Author:
mar637
Message:

implemented use of SDLog

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/SDAttr.cc

    r507 r717  
    3030//#---------------------------------------------------------------------------
    3131
    32 #include "SDAttr.h"
    3332#include <casa/aips.h>
    3433#include <casa/Arrays/Vector.h>
    … …  
    4342#include <scimath/Mathematics/InterpolateArray1D.h>
    4443
     44#include "SDAttr.h"
     45
    4546using namespace casa;
    4647using namespace asap;
    4748
    48 SDAttr::SDAttr ()
     49SDAttr::SDAttr()
    4950{
    5051   initData();
    … …  
    5354SDAttr::SDAttr(const SDAttr& other)
    5455{
    55    initData();                                 // state just private 'static' data
     56   initData();                     // state just private 'static' data
    5657}
    5758
    … …  
    5960{
    6061  if (this != &other) {
    61     ;                                      // state just private 'static' data
     62    ;                             // state just private 'static' data
    6263  }
    6364  return *this;
    … …  
    6869
    6970
    70 Float SDAttr::diameter (Instrument inst)  const
     71Float SDAttr::diameter(Instrument inst)  const
    7172{
    7273   Float D = 1.0;
    … …  
    103104        }
    104105   }
    105 //
    106106   return D;
    107107}
    108108
    109 Vector<Float> SDAttr::beamEfficiency (Instrument inst, const MEpoch& dateObs, const Vector<Float>& freqs) const
    110 {
    111 
    112 // Look at date where appropriate
    113 
     109Vector<Float> SDAttr::beamEfficiency(Instrument inst, const MEpoch& dateObs,
     110                                      const Vector<Float>& freqs) const
     111{
     112
     113  // Look at date where appropriate
    114114   MVTime t(dateObs.getValue());
    115115   uInt year = t.year();
    116 //
     116
    117117   Vector<Float> facs(freqs.nelements(),1.0);
    118118   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 
    142 Vector<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         }
     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     }
    172138   }
    173139   return facs;
    174140}
    175141
    176 Vector<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 
     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
    186187   Vector<Float> facs(freqs.nelements(),1.0);
    187188   for (uInt i=0; i<freqs.nelements(); i++) {
    188       facs(i) = SDAttr::findJyPerK (etaAp(i), D);
     189     facs(i) = SDAttr::findJyPerK(etaAp(i), D);
    189190   }
    190 //
    191191   return facs;
    192192}
    193193
    194194
    195 Float 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 
    207 Vector<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 FeedPolType 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 
     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}
    252248
    253249
    254250// Private
    255 
    256 Vector<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 
    272 void 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 
     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
    282277   MopEtaBeamX_.resize(3);
    283278   MopEtaBeamX_(0) = 86.0;
    284279   MopEtaBeamX_(1) = 100.0;
    285280   MopEtaBeamX_(2) = 115.0;
    286 //
     281
    287282   MopEtaBeam2003Y_.resize(3);
    288283   MopEtaBeam2003Y_(0) = 0.39;
    289284   MopEtaBeam2003Y_(1) = 0.37;
    290    MopEtaBeam2003Y_(2) = 0.37;                // replicated from (1)
    291 //
     285   MopEtaBeam2003Y_(2) = 0.37;          // replicated from (1)
     286   
    292287   MopEtaBeam2004Y_.resize(3);
    293288   MopEtaBeam2004Y_(0) = 0.49;
    … …  
    295290   MopEtaBeam2004Y_(2) = 0.42;
    296291
    297 // Aperture efficiency
    298 
     292   // Aperture efficiency
    299293   MopEtaApX_.resize(2);
    300294   MopEtaApX_(0) = 86.0;
    301295   MopEtaApX_(1) = 115.0;
    302 //
     296
    303297   MopEtaAp2004Y_.resize(2);
    304298   MopEtaAp2004Y_(0) = 0.33;
    305299   MopEtaAp2004Y_(1) = 0.24;
    306 //
     300
    307301   TidEtaApX_.resize(2);
    308302   TidEtaApY_.resize(2);
    … …  
    312306   TidEtaApY_(1) = 0.4848;
    313307
    314 // Gain elevation correction polynomial coefficients (for elevation in degrees)
     308   // Gain elevation correction polynomial coefficients (for elevation
     309   // in degrees)
    315310
    316311   TidGainElPoly_.resize(3);
    … …  
    324319                                     Bool throwIt)
    325320{
    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 
     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 TracChangeset for help on using the changeset viewer.