From 5b8eb9c6c41a95e93193c6173fdbf28a9cd6aec8 Mon Sep 17 00:00:00 2001
From: Ger van Diepen <diepen@astron.nl>
Date: Wed, 29 Sep 2010 13:55:25 +0000
Subject: [PATCH] bug 1436: Created a W-projection gridder that applies the
corrections
---
CEP/Imager/LofarGridder/CMakeLists.txt | 7 +-
.../include/LofarGridder/LofarGridder.h | 95 ++++++-
.../LofarGridder/LofarWProjectGridder.h | 102 ++++++++
CEP/Imager/LofarGridder/src/CMakeLists.txt | 2 +-
CEP/Imager/LofarGridder/src/LofarGridder.cc | 135 ++++++++--
.../LofarGridder/src/LofarWProjectGridder.cc | 231 ++++++++++++++++++
CEP/Imager/LofarGridder/src/Register.cc | 5 +-
CEP/Imager/LofarGridder/test/tLofarGridder.cc | 2 +-
8 files changed, 545 insertions(+), 34 deletions(-)
create mode 100644 CEP/Imager/LofarGridder/include/LofarGridder/LofarWProjectGridder.h
create mode 100644 CEP/Imager/LofarGridder/src/LofarWProjectGridder.cc
diff --git a/CEP/Imager/LofarGridder/CMakeLists.txt b/CEP/Imager/LofarGridder/CMakeLists.txt
index de66fe92d19..09119a11735 100644
--- a/CEP/Imager/LofarGridder/CMakeLists.txt
+++ b/CEP/Imager/LofarGridder/CMakeLists.txt
@@ -1,9 +1,12 @@
# $Id$
-lofar_package(LofarGridder 0.1 DEPENDS MWImager)
+#lofar_package(LofarGridder 0.1 DEPENDS MWImager)
+#lofar_package(LofarGridder 0.1 DEPENDS BBSKernel)
lofar_find_package(AskapSoft REQUIRED)
-lofar_find_package(Casacore COMPONENTS images REQUIRED)
+lofar_find_package(LOFAR COMPONENTS bbskernel parmdb mwcommon blob common REQUIRED)
+lofar_find_package(Casacore COMPONENTS ms images REQUIRED)
+lofar_find_package(Boost REQUIRED)
add_subdirectory(include/LofarGridder)
add_subdirectory(src)
diff --git a/CEP/Imager/LofarGridder/include/LofarGridder/LofarGridder.h b/CEP/Imager/LofarGridder/include/LofarGridder/LofarGridder.h
index 6991ec0dd2b..460b94b7331 100644
--- a/CEP/Imager/LofarGridder/include/LofarGridder/LofarGridder.h
+++ b/CEP/Imager/LofarGridder/include/LofarGridder/LofarGridder.h
@@ -1,4 +1,4 @@
-//# LofarGridder.h: Test visibility gridder.
+//# LofarGridder.h: Gridder for LOFAR data correcting for DD effects
//#
//# Copyright (C) 2009
//# ASTRON (Netherlands Institute for Radio Astronomy)
@@ -22,23 +22,26 @@
//#
//# @author Ger van Diepen <diepen at astron dot nl>
-#ifndef LOFAR_LOFARGRIDDER_TESTGRIDDER_H
-#define LOFAR_LOFARGRIDDER_TESTGRIDDER_H
+#ifndef LOFAR_LOFARGRIDDER_LOFARGRIDDER_H
+#define LOFAR_LOFARGRIDDER_LOFARGRIDDER_H
#include <gridding/TableVisGridder.h>
#include <Common/ParameterSet.h>
+#include <Common/lofar_vector.h>
+#include <Common/lofar_string.h>
namespace LOFAR
{
- // @brief Gridder to test dynamic loading of gridders
+ // @brief Gridder for LOFAR data correcting for DD effects
//
// @ingroup testgridder
- class LofarGridder : public askap::synthesis::TableVisGridder
+
+ class LofarGridder : public askap::synthesis::IVisGridder
{
public:
- // Standard two dimensional box gridding
- LofarGridder();
+ // Construct from the given parset.
+ explicit LofarGridder (const ParameterSet&);
// Clone this Gridder
virtual askap::synthesis::IVisGridder::ShPtr clone();
@@ -55,19 +58,89 @@ namespace LOFAR
// @brief Register the gridder create function with its name.
static void registerGridder();
+ // @brief Initialise the gridding
+ // @param axes axes specifications
+ // @param shape Shape of output image: u,v,pol,chan
+ // @param dopsf Make the psf?
+ virtual void initialiseGrid (const askap::scimath::Axes& axes,
+ const casa::IPosition& shape,
+ const bool dopsf=true);
+
+ // @brief Grid the visibility data.
+ // @param acc const data accessor to work with
+ // @note a non-const adapter is created behind the scene. If no on-the-fly
+ // visibility correction is performed, this adapter is equivalent to the
+ // original const data accessor.
+ virtual void grid (askap::synthesis::IConstDataAccessor& acc);
+
+ // Form the final output image
+ // @param out Output double precision image or PSF
+ virtual void finaliseGrid (casa::Array<double>& out);
+
+ // @brief Calculate weights image
+ // @details Form the sum of the convolution function squared,
+ // multiplied by the weights for each different convolution
+ // function. This is used in the evaluation of the position
+ // dependent sensitivity
+ // @param out Output double precision sum of weights images
+ virtual void finaliseWeights (casa::Array<double>& out);
+
+ // @brief Initialise the degridding
+ // @param axes axes specifications
+ // @param image Input image: cube: u,v,pol,chan
+ virtual void initialiseDegrid (const askap::scimath::Axes& axes,
+ const casa::Array<double>& image);
+
+ // @brief Make context-dependant changes to the gridder behaviour
+ // @param context context
+ virtual void customiseForContext (casa::String context);
+
+ // @brief assign weights
+ // @param viswt shared pointer to visibility weights
+ virtual void initVisWeights (askap::synthesis::IVisWeights::ShPtr viswt);
+
+ // @brief Degrid the visibility data.
+ // @param[in] acc non-const data accessor to work with
+ virtual void degrid (askap::synthesis::IDataAccessor& acc);
+
+ // @brief Finalise degridding.
+ virtual void finaliseDegrid();
+
// @brief Initialise the indices
// @param[in] acc const data accessor to work with
- virtual void initIndices(const askap::synthesis::IConstDataAccessor& acc);
+ /// virtual void initIndices (const askap::synthesis::IConstDataAccessor& acc);
// @brief Correct for gridding convolution function
// @param image image to be corrected
- virtual void correctConvolution(casa::Array<double>& image);
+ /// virtual void correctConvolution(casa::Array<double>& image);
protected:
// Initialize convolution function
// @param[in] acc const data accessor to work with
- virtual void initConvolutionFunction
- (const askap::synthesis::IConstDataAccessor& acc);
+ /// virtual void initConvolutionFunction
+ /// (const askap::synthesis::IConstDataAccessor& acc);
+
+ private:
+ // Initialize the corrections.
+ void initCorrections (const askap::synthesis::IConstDataAccessor& acc);
+
+
+ //# Data members.
+ // Is the gridder correction initialized?
+ bool itsInitialized;
+ // Nr of frequency channels to average.
+ uint itsFreqAvg;
+ // Nr of time slots to average.
+ uint itsTimeAvg;
+ // The effects to correct for.
+ vector<string> itsCorrect;
+ // The ParmDB to use.
+ string itsParmDBName;
+ // The gridder to use after correcting and averaging.
+ //# Unfortunatelty we need to use TableVisGridder because the
+ //# function getImageCntre is available only there.
+ askap::synthesis::IVisGridder::ShPtr itsIGridder;
+ askap::synthesis::TableVisGridder* itsGridder;
};
} //# end namespace
diff --git a/CEP/Imager/LofarGridder/include/LofarGridder/LofarWProjectGridder.h b/CEP/Imager/LofarGridder/include/LofarGridder/LofarWProjectGridder.h
new file mode 100644
index 00000000000..dfcf78b5e10
--- /dev/null
+++ b/CEP/Imager/LofarGridder/include/LofarGridder/LofarWProjectGridder.h
@@ -0,0 +1,102 @@
+//# LofarWProjectGridder.h: Gridder for LOFAR data correcting for DD effects
+//#
+//# Copyright (C) 2009
+//# ASTRON (Netherlands Institute for Radio Astronomy)
+//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//#
+//# This file is part of the LOFAR software suite.
+//# The LOFAR software suite is free software: you can redistribute it and/or
+//# modify it under the terms of the GNU General Public License as published
+//# by the Free Software Foundation, either version 3 of the License, or
+//# (at your option) any later version.
+//#
+//# The LOFAR software suite is distributed in the hope that it will be useful,
+//# but WITHOUT ANY WARRANTY; without even the implied warranty of
+//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+//# GNU General Public License for more details.
+//#
+//# You should have received a copy of the GNU General Public License along
+//# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//#
+//# $Id$
+//#
+//# @author Ger van Diepen <diepen at astron dot nl>
+
+#ifndef LOFAR_LOFARGRIDDER_LOFARWPROJECTGRIDDER_H
+#define LOFAR_LOFARGRIDDER_LOFARWPROJECTGRIDDER_H
+
+//# LOFAR includes
+#include <BBSKernel/StationResponse.h>
+#include <Common/ParameterSet.h>
+#include <Common/lofar_vector.h>
+#include <Common/lofar_string.h>
+
+//# askap includes
+#include <gridding/WProjectVisGridder.h>
+
+//# casacore includes
+#include <ms/MeasurementSets/MeasurementSet.h>
+#include <casa/Quanta/MVDirection.h>
+
+namespace LOFAR
+{
+ // @brief Gridder for LOFAR data correcting for DD effects
+ //
+ // @ingroup testgridder
+
+ class LofarWProjectGridder : public askap::synthesis::WProjectVisGridder
+ {
+ public:
+
+ // Construct from the given parset.
+ explicit LofarWProjectGridder (const ParameterSet&);
+
+ // Clone this Gridder
+ virtual askap::synthesis::IVisGridder::ShPtr clone();
+
+ virtual ~LofarWProjectGridder();
+
+ // @brief Function to create the gridder from a parset.
+ // This function will be registered in the gridder registry.
+ static askap::synthesis::IVisGridder::ShPtr makeGridder
+ (const ParameterSet&);
+
+ // @brief Return the (unique) name of the gridder.
+ static const std::string& gridderName();
+
+ // @brief Register the gridder create function with its name.
+ static void registerGridder();
+
+ private:
+ // Correct the visibilities before gridding or after degridding.
+ void correctVisibilities (askap::synthesis::IDataAccessor& acc,
+ bool forward);
+
+ // Initialize the corrections.
+ void initCorrections (const askap::synthesis::IConstDataAccessor& acc);
+
+ // Form an Instrument object.
+ BBS::Instrument makeInstrument (const casa::MeasurementSet& ms);
+
+
+ //# Data members.
+ // Is the gridder correction initialized?
+ bool itsInitialized;
+ // The ParmDB to use.
+ string itsParmDBName;
+ // The config file to use.
+ string itsConfigName;
+ string itsConfigPath;
+ // The StationResponse object to calculate the corrections.
+ BBS::StationResponse::Ptr itsResponse;
+ // Time and frequency interval in MS.
+ double itsTimeInterval;
+ double itsFreqInterval;
+ double itsStartFreq;
+ // The last facet center used.
+ casa::MVDirection itsLastFacetCenter;
+ };
+
+} //# end namespace
+
+#endif
diff --git a/CEP/Imager/LofarGridder/src/CMakeLists.txt b/CEP/Imager/LofarGridder/src/CMakeLists.txt
index c4da0604a08..cffa8da8d90 100644
--- a/CEP/Imager/LofarGridder/src/CMakeLists.txt
+++ b/CEP/Imager/LofarGridder/src/CMakeLists.txt
@@ -3,7 +3,7 @@
include(LofarPackageVersion)
lofar_add_library(lofargridder
- Package__Version.cc Register.cc LofarGridder.cc
+ Package__Version.cc Register.cc LofarWProjectGridder.cc
)
lofar_add_bin_program(versionlofargridder versionlofargridder.cc)
diff --git a/CEP/Imager/LofarGridder/src/LofarGridder.cc b/CEP/Imager/LofarGridder/src/LofarGridder.cc
index a10700cc2f0..cb949611263 100644
--- a/CEP/Imager/LofarGridder/src/LofarGridder.cc
+++ b/CEP/Imager/LofarGridder/src/LofarGridder.cc
@@ -1,4 +1,4 @@
-//# TestGridder.cc: Test visibility gridder.
+//# LofarGridder.cc: Gridder for LOFAR data correcting for DD effects
//#
//# Copyright (C) 2009
//# ASTRON (Netherlands Institute for Radio Astronomy)
@@ -22,31 +22,53 @@
//#
//# @author Ger van Diepen <diepen at astron dot nl>
+#include <lofar_config.h>
+#include <ParmDB/Grid.h>
+//# ASKAP includes
#include <LofarGridder/LofarGridder.h>
#include <gridding/VisGridderFactory.h>
+#include <dataaccess/TableConstDataAccessor.h>
+#include <dataaccess/TableConstDataIterator.h>
+#include <dataaccess/OnDemandBufferDataAccessor.h>
+//#casacore includes
+#include <tables/Tables/TableRecord.h>
+
+using namespace LOFAR::BBS;
+using namespace casa;
+using namespace askap;
using namespace askap::synthesis;
namespace LOFAR
{
- LofarGridder::LofarGridder()
- {}
+ LofarGridder::LofarGridder (const ParameterSet& parset)
+ : itsInitialized (false)
+ {
+ itsTimeAvg = parset.getInt ("average.timestep", 1);
+ itsFreqAvg = parset.getInt ("average.freqstep", 1);
+ itsCorrect = parset.getStringVector ("correct");
+ itsParmDBName = parset.getString ("correct.parmdb");
+ string gridderName = parset.getString ("name");
+ itsIGridder = VisGridderFactory::createGridder (gridderName, parset);
+ itsGridder = dynamic_cast<TableVisGridder*>(itsIGridder.get());
+ ASSERT (itsGridder != 0);
+ }
LofarGridder::~LofarGridder()
{}
// Clone a copy of this Gridder
- IVisGridder::ShPtr LofarGridder::clone()
+ IVisGridder::ShPtr LofarGridder::clone()
{
return IVisGridder::ShPtr(new LofarGridder(*this));
}
- IVisGridder::ShPtr LofarGridder::makeGridder (const ParameterSet&)
+ IVisGridder::ShPtr LofarGridder::makeGridder (const ParameterSet& parset)
{
std::cout << "LofarGridder::makeGridder" << std::endl;
- return IVisGridder::ShPtr(new LofarGridder());
+ return IVisGridder::ShPtr(new LofarGridder(parset));
}
const std::string& LofarGridder::gridderName()
@@ -60,24 +82,103 @@ namespace LOFAR
VisGridderFactory::registerGridder (gridderName(), &makeGridder);
}
- void LofarGridder::initIndices(const IConstDataAccessor&)
+ void LofarGridder::initialiseGrid(const scimath::Axes& axes,
+ const casa::IPosition& shape,
+ const bool dopsf)
{
+ itsGridder->initialiseGrid (axes, shape, dopsf);
}
- void LofarGridder::initConvolutionFunction(const IConstDataAccessor&)
+ void LofarGridder::grid(IConstDataAccessor& acc)
{
- itsSupport=0;
- itsOverSample=1;
- itsCSize=2*(itsSupport+1)*itsOverSample+1; // 3
- itsCCenter=(itsCSize-1)/2; // 1
- itsConvFunc.resize(1);
- itsConvFunc[0].resize(itsCSize, itsCSize); // 3, 3, 1
- itsConvFunc[0].set(0.0);
- itsConvFunc[0](itsCCenter,itsCCenter)=1.0; // 1,1,0 = 1
+ if (! itsInitialized) {
+ initCorrections (acc);
+ }
+ // do correction, uvw rotation, and averaging here
+ // For time being only correction is done.
+ // It is possible to do the following:
+ /// Array<Complex> data (acc.visibility());
+ // which makes a reference copy of the array, thus changs it later.
+ // This is an optimization we won't make yet.
+ // Can use PolConverter::isLinear() to check if data are XX,XY,YX,YY
+ OnDemandBufferDataAccessor acc2(acc);
+ Array<Complex>& data = acc2.rwVisibility();
+ itsGridder->grid (acc2);
+ // Get facet center in J2000.
+ casa::MVDirection center = itsGridder->getImageCentre();
+ double time = acc.time(); // time in sec (do I need exposure?)
+ const Vector<double>& freq = acc.frequency();
+ Axis::ShPtr timeAxis(new RegularAxis ());
+ Axis::ShPtr freqAxis(new RegularAxis ());
}
+
+ void LofarGridder::finaliseGrid(casa::Array<double>& out)
+ {
+ // If averaging is done, it has to grid the last timeslots.
+ itsGridder->finaliseGrid (out);
+ }
+
+ void LofarGridder::finaliseWeights(casa::Array<double>& out)
+ {
+ itsGridder->finaliseWeights (out);
+ }
+
+ void LofarGridder::initialiseDegrid(const scimath::Axes& axes,
+ const casa::Array<double>& image)
+ {
+ itsGridder->initialiseDegrid (axes, image);
+ }
+
+ void LofarGridder::customiseForContext(casa::String context)
+ {
+ itsGridder->customiseForContext (context);
+ }
+
+ void LofarGridder::initVisWeights(IVisWeights::ShPtr viswt)
+ {
+ itsGridder->initVisWeights (viswt);
+ }
+
+ void LofarGridder::degrid(IDataAccessor& acc)
+ {
+ if (! itsInitialized) {
+ initCorrections (acc);
+ }
+ // do correction, uvw rotation, and averaging here
+ // For time being only correction is done.
+ OnDemandBufferDataAccessor acc2(acc);
+ Array<Complex>& data = acc2.rwVisibility();
+ itsGridder->degrid (acc2);
+ }
+
+ void LofarGridder::finaliseDegrid()
+ {
+ itsGridder->finaliseDegrid();
+ }
+
+// void LofarGridder::initIndices(const IConstDataAccessor& acc)
+// {
+// itsGridder->initIndices (acc);
+// }
+
+// void LofarGridder::initConvolutionFunction(const IConstDataAccessor& acc)
+// {
+// itsGridder->initConvolutionFunction (acc);
+// }
- void LofarGridder::correctConvolution(casa::Array<double>& image)
+// void LofarGridder::correctConvolution(casa::Array<double>& image)
+// {
+// itsGridder->correctConvolution (image);
+// }
+
+ void LofarGridder::initCorrections (const IConstDataAccessor& acc)
{
+ const TableConstDataAccessor& tacc =
+ dynamic_cast<const TableConstDataAccessor&>(acc);
+ const TableConstDataIterator& titer = tacc.iterator();
+ Table ms = titer.table();
+ Table antTab(ms.keywordSet().asTable("ANTENNA"));
+ itsInitialized = true;
}
} //# end namespace
diff --git a/CEP/Imager/LofarGridder/src/LofarWProjectGridder.cc b/CEP/Imager/LofarGridder/src/LofarWProjectGridder.cc
new file mode 100644
index 00000000000..f0b68ab3096
--- /dev/null
+++ b/CEP/Imager/LofarGridder/src/LofarWProjectGridder.cc
@@ -0,0 +1,231 @@
+//# LofarWProjectGridder.cc: WProjection gridder for LOFAR data correcting for DD effects
+//#
+//# Copyright (C) 2010
+//# ASTRON (Netherlands Institute for Radio Astronomy)
+//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//#
+//# This file is part of the LOFAR software suite.
+//# The LOFAR software suite is free software: you can redistribute it and/or
+//# modify it under the terms of the GNU General Public License as published
+//# by the Free Software Foundation, either version 3 of the License, or
+//# (at your option) any later version.
+//#
+//# The LOFAR software suite is distributed in the hope that it will be useful,
+//# but WITHOUT ANY WARRANTY; without even the implied warranty of
+//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+//# GNU General Public License for more details.
+//#
+//# You should have received a copy of the GNU General Public License along
+//# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//#
+//# $Id$
+//#
+//# @author Ger van Diepen <diepen at astron dot nl>
+
+#include <lofar_config.h>
+#include <ParmDB/Grid.h>
+
+//# ASKAP includes
+#include <LofarGridder/LofarWProjectGridder.h>
+#include <gridding/VisGridderFactory.h>
+#include <dataaccess/TableConstDataAccessor.h>
+#include <dataaccess/TableConstDataIterator.h>
+
+//# casacore includes
+#include <tables/Tables/TableRecord.h>
+#include <ms/MeasurementSets/MSColumns.h>
+#include <ms/MeasurementSets/MSSpWindowColumns.h>
+#include <ms/MeasurementSets/MSAntennaColumns.h>
+#include <ms/MeasurementSets/MSObsColumns.h>
+#include <measures/Measures/MeasTable.h>
+#include <casa/OS/Path.h>
+
+using namespace LOFAR::BBS;
+using namespace casa;
+using namespace askap;
+using namespace askap::synthesis;
+
+namespace LOFAR
+{
+
+ LofarWProjectGridder::LofarWProjectGridder (const ParameterSet& parset)
+ : WProjectVisGridder (*(dynamic_cast<const WProjectVisGridder*>
+ (WProjectVisGridder::createGridder(parset).get()))),
+ itsInitialized (false)
+ {
+ // Currently only beam can be corrected for.
+ /// itsTimeAvg = parset.getInt ("average.timestep", 1);
+ /// itsFreqAvg = parset.getInt ("average.freqstep", 1);
+ /// itsCorrect = parset.getStringVector ("correct");
+ itsParmDBName = parset.getString ("ParmDB.Instrument", "");
+ itsConfigName = parset.getString ("Beam.StationConfig.Name");
+ itsConfigPath = parset.getString ("Beam.StationConfig.Path");
+ }
+
+ LofarWProjectGridder::~LofarWProjectGridder()
+ {}
+
+ // Clone a copy of this Gridder
+ IVisGridder::ShPtr LofarWProjectGridder::clone()
+ {
+ return IVisGridder::ShPtr(new LofarWProjectGridder(*this));
+ }
+
+ IVisGridder::ShPtr LofarWProjectGridder::makeGridder
+ (const ParameterSet& parset)
+ {
+ std::cout << "LofarWProjectGridder::makeGridder" << std::endl;
+ return IVisGridder::ShPtr(new LofarWProjectGridder(parset));
+ }
+
+ const std::string& LofarWProjectGridder::gridderName()
+ {
+ static std::string name("LofarGridder.WProject");
+ return name;
+ }
+
+ void LofarWProjectGridder::registerGridder()
+ {
+ VisGridderFactory::registerGridder (gridderName(), &makeGridder);
+ }
+
+ void LofarWProjectGridder::correctVisibilities (IDataAccessor& acc,
+ bool forward)
+ {
+ if (!itsInitialized) {
+ initCorrections (acc);
+ }
+ // If changed, set the facet center as direction of interest.
+ MVDirection facetCenter = getImageCentre();
+ if (facetCenter != itsLastFacetCenter) {
+ itsLastFacetCenter = facetCenter;
+ itsResponse->setDirection (MDirection(facetCenter, MDirection::J2000));
+ }
+ // Get read/write access to the visibilities. It makes a copy if needed.
+ Array<Complex>& visData = acc.rwVisibility();
+ ASSERT (visData.contiguousStorage());
+ Complex* data = visData.data();
+ uInt nfreq = visData.shape()[1];
+ double time = acc.time();
+ // Set the axes for the BBS request.
+ Axis::ShPtr freqAxis
+ (new RegularAxis(itsStartFreq-itsFreqInterval/2, itsFreqInterval, nfreq));
+ Axis::ShPtr timeAxis
+ (new RegularAxis(time-itsTimeInterval/2, itsTimeInterval, 1));
+ Grid grid(timeAxis, freqAxis);
+ itsResponse->setEvalGrid (grid);
+ const casa::Vector<uInt>& ant1 = acc.antenna1();
+ const casa::Vector<uInt>& ant2 = acc.antenna2();
+ for (uInt i=0; i<ant1.size(); ++i) {
+ JonesMatrix::View j1 = itsResponse->evaluate (ant1[i]);
+ JonesMatrix::View j2 = itsResponse->evaluate (ant2[i]);
+ for (uInt j=0; j<nfreq; ++j) {
+ dcomplex cl0 = j1(j,0).getDComplex(0,0);
+ dcomplex cl1 = j1(j,0).getDComplex(1,0);
+ dcomplex cl2 = j1(j,0).getDComplex(0,1);
+ dcomplex cl3 = j1(j,0).getDComplex(1,1);
+ dcomplex cr0 = conj (j2(j,0).getDComplex(0,0));
+ dcomplex cr1 = conj (j2(j,0).getDComplex(1,0));
+ dcomplex cr2 = conj (j2(j,0).getDComplex(0,1));
+ dcomplex cr3 = conj (j2(j,0).getDComplex(1,1));
+ dcomplex d0 = data[0];
+ dcomplex d1 = data[1];
+ dcomplex d2 = data[2];
+ dcomplex d3 = data[3];
+ if (forward) {
+ // For gridding use inverse of responses.
+ // data = inv(res(i)) * data * inv(conjtranspose(res(j)))
+ // Note that inv(conjtranspose(A)) == conjtranspose(inv(A))
+ // inv(a b) == ( d -b) / (ad-bc)
+ // (c d) (-c a)
+ dcomplex tmp0 = cl3*d0 - cl1*d2;
+ dcomplex tmp1 = cl3*d1 - cl1*d3;
+ dcomplex tmp2 = cl0*d2 - cl2*d0;
+ dcomplex tmp3 = cl0*d3 - cl2*d1;
+ dcomplex factor = 1. / ((cl0*cl3 - cl1*cl2) * (cr0*cr3 - cr1*cr2));
+ data[0] = factor * (tmp0*cr3 - tmp1*cr1);
+ data[1] = factor * (tmp1*cr0 - tmp0*cr2);
+ data[2] = factor * (tmp2*cr3 - tmp3*cr1);
+ data[3] = factor * (tmp3*cr0 - tmp2*cr2);
+ } else {
+ // Degridding; data = res(i) * data * conjtranspose(res(j))
+ dcomplex tmp0 = cl0*d0 + cl1*d2;
+ dcomplex tmp1 = cl0*d1 + cl1*d3;
+ dcomplex tmp2 = cl2*d0 + cl3*d2;
+ dcomplex tmp3 = cl2*d1 + cl3*d3;
+ data[0] = tmp0*cr0 + tmp1*cr1;
+ data[1] = tmp0*cr2 + tmp1*cr3;
+ data[2] = tmp2*cr0 + tmp3*cr1;
+ data[3] = tmp2*cr2 + tmp3*cr3;
+ }
+ data += 4;
+ }
+ }
+ }
+
+ void LofarWProjectGridder::initCorrections (const IConstDataAccessor& acc)
+ {
+ // Ensure there are polarizations XX,XY,YX,YY.
+ ASSERT (acc.stokes().size() == 4);
+ ASSERT (acc.stokes()[0] == Stokes::XX);
+ ASSERT (acc.stokes()[1] == Stokes::XY);
+ ASSERT (acc.stokes()[2] == Stokes::YX);
+ ASSERT (acc.stokes()[3] == Stokes::YY);
+ const TableConstDataAccessor& tacc =
+ dynamic_cast<const TableConstDataAccessor&>(acc);
+ const TableConstDataIterator& titer = tacc.iterator();
+ MeasurementSet ms(titer.table());
+ // Get the time interval.
+ ROMSColumns msCol(ms);
+ // Get the reference frequency and width.
+ // It assumes that the interval is constant for this band.
+ ROMSSpWindowColumns spwCols(ms.spectralWindow());
+ double refFreq = spwCols.refFrequency()(0);
+ itsFreqInterval = spwCols.chanWidth()(0).data()[0];
+ casa::Vector<double> freqs = spwCols.chanFreq()(0);
+ itsStartFreq = freqs[0];
+ for (uInt i=1; i<freqs.size(); ++i) {
+ ASSERTSTR (casa::near(freqs[i], itsStartFreq + i*itsFreqInterval),
+ "Frequency channels are not regularly spaced");
+ }
+ // Form the StationResponse object.
+ itsResponse = StationResponse::Ptr (new StationResponse
+ (makeInstrument(ms),
+ itsConfigName, Path(itsConfigPath),
+ refFreq));
+ // Set the pointing direction (for beamforming).
+ // Use first value of MDirection array in first row in FIELD subtable.
+ ROMSFieldColumns fieldCols(ms.field());
+ itsResponse->setPointing (fieldCols.delayDirMeasCol()(0).data()[0]);
+ itsTimeInterval = msCol.interval()(0);
+ // Initialization is done.
+ itsInitialized = true;
+ }
+
+ BBS::Instrument LofarWProjectGridder::makeInstrument
+ (const MeasurementSet& ms)
+ {
+ // Get all station positions.
+ ROMSAntennaColumns antCols(ms.antenna());
+ vector<Station> stations;
+ stations.reserve (ms.antenna().nrow());
+ for (uInt i=0; i<ms.antenna().nrow(); ++i) {
+ stations.push_back (Station(antCols.name()(i),
+ antCols.positionMeas()(i)));
+ }
+ // Find observatory position.
+ // If not found, set it to the position of the middle station.
+ MPosition arrayPos;
+ Bool fndObs = False;
+ if (ms.observation().nrow() > 0) {
+ ROMSObservationColumns obsCols(ms.observation());
+ String telescope = obsCols.telescopeName()(0);
+ fndObs = MeasTable::Observatory (arrayPos, telescope);
+ }
+ if (!fndObs && stations.size() > 0) {
+ arrayPos = stations[stations.size()/2].position();
+ }
+ return BBS::Instrument ("LOFAR", arrayPos, stations);
+ }
+
+} //# end namespace
diff --git a/CEP/Imager/LofarGridder/src/Register.cc b/CEP/Imager/LofarGridder/src/Register.cc
index 51b68814d2c..a8c1fd99429 100644
--- a/CEP/Imager/LofarGridder/src/Register.cc
+++ b/CEP/Imager/LofarGridder/src/Register.cc
@@ -23,10 +23,11 @@
//# @author Ger van Diepen <diepen at astron dot nl>
//# Includes
+#include <lofar_config.h>
#include <LofarGridder/Register.h>
-#include <LofarGridder/LofarGridder.h>
+#include <LofarGridder/LofarWProjectGridder.h>
void register_lofargridder()
{
- LOFAR::LofarGridder::registerGridder();
+ LOFAR::LofarWProjectGridder::registerGridder();
}
diff --git a/CEP/Imager/LofarGridder/test/tLofarGridder.cc b/CEP/Imager/LofarGridder/test/tLofarGridder.cc
index 23e3cb68c94..1abe5b717b1 100644
--- a/CEP/Imager/LofarGridder/test/tLofarGridder.cc
+++ b/CEP/Imager/LofarGridder/test/tLofarGridder.cc
@@ -37,7 +37,7 @@ int main (int argc, const char** argv)
// Create a LofarGridder instance.
ParameterSet parset;
- parset.add ("gridder", "LofarGridder");
+ parset.add ("gridder", "LofarGridder.WProject");
IVisGridder::ShPtr gridder = VisGridderFactory::make (parset);
} catch (const std::exception& x) {
--
GitLab