Line data    Source code

       1             : //# CTIter.cc: Implementation of CTIter.h
       2             : //# Copyright (C) 2011 
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: aips2-request@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //----------------------------------------------------------------------------
      27             : 
      28             : #include <casacore/measures/Measures/Muvw.h>
      29             : #include <casacore/measures/Measures/MCBaseline.h>
      30             : #include <msvis/MSVis/ViImplementation2.h>
      31             : #include <synthesis/CalTables/CTIter.h>
      32             : #include <casacore/tables/Tables/ScalarColumn.h>
      33             : #include <casacore/casa/Arrays.h>
      34             : #include <casacore/casa/OS/Timer.h>
      35             : 
      36             : using namespace casacore;
      37             : namespace casa { //# NAMESPACE CASA - BEGIN
      38             : 
      39             : //----------------------------------------------------------------------------
      40             : 
      41         880 : ROCTIter::ROCTIter(NewCalTable tab, const Block<String>& sortcol) :
      42             :   sortCols_(sortcol.begin( ),sortcol.end( )),
      43             :   singleSpw_(false),
      44             :   parentNCT_(tab),
      45             :   calCol_(NULL),
      46             :   ti_(NULL),
      47             :   inct_(NULL),
      48             :   iROCTMainCols_(NULL),
      49             :   init_epoch_(false),
      50         880 :   init_uvw_(false)
      51             : {
      52         880 :   calCol_=new ROCTColumns(tab);
      53         880 :   ti_=new TableIterator(tab,sortcol);
      54             :   // Attach initial accessors:
      55         880 :   attach();
      56             : 
      57             :   // If SPW a sort column, then 
      58         880 :   singleSpw_=anyEQ(sortCols_,String("SPECTRAL_WINDOW_ID"));
      59             : 
      60             :   /*
      61             :   cout << "singleSpw_ = " << boolalpha << singleSpw_ << endl;
      62             :   cout << "inct_->nrow() = " << inct_->nrow() << endl;
      63             :   cout << "this->nrow() = " << this->nrow() << endl;
      64             :   cout << "iROCTMainCols_->spwId() = " << iROCTMainCols_->spwId().getColumn() << endl;
      65             :   cout << "iROCTMainCols_->spwId()(0) = " << iROCTMainCols_->spwId()(0) << endl;
      66             :   cout << "thisSpw() = " << this->thisSpw() << endl;
      67             : 
      68             :   cout << "done." << endl;
      69             :   */
      70         880 : };
      71             : 
      72             : //----------------------------------------------------------------------------
      73             : 
      74         880 : ROCTIter::~ROCTIter()
      75             : {
      76         880 :   if (calCol_!=NULL) delete calCol_;
      77         880 :   if (ti_!=NULL) delete ti_;
      78         880 :   if (iROCTMainCols_!=NULL) delete iROCTMainCols_;
      79         880 :   if (inct_!=NULL) delete inct_;
      80         880 : };
      81             : 
      82       54155 : void ROCTIter::next() { 
      83             :   // Advance the TableIterator
      84       54155 :   ti_->next();
      85             : 
      86             :   // attach accessors to new iteration
      87       54155 :   this->attach();
      88       54155 : };
      89             : 
      90           0 : void ROCTIter::next0() { 
      91             :   // Advance the TableIterator
      92           0 :   ti_->next();
      93           0 : };
      94             : 
      95           0 : void ROCTIter::setCTColumns(const NewCalTable& tab) {
      96             :   // Set subtable columns from another table
      97           0 :   delete calCol_;
      98           0 :   calCol_ = new ROCTColumns(tab);
      99           0 : }
     100             : 
     101        8897 : Double ROCTIter::thisTime() const { return iROCTMainCols_->time()(0); };
     102       10770 : Vector<Double> ROCTIter::time() const { return iROCTMainCols_->time().getColumn(); };
     103           0 : void ROCTIter::time(Vector<Double>& v) const { iROCTMainCols_->time().getColumn(v); };
     104             : 
     105       16806 : Int ROCTIter::thisField() const { return iROCTMainCols_->fieldId()(0); };
     106       38055 : Vector<Int> ROCTIter::field() const { return iROCTMainCols_->fieldId().getColumn(); };
     107           0 : void ROCTIter::field(Vector<Int>& v) const { iROCTMainCols_->fieldId().getColumn(v); };
     108             : 
     109       49151 : Int ROCTIter::thisSpw() const { return iROCTMainCols_->spwId()(0); };
     110       10383 : Vector<Int> ROCTIter::spw() const { return iROCTMainCols_->spwId().getColumn(); };
     111           0 : void ROCTIter::spw(Vector<Int>& v) const { iROCTMainCols_->spwId().getColumn(v); };
     112             : 
     113          60 : Int ROCTIter::thisScan() const { return iROCTMainCols_->scanNo()(0); };
     114       10383 : Vector<Int> ROCTIter::scan() const { return iROCTMainCols_->scanNo().getColumn(); };
     115           0 : void ROCTIter::scan(Vector<Int>& v) const { iROCTMainCols_->scanNo().getColumn(v); };
     116             : 
     117           0 : Int ROCTIter::thisObs() const { return iROCTMainCols_->obsId()(0); };
     118       10383 : Vector<Int> ROCTIter::obs() const { return iROCTMainCols_->obsId().getColumn(); };
     119           0 : void ROCTIter::obs(Vector<Int>& v) const { iROCTMainCols_->obsId().getColumn(v); };
     120             : 
     121       20125 : Int ROCTIter::thisAntenna1() const { return iROCTMainCols_->antenna1()(0); };
     122       16302 : Vector<Int> ROCTIter::antenna1() const { return iROCTMainCols_->antenna1().getColumn(); };
     123           0 : void ROCTIter::antenna1(Vector<Int>& v) const { iROCTMainCols_->antenna1().getColumn(v); };
     124         240 : Int ROCTIter::thisAntenna2() const { return iROCTMainCols_->antenna2()(0); };
     125       12584 : Vector<Int> ROCTIter::antenna2() const { return iROCTMainCols_->antenna2().getColumn(); };
     126           0 : void ROCTIter::antenna2(Vector<Int>& v) const { iROCTMainCols_->antenna2().getColumn(v); };
     127             : 
     128           0 : Double ROCTIter::thisInterval() const { return iROCTMainCols_->interval()(0); };
     129           0 : Vector<Double> ROCTIter::interval() const { return iROCTMainCols_->interval().getColumn(); };
     130           0 : void ROCTIter::interval(Vector<Double>& v) const { iROCTMainCols_->interval().getColumn(v); };
     131             : 
     132       23875 : Cube<Complex> ROCTIter::cparam() const { return iROCTMainCols_->cparam().getColumn(); };
     133        5286 : void ROCTIter::cparam(Cube<Complex>& c) const { iROCTMainCols_->cparam().getColumn(c); };
     134         919 : Cube<Float> ROCTIter::fparam() const { return iROCTMainCols_->fparam().getColumn(); };
     135           0 : void ROCTIter::fparam(Cube<Float>& f) const { iROCTMainCols_->fparam().getColumn(f); };
     136             : 
     137             : // Complex as Float
     138         348 : Cube<Float> ROCTIter::casfparam(String what) const { return iROCTMainCols_->fparamArray(what); };
     139           0 : void ROCTIter::casfparam(Cube<Float>& f,String what) const { iROCTMainCols_->fparamArray(f,what); };
     140             : 
     141           0 : Cube<Float> ROCTIter::paramErr() const { return iROCTMainCols_->paramerr().getColumn(); };
     142        5265 : void ROCTIter::paramErr(Cube<Float>& c) const { iROCTMainCols_->paramerr().getColumn(c); };
     143             : 
     144         236 : Cube<Float> ROCTIter::snr() const { return iROCTMainCols_->snr().getColumn(); };
     145        5265 : void ROCTIter::snr(Cube<Float>& c) const { iROCTMainCols_->snr().getColumn(c); };
     146           0 : Cube<Float> ROCTIter::wt() const { return iROCTMainCols_->weight().getColumn(); };
     147           0 : void ROCTIter::wt(Cube<Float>& c) const { iROCTMainCols_->weight().getColumn(c); };
     148             : 
     149       38735 : Cube<Bool> ROCTIter::flag() const { return iROCTMainCols_->flag().getColumn(); };
     150           0 : void ROCTIter::flag(Cube<Bool>& c) const { iROCTMainCols_->flag().getColumn(c); };
     151             : 
     152           0 : Vector<Int> ROCTIter::chan() const {
     153           0 :   Vector<Int> chans;
     154           0 :   this->chan(chans);
     155           0 :   return chans;
     156             : }
     157             : 
     158        5265 : Int ROCTIter::nchan() const {
     159        5265 :   if (singleSpw_)
     160        5265 :     return calCol_->spectralWindow().numChan()(this->thisSpw());
     161             :   else
     162             :     // more than one spw per iteration...
     163           0 :     throw(AipsError("Please sort by spw."));
     164             : }
     165             : 
     166           0 : void ROCTIter::chan(Vector<Int>& v) const {
     167           0 :   if (singleSpw_) {
     168           0 :     v.resize(calCol_->spectralWindow().numChan()(this->thisSpw()));
     169             :     // TBD: observe channel selection here:
     170           0 :     indgen(v);
     171             :   }
     172             :   else
     173             :     // more than one spw per iteration...
     174           0 :     throw(AipsError("Please sort by spw."));
     175           0 : }
     176             : 
     177         116 : Vector<Double> ROCTIter::freq() const {
     178         116 :   Vector<Double> freqs;
     179         116 :   this->freq(freqs);
     180         116 :   return freqs;
     181             : }
     182             : 
     183         194 : void ROCTIter::freq(Vector<Double>& v) const {
     184         194 :   if (singleSpw_) {
     185         194 :     v.resize();
     186         194 :     calCol_->spectralWindow().chanFreq().get(this->thisSpw(),v);
     187             :   }
     188             :   else
     189             :     // more than one spw per iteration...
     190           0 :     throw(AipsError("Please sort by spw."));
     191         194 : }
     192             : 
     193           0 : int ROCTIter::freqFrame(int spwId) const {
     194           0 :   int frame = calCol_->spectralWindow().measFreqRef()(spwId);
     195           0 :   return frame;
     196             : }
     197             : 
     198       55656 : void ROCTIter::attach() {
     199             :   // Attach accessors for current iteration:
     200       55656 :   if (iROCTMainCols_!=NULL) delete iROCTMainCols_;
     201       55656 :   if (inct_!=NULL) delete inct_;
     202             : 
     203       55656 :   inct_= new NewCalTable(ti_->table());
     204       55656 :   iROCTMainCols_ = new ROCTMainColumns(*inct_);
     205       55656 : }
     206             : 
     207           0 : casacore::MDirection ROCTIter::azel0(casacore::Double time) {
     208           0 :   if (!init_epoch_) {
     209           0 :     initEpoch();
     210             :   }
     211             : 
     212             :   try {
     213           0 :     updatePhaseCenter();
     214           0 :   } catch (const casacore::AipsError& err) {
     215           0 :     throw(casacore::AipsError("azel: " + err.getMesg()));
     216             :   }
     217             : 
     218           0 :   casacore::MSDerivedValues msd;
     219           0 :   msd.setAntennas(calCol_->antenna());
     220           0 :   msd.setFieldCenter(phaseCenter_);
     221             : 
     222           0 :   casacore::MDirection azel;
     223           0 :   vi::ViImplementation2::azel0Calculate(time, msd, azel, epoch_);
     224           0 :   return azel;
     225             : }
     226             : 
     227           0 : casacore::Double ROCTIter::hourang(casacore::Double time) {
     228           0 :   if (!init_epoch_) {
     229           0 :     initEpoch();
     230             :   }
     231             : 
     232             :   try {
     233           0 :     updatePhaseCenter();
     234           0 :   } catch (const casacore::AipsError& err) {
     235           0 :     throw(casacore::AipsError("hourang: " + err.getMesg()));
     236             :   }
     237             : 
     238           0 :   casacore::MSDerivedValues msd;
     239           0 :   msd.setAntennas(calCol_->antenna());
     240           0 :   msd.setFieldCenter(phaseCenter_);
     241             : 
     242           0 :   return vi::ViImplementation2::hourangCalculate(time, msd, epoch_);
     243             : }
     244             : 
     245           0 : casacore::Float ROCTIter::parang0(casacore::Double time) {
     246           0 :   if (!init_epoch_) {
     247           0 :     initEpoch();
     248             :   }
     249             : 
     250             :   try {
     251           0 :     updatePhaseCenter();
     252           0 :   } catch (const casacore::AipsError& err) {
     253           0 :     throw(casacore::AipsError("parang: " + err.getMesg()));
     254             :   }
     255             : 
     256           0 :   casacore::MSDerivedValues msd;
     257           0 :   msd.setAntennas(calCol_->antenna());
     258           0 :   msd.setFieldCenter(phaseCenter_);
     259             : 
     260           0 :   return vi::ViImplementation2::parang0Calculate(time, msd, epoch_);
     261             : }
     262             : 
     263           0 : casacore::Matrix<casacore::Double> ROCTIter::uvw() {
     264           0 :   casacore::Vector<casacore::Int> ant2 = antenna2();
     265           0 :   if (ant2(0) == -1) {
     266           0 :     throw(AipsError("UVW axes cannot be calculated for antenna-based calibration tables."));  
     267             :   }
     268             : 
     269           0 :   updateAntennaUVW();
     270             : 
     271           0 :   casacore::Vector<casacore::Int> ant1 = antenna1();
     272           0 :   auto nbaseline = ant1.size();
     273           0 :   casacore::Matrix<casacore::Double> uvw(3, nbaseline);
     274             : 
     275           0 :   for (uInt i = 0; i < nbaseline; ++i) {
     276           0 :     uvw.column(i) = antennaUVW_[ant2(i)] - antennaUVW_[ant1(i)];
     277             :   }
     278             : 
     279           0 :   return uvw;
     280             : }
     281             : 
     282           0 : void ROCTIter::updatePhaseCenter() {
     283             :   // Set MSDerivedValues phase center for field
     284           0 :   if ((thisTime() == 0) || (thisField() < 0)) {
     285           0 :     throw(AipsError("cannot calculate this value with no valid timestamp or field ID."));
     286             :   }
     287             : 
     288           0 :   phaseCenter_ = calCol_->field().phaseDirMeas(thisField(), thisTime());
     289           0 : }
     290             : 
     291           0 : void ROCTIter::initEpoch() {
     292           0 :   epoch_ = iROCTMainCols_->timeMeas()(0);
     293           0 :   init_epoch_ = true;
     294           0 : }
     295             : 
     296           0 : void ROCTIter::initUVW() {
     297             :   // Calculate relative positions of antennas
     298           0 :   nAnt_ = calCol_->antenna().nrow();
     299           0 :   auto antPosMeas = calCol_->antenna().positionMeas();
     300           0 :   refAntPos_ = antPosMeas(0); // use first antenna for reference
     301           0 :   auto refAntPosValue = refAntPos_.getValue();
     302             : 
     303             :   // Set up baseline and uvw types
     304           0 :   MBaseline::getType(baseline_type_, MPosition::showType(refAntPos_.getRef().getType()));
     305             : 
     306           0 :   mvbaselines_.resize(nAnt_);
     307             : 
     308           0 :   for (Int ant = 0; ant < nAnt_; ++ant) {
     309             :     // MVBaselines are basically xyz Vectors, not Measures
     310           0 :     mvbaselines_[ant] = MVBaseline(refAntPosValue, antPosMeas(ant).getValue());    
     311             :   }
     312             : 
     313           0 :   init_uvw_ = true;
     314           0 : }
     315             : 
     316           0 : void ROCTIter::updateAntennaUVW() {
     317             :   // Set antennaUVW_ for current iteration
     318           0 :   if (!init_uvw_) {
     319           0 :     initUVW();
     320             :   }
     321             : 
     322           0 :   updatePhaseCenter();
     323             : 
     324           0 :   MeasFrame measFrame(refAntPos_, epoch_, phaseCenter_);
     325             : 
     326             :   // Antenna frame
     327           0 :   MBaseline::Ref baselineRef(baseline_type_);
     328           0 :   MVBaseline mvbaseline;
     329           0 :   MBaseline baselineMeas;
     330           0 :   baselineMeas.set(mvbaseline, baselineRef);
     331           0 :   baselineMeas.getRefPtr()->set(measFrame);
     332             : 
     333             :   // Conversion engine to phasedir type
     334             :   casacore::MBaseline::Types phasedir_type;
     335           0 :   MBaseline::getType(phasedir_type, MDirection::showType(phaseCenter_.getRef().getType()));
     336           0 :   MBaseline::Ref uvwRef(phasedir_type);
     337           0 :   MBaseline::Convert baselineConv(baselineMeas, uvwRef);
     338             : 
     339             :   // WSRT convention: phase opposite to VLA (l increases toward increasing RA)
     340           0 :   Int obsId = thisObs();
     341           0 :   if (obsId < 0) obsId = 0;
     342           0 :   casacore::String telescope = calCol_->observation().telescopeName()(obsId);
     343           0 :   bool wsrtConvention = (telescope == "WSRT");
     344             : 
     345           0 :   antennaUVW_.resize(nAnt_);
     346             : 
     347           0 :   for (int i = 0; i < nAnt_; ++i) {
     348           0 :     baselineMeas.set(mvbaselines_[i], baselineRef);
     349           0 :     MBaseline baselineOutFrame = baselineConv(baselineMeas);
     350           0 :     MVuvw uvwOutFrame(baselineOutFrame.getValue(), phaseCenter_.getValue());
     351             : 
     352           0 :     if (wsrtConvention) {
     353           0 :       antennaUVW_[i] = -uvwOutFrame.getValue();
     354             :     } else {
     355           0 :       antennaUVW_[i] = uvwOutFrame.getValue();
     356             :     }
     357             :   }
     358           0 : }
     359             : 
     360             : 
     361             : // CTIter
     362             : 
     363         419 : CTIter::CTIter(NewCalTable tab, const Block<String>& sortcol) :
     364             :   ROCTIter(tab,sortcol),
     365             :   irwnct_(NULL),
     366         419 :   iRWCTMainCols_(NULL)
     367             : {
     368             :   // Attach first iteration
     369             :   //  TBD: this unnecessarily redoes the ROCTIter attach...
     370         419 :   attach();
     371         419 : }
     372             : 
     373         621 : CTIter::~CTIter() {
     374         419 :   if (iRWCTMainCols_!=NULL) delete iRWCTMainCols_;
     375         419 :   if (irwnct_!=NULL) delete irwnct_;
     376         621 : }
     377             : 
     378             : // Set fieldid
     379        7768 : void CTIter::setfield(Int fieldid) {
     380        7768 :   iRWCTMainCols_->fieldId().fillColumn(fieldid); 
     381        7768 : }
     382             : 
     383             : // Set scan number
     384        7768 : void CTIter::setscan(Int scan) {
     385        7768 :   iRWCTMainCols_->scanNo().fillColumn(scan); 
     386        7768 : }
     387             : 
     388             : // Set obsid
     389        7768 : void CTIter::setobs(Int obs) {
     390        7768 :   iRWCTMainCols_->obsId().fillColumn(obs); 
     391        7768 : }
     392             : 
     393             : // Set antenna2 (e.g., used for setting refant)
     394        8034 : void CTIter::setantenna2(const Vector<Int>& a2) {
     395        8034 :   iRWCTMainCols_->antenna2().putColumn(a2); 
     396        8034 : }
     397             : 
     398        6879 : void CTIter::setflag(const Cube<Bool>& fl) {
     399        6879 :   iRWCTMainCols_->flag().putColumn(fl); 
     400        6879 : }
     401             : 
     402         227 : void CTIter::setfparam(const Cube<Float>& f) {
     403         227 :   iRWCTMainCols_->fparam().putColumn(f); 
     404         227 : };
     405             : 
     406        8956 : void CTIter::setcparam(const Cube<Complex>& c) {
     407        8956 :   iRWCTMainCols_->cparam().putColumn(c); 
     408        8956 : };
     409             : 
     410       29380 : void CTIter::attach() {
     411             : 
     412             :   // Attach readonly access
     413       29380 :   ROCTIter::attach();
     414             : 
     415             :   // Attach writable access
     416       29380 :   if (iRWCTMainCols_!=NULL) delete iRWCTMainCols_;
     417       29380 :   if (irwnct_!=NULL) delete irwnct_;
     418       29380 :   irwnct_= new NewCalTable(this->table());
     419       29380 :   iRWCTMainCols_ = new CTMainColumns(*irwnct_);
     420       29380 : }
     421             : 
     422             : 
     423             : 
     424             : } //# NAMESPACE CASA - END