From b8b818e9c084505c7b5bf5d9f7fd802e3c8a92cb Mon Sep 17 00:00:00 2001 From: Joris van Zwieten <zwieten@astron.nl> Date: Thu, 28 Feb 2008 12:14:21 +0000 Subject: [PATCH] Bug 1163: - Refactored BaseLinPS and BaseDFTPS into a proper phase shift and a point source coherence. - Added gaussian source coherence. --- .../BBSKernel/MNS/MeqGaussianCoherency.h | 85 ++++++ .../{MeqGaussSource.h => MeqGaussianSource.h} | 32 +-- .../include/BBSKernel/MNS/MeqJonesCMul2.h | 4 +- .../include/BBSKernel/MNS/MeqJonesCMul3.h | 4 +- .../include/BBSKernel/MNS/MeqJonesMul.h | 68 +++++ .../MNS/{MeqBaseDFTPS.h => MeqPhaseShift.h} | 32 +-- .../{MeqBaseLinPS.h => MeqPointCoherency.h} | 20 +- .../include/BBSKernel/MNS/MeqPointSource.h | 8 +- .../BBSKernel/include/BBSKernel/Makefile.am | 8 +- CEP/BB/BBSKernel/src/MNS/MeqBaseDFTPS.cc | 198 ------------- CEP/BB/BBSKernel/src/MNS/MeqBaseLinPS.cc | 137 --------- .../BBSKernel/src/MNS/MeqGaussianCoherency.cc | 267 ++++++++++++++++++ ...MeqGaussSource.cc => MeqGaussianSource.cc} | 44 +-- CEP/BB/BBSKernel/src/MNS/MeqJonesMul.cc | 159 +++++++++++ CEP/BB/BBSKernel/src/MNS/MeqPhaseShift.cc | 219 ++++++++++++++ CEP/BB/BBSKernel/src/MNS/MeqPointCoherency.cc | 149 ++++++++++ CEP/BB/BBSKernel/src/MNS/MeqSourceList.cc | 8 +- CEP/BB/BBSKernel/src/Makefile.am | 8 +- CEP/BB/BBSKernel/src/Model.cc | 28 +- 19 files changed, 1055 insertions(+), 423 deletions(-) create mode 100644 CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqGaussianCoherency.h rename CEP/BB/BBSKernel/include/BBSKernel/MNS/{MeqGaussSource.h => MeqGaussianSource.h} (78%) create mode 100644 CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqJonesMul.h rename CEP/BB/BBSKernel/include/BBSKernel/MNS/{MeqBaseDFTPS.h => MeqPhaseShift.h} (66%) rename CEP/BB/BBSKernel/include/BBSKernel/MNS/{MeqBaseLinPS.h => MeqPointCoherency.h} (75%) delete mode 100644 CEP/BB/BBSKernel/src/MNS/MeqBaseDFTPS.cc delete mode 100644 CEP/BB/BBSKernel/src/MNS/MeqBaseLinPS.cc create mode 100644 CEP/BB/BBSKernel/src/MNS/MeqGaussianCoherency.cc rename CEP/BB/BBSKernel/src/MNS/{MeqGaussSource.cc => MeqGaussianSource.cc} (59%) create mode 100644 CEP/BB/BBSKernel/src/MNS/MeqJonesMul.cc create mode 100644 CEP/BB/BBSKernel/src/MNS/MeqPhaseShift.cc create mode 100644 CEP/BB/BBSKernel/src/MNS/MeqPointCoherency.cc diff --git a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqGaussianCoherency.h b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqGaussianCoherency.h new file mode 100644 index 00000000000..da27c806cd9 --- /dev/null +++ b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqGaussianCoherency.h @@ -0,0 +1,85 @@ +//# MeqGaussianCoherency.h: Spatial coherence function of an elliptical +//# gaussian source. +//# +//# Copyright (C) 2005 +//# ASTRON (Netherlands Foundation for Research in Astronomy) +//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands, seg@astron.nl +//# +//# This program 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 2 of the License, or +//# (at your option) any later version. +//# +//# This program 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 this program; if not, write to the Free Software +//# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +//# +//# $Id$ + +#ifndef MNS_MEQGAUSSIANCOHERENCY_H +#define MNS_MEQGAUSSIANCOHENRECY_H + +// \file +// Spatial coherence function of an elliptical gaussian source. + +//# Includes +#include <BBSKernel/MNS/MeqJonesExpr.h> +#include <BBSKernel/MNS/MeqGaussianSource.h> +#include <BBSKernel/MNS/MeqStatUVW.h> + +namespace LOFAR +{ +namespace BBS +{ + +// \ingroup BBSKernel +// \addtogroup MNS +// @{ + +//# Forward Declarations + + +class MeqGaussianCoherency: public MeqJonesExprRep +{ +public: + enum + { + IN_I, + IN_Q, + IN_U, + IN_V, + IN_MAJOR, + IN_MINOR, + IN_PHI, + N_InputPort + } InputPort; + + MeqGaussianCoherency(const MeqGaussianSource *source, MeqStatUVW *station1, MeqStatUVW *station2); + ~MeqGaussianCoherency(); + + // Calculate the results for the given domain. + virtual MeqJonesResult getJResult(const MeqRequest &request); + +private: +#ifdef EXPR_GRAPH + virtual std::string getLabel(); +#endif + + MeqMatrix computeCoherence(const MeqRequest &request, + const MeqMatrix &uBaseline, const MeqMatrix &vBaseline, + const MeqMatrix &major, const MeqMatrix &minor, const MeqMatrix &phi); + + const MeqGaussianSource *itsSource; + MeqStatUVW *itsStation1, *itsStation2; +}; + +// @} + +} // namespace BBS +} // namespace LOFAR +#endif diff --git a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqGaussSource.h b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqGaussianSource.h similarity index 78% rename from CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqGaussSource.h rename to CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqGaussianSource.h index 2baa4c9ce5a..6d177397921 100644 --- a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqGaussSource.h +++ b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqGaussianSource.h @@ -1,4 +1,4 @@ -//# MeqGaussSource.h: Class holding a gaussian source +//# MeqGaussianSource.h: Class holding a gaussian source //# //# Copyright (C) 2002 //# ASTRON (Netherlands Foundation for Research in Astronomy) @@ -20,8 +20,8 @@ //# //# $Id$ -#ifndef MNS_MEQGAUSSSOURCE_H -#define MNS_MEQGAUSSSOURCE_H +#ifndef MNS_MEQGAUSSIANSOURCE_H +#define MNS_MEQGAUSSIANSOURCE_H // \file // Class holding a gaussian source @@ -40,33 +40,33 @@ namespace BBS // @{ -class MeqGaussSource: public MeqSource +class MeqGaussianSource: public MeqSource { public: - MeqGaussSource (const string& name, + MeqGaussianSource(const string& name, const string& groupName, const MeqExpr& fluxI, const MeqExpr& fluxQ, const MeqExpr& fluxU, const MeqExpr& fluxV, const MeqExpr& ra, const MeqExpr& dec, - const MeqExpr& minor, const MeqExpr& major, + const MeqExpr& major, const MeqExpr& minor, const MeqExpr& phi); - virtual ~MeqGaussSource(); + virtual ~MeqGaussianSource(); - MeqExpr& getI() + MeqExpr getI() const { return itsI; } - MeqExpr& getQ() + MeqExpr getQ() const { return itsQ; } - MeqExpr& getU() + MeqExpr getU() const { return itsU; } - MeqExpr& getV() + MeqExpr getV() const { return itsV; } - MeqExpr& getMinor() - { return itsMinor; } - MeqExpr& getMajor() + MeqExpr getMajor() const { return itsMajor; } - MeqExpr& getPhi() + MeqExpr getMinor() const + { return itsMinor; } + MeqExpr getPhi() const { return itsPhi; } private: @@ -74,8 +74,8 @@ private: MeqExpr itsQ; MeqExpr itsU; MeqExpr itsV; - MeqExpr itsMinor; MeqExpr itsMajor; + MeqExpr itsMinor; MeqExpr itsPhi; }; diff --git a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqJonesCMul2.h b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqJonesCMul2.h index 60aab4f5902..a8fa5caa402 100644 --- a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqJonesCMul2.h +++ b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqJonesCMul2.h @@ -23,7 +23,7 @@ #define MNS_MEQJONESCMUL2_H // \file -// Calculate left*conj(right) +// Calculate left*conj(transpose(right)) //# Includes #include <BBSKernel/MNS/MeqJonesExpr.h> @@ -38,7 +38,7 @@ namespace BBS // @{ -// Calculate left*conj(right). +// Calculate left*conj(transpose(right)). class MeqJonesCMul2: public MeqJonesExprRep { diff --git a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqJonesCMul3.h b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqJonesCMul3.h index f45cf222f31..79f513e19da 100644 --- a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqJonesCMul3.h +++ b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqJonesCMul3.h @@ -24,7 +24,7 @@ #define MNS_MEQJONESCMUL3_H // \file -// Calculate left*mid*conj(right) +// Calculate left*mid*conj(transpose(right)) //# Includes #include <BBSKernel/MNS/MeqJonesExpr.h> @@ -39,7 +39,7 @@ namespace BBS // @{ -// Calculate left*mid*conj(right). +// Calculate left*mid*conj(transpose(right)). class MeqJonesCMul3: public MeqJonesExprRep { diff --git a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqJonesMul.h b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqJonesMul.h new file mode 100644 index 00000000000..592e7b7f82f --- /dev/null +++ b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqJonesMul.h @@ -0,0 +1,68 @@ +//# MeqJonesMul.h: Multiply each component of a MeqJonesExpr with a single +//# MeqExpr. +//# +//# Copyright (C) 2007 +//# ASTRON (Netherlands Foundation for Research in Astronomy) +//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands, seg@astron.nl +//# +//# This program 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 2 of the License, or +//# (at your option) any later version. +//# +//# This program 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 this program; if not, write to the Free Software +//# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +//# +//# $Id$ + +#ifndef MNS_MEQJONESMUL_H +#define MNS_MEQJONESMUL_H + +// \file Multiply each component of a MeqJonesExpr with a single MeqExpr. + +#include <BBSKernel/MNS/MeqExpr.h> +#include <BBSKernel/MNS/MeqJonesExpr.h> +#include <BBSKernel/MNS/MeqJonesResult.h> + +namespace LOFAR +{ +namespace BBS +{ + +// \ingroup BBSKernel +// \addtogroup MNS +// @{ + + +// Multiply each component of a MeqJonesExpr with a single MeqExpr. + +class MeqJonesMul: public MeqJonesExprRep +{ +public: + MeqJonesMul(const MeqJonesExpr &left, const MeqExpr &right); + ~MeqJonesMul(); + + // Get the result of the expression for the given domain. + MeqJonesResult getJResult(const MeqRequest &request); + +private: +#ifdef EXPR_GRAPH + virtual std::string getLabel(); +#endif + + MeqJonesExpr itsLeft; + MeqExpr itsRight; +}; + +// @} + +} // namespace BBS +} // namespace LOFAR + +#endif diff --git a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqBaseDFTPS.h b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqPhaseShift.h similarity index 66% rename from CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqBaseDFTPS.h rename to CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqPhaseShift.h index f0eeff6c04f..0670ab0afd2 100644 --- a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqBaseDFTPS.h +++ b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqPhaseShift.h @@ -1,4 +1,5 @@ -//# MeqBaseDFTPS.h: Baseline DFT of a point source +//# MeqPhaseShift.h: Phase delay due to baseline geometry with respect to source +//# direction. //# //# Copyright (C) 2005 //# ASTRON (Netherlands Foundation for Research in Astronomy) @@ -20,11 +21,11 @@ //# //# $Id$ -#ifndef MNS_MEQBASEDFTPS_H -#define MNS_MEQBASEDFTPS_H +#ifndef MNS_MEQPHASESHIFT_H +#define MNS_MEQPHASESHIFT_H // \file -// Baseline DFT of a point source. +// Phase delay due to baseline geometry with respect to source direction. //# Includes #include <BBSKernel/MNS/MeqExpr.h> @@ -38,25 +39,22 @@ namespace BBS // \addtogroup MNS // @{ -class MeqBaseDFTPS: public MeqExprRep +class MeqPhaseShift: public MeqExprRep { public: - MeqBaseDFTPS (const MeqExpr& left, const MeqExpr& right, - const MeqExpr& lmn); + MeqPhaseShift (const MeqExpr& left, const MeqExpr& right); + ~MeqPhaseShift(); - ~MeqBaseDFTPS(); - - // Calculate the results for the given domain. - virtual MeqResult getResult (const MeqRequest&); + // Calculate the results for the given domain. + virtual MeqResult getResult(const MeqRequest &request); private: -#ifdef EXPR_GRAPH - virtual std::string getLabel(); -#endif + #ifdef EXPR_GRAPH + virtual std::string getLabel(); + #endif - MeqExpr itsLeft; - MeqExpr itsRight; - MeqExpr itsLMN; + MeqExpr itsLeft; + MeqExpr itsRight; }; // @} diff --git a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqBaseLinPS.h b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqPointCoherency.h similarity index 75% rename from CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqBaseLinPS.h rename to CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqPointCoherency.h index 650a721a7f0..4d71c18d429 100644 --- a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqBaseLinPS.h +++ b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqPointCoherency.h @@ -1,4 +1,4 @@ -//# MeqBaseLinPS.h: Baseline prediction of a point source with linear polarisation +//# MeqPointCoherency.h: Spatial coherence function of a point source. //# //# Copyright (C) 2005 //# ASTRON (Netherlands Foundation for Research in Astronomy) @@ -20,11 +20,11 @@ //# //# $Id$ -#ifndef MNS_MEQBASELINPS_H -#define MNS_MEQBASELINPS_H +#ifndef MNS_MEQPOINTCOHERENCY_H +#define MNS_MEQPOINTCOHENRECY_H // \file -// Baseline prediction of a point source with linear polarisation. +// Spatial coherence function of a point source. //# Includes #include <BBSKernel/MNS/MeqJonesExpr.h> @@ -44,23 +44,21 @@ namespace BBS //# Forward Declarations -class MeqBaseLinPS: public MeqJonesExprRep +class MeqPointCoherency: public MeqJonesExprRep { public: - MeqBaseLinPS (const MeqExpr& dft, MeqPointSource* src); - - ~MeqBaseLinPS(); + MeqPointCoherency(const MeqPointSource *source); + ~MeqPointCoherency(); // Calculate the results for the given domain. - virtual MeqJonesResult getJResult (const MeqRequest&); + virtual MeqJonesResult getJResult(const MeqRequest &request); private: #ifdef EXPR_GRAPH virtual std::string getLabel(); #endif - MeqExpr itsDFT; - MeqPointSource* itsSource; + const MeqPointSource *itsSource; }; // @} diff --git a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqPointSource.h b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqPointSource.h index d9b527ce1a8..fccfc9fd5aa 100644 --- a/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqPointSource.h +++ b/CEP/BB/BBSKernel/include/BBSKernel/MNS/MeqPointSource.h @@ -51,13 +51,13 @@ public: virtual ~MeqPointSource(); - MeqExpr& getI() + MeqExpr getI() const { return itsI; } - MeqExpr& getQ() + MeqExpr getQ() const { return itsQ; } - MeqExpr& getU() + MeqExpr getU() const { return itsU; } - MeqExpr& getV() + MeqExpr getV() const { return itsV; } private: diff --git a/CEP/BB/BBSKernel/include/BBSKernel/Makefile.am b/CEP/BB/BBSKernel/include/BBSKernel/Makefile.am index 54e2fffcd54..9c1052d4695 100644 --- a/CEP/BB/BBSKernel/include/BBSKernel/Makefile.am +++ b/CEP/BB/BBSKernel/include/BBSKernel/Makefile.am @@ -15,19 +15,19 @@ pkginclude_HEADERS = Package__Version.h \ pkginclude_MNS_HEADERS = \ MNS/MeqAzEl.h \ - MNS/MeqBaseDFTPS.h \ - MNS/MeqBaseLinPS.h \ MNS/MeqDFTPS.h \ MNS/MeqDiag.h \ MNS/MeqDipoleBeam.h \ MNS/MeqDomain.h \ MNS/MeqExpr.h \ MNS/MeqFunklet.h \ - MNS/MeqGaussSource.h \ + MNS/MeqGaussianCoherency.h \ + MNS/MeqGaussianSource.h \ MNS/MeqJonesCMul2.h \ MNS/MeqJonesCMul3.h \ MNS/MeqJonesExpr.h \ MNS/MeqJonesInvert.h \ + MNS/MeqJonesMul.h \ MNS/MeqJonesMul2.h \ MNS/MeqJonesMul3.h \ MNS/MeqJonesNode.h \ @@ -47,6 +47,8 @@ pkginclude_MNS_HEADERS = \ MNS/MeqParm.h \ MNS/MeqParmSingle.h \ MNS/MeqPhaseRef.h \ + MNS/MeqPhaseShift.h \ + MNS/MeqPointCoherency.h \ MNS/MeqPointSource.h \ MNS/MeqPolc.h \ MNS/MeqPolcLog.h \ diff --git a/CEP/BB/BBSKernel/src/MNS/MeqBaseDFTPS.cc b/CEP/BB/BBSKernel/src/MNS/MeqBaseDFTPS.cc deleted file mode 100644 index dc78008254d..00000000000 --- a/CEP/BB/BBSKernel/src/MNS/MeqBaseDFTPS.cc +++ /dev/null @@ -1,198 +0,0 @@ -//# MeqBaseDFTPS.cc: Baseline prediction of a point source with linear polarisation -//# -//# Copyright (C) 2005 -//# ASTRON (Netherlands Foundation for Research in Astronomy) -//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands, seg@astron.nl -//# -//# This program 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 2 of the License, or -//# (at your option) any later version. -//# -//# This program 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 this program; if not, write to the Free Software -//# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -//# -//# $Id$ - -#include <lofar_config.h> - -#include <BBSKernel/MNS/MeqBaseDFTPS.h> -#include <BBSKernel/MNS/MeqDFTPS.h> -#include <BBSKernel/MNS/MeqStatUVW.h> -#include <BBSKernel/MNS/MeqMatrixTmp.h> -#include <Common/LofarLogger.h> -#include <Common/lofar_iomanip.h> - -using namespace casa; - -namespace LOFAR -{ -namespace BBS -{ -using LOFAR::dcomplex; -using LOFAR::conj; - -MeqBaseDFTPS::MeqBaseDFTPS (const MeqExpr& left, const MeqExpr& right, - const MeqExpr& lmn) - : itsLeft (left), - itsRight (right), - itsLMN (lmn) -{ - addChild (itsLeft); - addChild (itsRight); - addChild (itsLMN); -} - -MeqBaseDFTPS::~MeqBaseDFTPS() -{} - - -MeqResult MeqBaseDFTPS::getResult (const MeqRequest& request) -{ - MeqResult result(request.nspid()); - - // Get N (for the division). - // Assert it is a scalar value. - MeqResultVec lmnRes; - const MeqResultVec& lmn = itsLMN.getResultVecSynced (request, lmnRes); - const MeqResult& nk = lmn[2]; - DBGASSERT (nk.getValue().nelements() == request.ny() || - nk.getValue().nelements() == 1); - // Calculate the left and right station Jones matrix elements. - // A delta in the station source predict is only available if multiple - // frequency channels are used. - bool multFreq = request.nx() > 1; - MeqResultVec reslBuf, resrBuf; - const MeqResultVec& resl = itsLeft.getResultVecSynced (request, resrBuf); - const MeqResultVec& resr = itsRight.getResultVecSynced (request, reslBuf); - const MeqResult& left = resl[0]; - const MeqResult& right = resr[0]; - const MeqResult& leftDelta = resl[1]; - const MeqResult& rightDelta = resr[1]; - - // It is tried to compute the DFT as efficient as possible. - // Therefore the baseline contribution is split into its antenna parts. - // dft = exp(2i.pi(ul+vm+wn)) / n (with u,v,w in wavelengths) - // = (exp(2i.pi((u1.l+v1.m+w1.n) - (u2.l+v2.m+w2.n))/wvl))/n (u,v,w in m) - // = ((exp(i(u1.l+v1.m+w1.n)) / exp(i(u2.l+v2.m+w2.m))) ^ (2.pi/wvl)) / n - // So left and right return the exp values independent of wavelength. - // Thereafter they are scaled to the freq domain by raising the values - // for each time to the appropriate powers. - // Alas the rule - // x^(a*b) = (x^a)^b - // which is valid for real numbers, is only valid for complex numbers - // if b is an integer number. - // Therefore the station calculations (in MeqStatSources) are done as - // follows, where it should be noted that the frequencies are regularly - // strided. - // f = f0 + k.df (k = 0 ... nchan-1) - // s1 = (u1.l+v1.m+w1.n).2i.pi/c - // s2 = (u2.l+v2.m+w2.n).2i.pi/c - // dft = exp(s1(f0+k.df)) / exp(s2(f0+k.df)) / n - // = (exp(s1.f0)/exp(s2.f0)) . (exp(s1.k.df)/exp(s2.k.df)) / n - // = (exp(s1.f0)/exp(s2.f0)) . (exp(s1.df)/exp(s2.df))^k / n - // In principle the power is expensive, but because the frequencies are - // regularly strided, it is possible to use multiplication. - // So it gets - // dft(f0) = (exp(s1.f0)/exp(s2.f0)) / n - // dft(fj) = dft(fj-1) * (exp(s1.df)/exp(s2.df)) - // Using a python script (tuvw.py) is is checked that this way of - // calculation is accurate enough. - // Another optimization can be achieved in the division of the two - // complex numbers which can be turned into a cheaper multiplication. - // exp(x)/exp(y) = (cos(x) + i.sin(x)) / (cos(y) + i.sin(y)) - // = (cos(x) + i.sin(x)) * (cos(y) - i.sin(y)) -/* { - - MeqDFTPS *dftpsl = dynamic_cast<MeqDFTPS*>(itsLeft.rep()); - MeqDFTPS *dftpsr = dynamic_cast<MeqDFTPS*>(itsRight.rep()); - - const MeqResult& ul = dftpsl->uvw()->getU(request); - const MeqResult& vl = dftpsl->uvw()->getV(request); - const MeqResult& wl = dftpsl->uvw()->getW(request); - const MeqResult& ur = dftpsr->uvw()->getU(request); - const MeqResult& vr = dftpsr->uvw()->getV(request); - const MeqResult& wr = dftpsr->uvw()->getW(request); - cout << "U: " << setprecision(20) << ur.getValue() - ul.getValue() - << endl; - cout << "V: " << setprecision(20) << vr.getValue() - vl.getValue() - << endl; - cout << "W: " << setprecision(20) << wr.getValue() - wl.getValue() - << endl;*/ -// LOG_TRACE_FLOW ("U: " << ur.getValue() - ul.getValue()); -// LOG_TRACE_FLOW ("V: " << vr.getValue() - vl.getValue()); -// LOG_TRACE_FLOW ("W: " << wr.getValue() - wl.getValue()); -// } - - MeqMatrix res(makedcomplex(0,0), request.nx(), request.ny(), false); - for (int iy=0; iy<request.ny(); ++iy) { - dcomplex tmpl = left.getValue().getDComplex(0,iy); - dcomplex tmpr = right.getValue().getDComplex(0,iy); - // We have to divide by N. - // However, we divide by 2N to get the factor 0.5 needed in (I+Q)/2, etc. - // in MeqBaseLinPS. -// double tmpnk = 2. * nk.getValue().getDouble(0,iy); - double tmpnk = 2.0; - dcomplex factor; - if (multFreq) { - dcomplex deltal = leftDelta.getValue().getDComplex(0,iy); - dcomplex deltar = rightDelta.getValue().getDComplex(0,iy); - factor = deltar * conj(deltal); - } - res.fillRowWithProducts (tmpr * conj(tmpl) / tmpnk, factor, iy); - } - result.setValue (res); - -// cout << "DFT:" << endl; -// cout << setprecision(20) << res << endl; - - // Evaluate (if needed) for the perturbed parameter values. - // Note that we do not have to test for perturbed values in nk, - // because the left and right value already depend on nk. - const MeqParmFunklet* perturbedParm; - for (int spinx=0; spinx<request.nspid(); spinx++) { - bool eval = false; - if (left.isDefined(spinx)) { - eval = true; - perturbedParm = left.getPerturbedParm(spinx); - } else if (right.isDefined(spinx)) { - eval = true; - perturbedParm = right.getPerturbedParm(spinx); - } - if (eval) { - MeqMatrix pres(makedcomplex(0,0), request.nx(), request.ny(), false); - for (int iy=0; iy<request.ny(); ++iy) { - dcomplex tmpl = left.getPerturbedValue(spinx).getDComplex(0,iy); - dcomplex tmpr = right.getPerturbedValue(spinx).getDComplex(0,iy); -// double tmpnk = 2. * nk.getPerturbedValue(spinx).getDouble(0,iy); - double tmpnk = 2.0; - dcomplex factor; - if (multFreq) { - dcomplex deltal = leftDelta.getPerturbedValue(spinx).getDComplex(0,iy); - dcomplex deltar = rightDelta.getPerturbedValue(spinx).getDComplex(0,iy); - factor = deltar * conj(deltal); - } - pres.fillRowWithProducts(tmpr * conj(tmpl) / tmpnk, factor, iy); - result.setPerturbedValue (spinx, pres); - result.setPerturbedParm (spinx, perturbedParm); - } - } - } - return result; -} - -#ifdef EXPR_GRAPH -std::string MeqBaseDFTPS::getLabel() -{ - return std::string("MeqBaseDFTPS\\nBaseline DFT of a point source"); -} -#endif - -} // namespace BBS -} // namespace LOFAR diff --git a/CEP/BB/BBSKernel/src/MNS/MeqBaseLinPS.cc b/CEP/BB/BBSKernel/src/MNS/MeqBaseLinPS.cc deleted file mode 100644 index 214e0109df9..00000000000 --- a/CEP/BB/BBSKernel/src/MNS/MeqBaseLinPS.cc +++ /dev/null @@ -1,137 +0,0 @@ -//# MeqBaseLinPS.cc: Baseline prediction of a point source with linear polarisation -//# -//# Copyright (C) 2005 -//# ASTRON (Netherlands Foundation for Research in Astronomy) -//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands, seg@astron.nl -//# -//# This program 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 2 of the License, or -//# (at your option) any later version. -//# -//# This program 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 this program; if not, write to the Free Software -//# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -//# -//# $Id$ - -#include <lofar_config.h> -//#include <Common/Timer.h> - -#include <BBSKernel/MNS/MeqBaseLinPS.h> -#include <BBSKernel/MNS/MeqMatrixTmp.h> -#include <Common/LofarLogger.h> - -using namespace casa; - -namespace LOFAR -{ -namespace BBS -{ - -MeqBaseLinPS::MeqBaseLinPS (const MeqExpr& dft, MeqPointSource* source) -: itsDFT (dft), - itsSource (source) -{ - addChild (itsDFT); - addChild (itsSource->getI()); - addChild (itsSource->getQ()); - addChild (itsSource->getU()); - addChild (itsSource->getV()); -} - -MeqBaseLinPS::~MeqBaseLinPS() -{} - -MeqJonesResult MeqBaseLinPS::getJResult (const MeqRequest& request) -{ - //static NSTimer timer("MeqBaseLinPS::getResult", true); - //timer.start(); - - // Allocate the result. - MeqJonesResult result(request.nspid()); - { - MeqResult& resXX = result.result11(); - MeqResult& resXY = result.result12(); - MeqResult& resYX = result.result21(); - MeqResult& resYY = result.result22(); - // Calculate the source fluxes. - MeqResult ikBuf, qkBuf, ukBuf, vkBuf; - const MeqResult& ik = itsSource->getI().getResultSynced (request, ikBuf); - const MeqResult& qk = itsSource->getQ().getResultSynced (request, qkBuf); - const MeqResult& uk = itsSource->getU().getResultSynced (request, ukBuf); - const MeqResult& vk = itsSource->getV().getResultSynced (request, vkBuf); - // Calculate the baseline DFT. - MeqResult dftBuf; - const MeqResult& dft = itsDFT.getResultSynced (request, dftBuf); - // Calculate the XX values, etc. - MeqMatrix uvk = tocomplex(uk.getValue(), vk.getValue()); - resXX.setValue ((ik.getValue() + qk.getValue()) * dft.getValue()); - resXY.setValue (uvk * dft.getValue()); - resYX.setValue (conj(uvk) * dft.getValue()); - resYY.setValue ((ik.getValue() - qk.getValue()) * dft.getValue()); - // Evaluate (if needed) for the perturbed parameter values. - const MeqParmFunklet* perturbedParm; - for (int spinx=0; spinx<request.nspid(); spinx++) { - bool eval1 = false; - bool eval2 = false; - if (dft.isDefined(spinx)) { - eval1 = true; - eval2 = true; - perturbedParm = dft.getPerturbedParm (spinx); - } else { - if (ik.isDefined(spinx)) { - eval1 = true; - perturbedParm = ik.getPerturbedParm (spinx); - } else if (qk.isDefined(spinx)) { - eval1 = true; - perturbedParm = qk.getPerturbedParm (spinx); - } - if (uk.isDefined(spinx)) { - eval2 = true; - perturbedParm = uk.getPerturbedParm (spinx); - } else if (vk.isDefined(spinx)) { - eval2 = true; - perturbedParm = vk.getPerturbedParm (spinx); - } - } - if (eval1) { - const MeqMatrix& ikp = ik.getPerturbedValue(spinx); - const MeqMatrix& qkp = qk.getPerturbedValue(spinx); - resXX.setPerturbedValue (spinx, - (ikp+qkp) * dft.getPerturbedValue(spinx)); - resYY.setPerturbedValue (spinx, - (ikp-qkp) * dft.getPerturbedValue(spinx)); - resXX.setPerturbedParm (spinx, perturbedParm); - resYY.setPerturbedParm (spinx, perturbedParm); - } - if (eval2) { - MeqMatrix uvk = tocomplex(uk.getPerturbedValue(spinx), - vk.getPerturbedValue(spinx)); - resXY.setPerturbedValue (spinx, - uvk * dft.getPerturbedValue(spinx)); - resYX.setPerturbedValue (spinx, - conj(uvk) * dft.getPerturbedValue(spinx)); - resXY.setPerturbedParm (spinx, perturbedParm); - resYX.setPerturbedParm (spinx, perturbedParm); - } - } - } - //timer.stop(); - return result; -} - -#ifdef EXPR_GRAPH -std::string MeqBaseLinPS::getLabel() -{ - return std::string("MeqBaseLinPS\\nPrediction of a linearly polarized point source\\n" + itsSource->getName() + " (" + itsSource->getGroupName() + ")"); -} -#endif - -} // namespace BBS -} // namespace LOFAR diff --git a/CEP/BB/BBSKernel/src/MNS/MeqGaussianCoherency.cc b/CEP/BB/BBSKernel/src/MNS/MeqGaussianCoherency.cc new file mode 100644 index 00000000000..d68069a3c52 --- /dev/null +++ b/CEP/BB/BBSKernel/src/MNS/MeqGaussianCoherency.cc @@ -0,0 +1,267 @@ +//# MeqGaussianCoherency.cc: Spatial coherence function of an elliptical +//# gaussian source. +//# +//# Copyright (C) 2005 +//# ASTRON (Netherlands Foundation for Research in Astronomy) +//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands, seg@astron.nl +//# +//# This program 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 2 of the License, or +//# (at your option) any later version. +//# +//# This program 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 this program; if not, write to the Free Software +//# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +//# +//# $Id$ + +#include <lofar_config.h> +#include <BBSKernel/MNS/MeqGaussianCoherency.h> +#include <BBSKernel/MNS/MeqMatrix.h> +#include <BBSKernel/MNS/MeqMatrixTmp.h> +#include <Common/lofar_complex.h> +#include <Common/lofar_math.h> +#include <Common/LofarLogger.h> + +using namespace casa; + +namespace LOFAR +{ +namespace BBS +{ +using LOFAR::exp; +using LOFAR::conj; + +MeqGaussianCoherency::MeqGaussianCoherency(const MeqGaussianSource *source, + MeqStatUVW *station1, + MeqStatUVW *station2) + : itsSource(source), + itsStation1(station1), + itsStation2(station2) +{ + addChild(itsSource->getI()); + addChild(itsSource->getQ()); + addChild(itsSource->getU()); + addChild(itsSource->getV()); + addChild(itsSource->getMajor()); + addChild(itsSource->getMinor()); + addChild(itsSource->getPhi()); +} + + +MeqGaussianCoherency::~MeqGaussianCoherency() +{ +} + + +MeqJonesResult MeqGaussianCoherency::getJResult(const MeqRequest &request) +{ + // Evaluate source parameters. + // Note: The result of any parameter is either scalar or it is 2D and + // conforms to the size of the request. + MeqResult ikBuf, qkBuf, ukBuf, vkBuf, majorBuf, minorBuf, phiBuf; + const MeqResult &ik = + getChild(MeqGaussianCoherency::IN_I).getResultSynced(request, ikBuf); + const MeqResult &qk = + getChild(MeqGaussianCoherency::IN_Q).getResultSynced(request, qkBuf); + const MeqResult &uk = + getChild(MeqGaussianCoherency::IN_U).getResultSynced(request, ukBuf); + const MeqResult &vk = + getChild(MeqGaussianCoherency::IN_V).getResultSynced(request, vkBuf); + const MeqResult &major = + getChild(MeqGaussianCoherency::IN_MAJOR).getResultSynced(request, majorBuf); + const MeqResult &minor = + getChild(MeqGaussianCoherency::IN_MINOR).getResultSynced(request, minorBuf); + const MeqResult &phi = + getChild(MeqGaussianCoherency::IN_PHI).getResultSynced(request, phiBuf); + + // Assume major, minor, phi are scalar. + DBGASSERT(!major.getValue().isArray()); + DBGASSERT(!minor.getValue().isArray()); + DBGASSERT(!phi.getValue().isArray()); + + // Note: The result of a MeqStatUVW node is always 1D in time. + const MeqResult &uStation1 = itsStation1->getU(request); + const MeqResult &vStation1 = itsStation1->getV(request); + const MeqResult &uStation2 = itsStation2->getU(request); + const MeqResult &vStation2 = itsStation2->getV(request); + + // Check pre-conditions. + DBGASSERT(uStation1.getValue().nx() == 1 && uStation1.getValue().nelements() == request.ny()); + DBGASSERT(vStation1.getValue().nx() == 1 && vStation1.getValue().nelements() == request.ny()); + DBGASSERT(uStation2.getValue().nx() == 1 && uStation2.getValue().nelements() == request.ny()); + DBGASSERT(vStation2.getValue().nx() == 1 && vStation1.getValue().nelements() == request.ny()); + + // Compute baseline uv-coordinates (1D in time). + MeqMatrix uBaseline = uStation2.getValue() - uStation1.getValue(); + MeqMatrix vBaseline = vStation2.getValue() - vStation1.getValue(); + + // Compute scaled Stokes vector. + MeqMatrix uvScaled = tocomplex(uk.getValue(), vk.getValue()) * 0.5; + MeqMatrix iScaled = ik.getValue() * 0.5; + MeqMatrix qScaled = qk.getValue() * 0.5; + + // Compute spatial coherence for a gaussian source. + MeqMatrix coherence = computeCoherence(request, uBaseline, vBaseline, major.getValue(), minor.getValue(), phi.getValue()); + + // Allocate the result. + MeqJonesResult result(request.nspid()); + MeqResult &resXX = result.result11(); + MeqResult &resXY = result.result12(); + MeqResult &resYX = result.result21(); + MeqResult &resYY = result.result22(); + + resXX.setValue((iScaled + qScaled) * coherence); + resXY.setValue(uvScaled * coherence); + resYX.setValue(conj(uvScaled) * coherence); + resYY.setValue((iScaled - qScaled) * coherence); + + // Evaluate (if needed) for the perturbed parameter values. + const MeqParmFunklet* perturbedParm = 0; + MeqMatrix coherencePerturbed; + MeqMatrix iScaledPerturbed, qScaledPerturbed, uvScaledPerturbed; + + for(int spinx = 0; spinx < request.nspid(); ++spinx) + { + bool evalCoherence = false; + bool evalIQ = false; + bool evalUV = false; + + coherencePerturbed = coherence; + iScaledPerturbed = iScaled; + qScaledPerturbed = qScaled; + uvScaledPerturbed = uvScaled; + + if(major.isDefined(spinx)) + { + evalCoherence = true; + perturbedParm = major.getPerturbedParm (spinx); + } + else if(minor.isDefined(spinx)) + { + evalCoherence = true; + perturbedParm = minor.getPerturbedParm (spinx); + } + else if(phi.isDefined(spinx)) + { + evalCoherence = true; + perturbedParm = phi.getPerturbedParm (spinx); + } + + if(ik.isDefined(spinx)) + { + evalIQ = true; + perturbedParm = ik.getPerturbedParm (spinx); + } + else if(qk.isDefined(spinx)) + { + evalIQ = true; + perturbedParm = qk.getPerturbedParm (spinx); + } + + if(uk.isDefined(spinx)) + { + evalUV = true; + perturbedParm = uk.getPerturbedParm (spinx); + } + else if(vk.isDefined(spinx)) + { + evalUV = true; + perturbedParm = vk.getPerturbedParm (spinx); + } + + if(evalCoherence) + { + coherencePerturbed = computeCoherence(request, uBaseline, vBaseline, + major.getPerturbedValue(spinx), + minor.getPerturbedValue(spinx), + phi.getPerturbedValue(spinx)); + } + + if(evalIQ) + { + iScaledPerturbed = ik.getPerturbedValue(spinx) * 0.5; + qScaledPerturbed = qk.getPerturbedValue(spinx) * 0.5; + } + + if(evalUV) + { + uvScaledPerturbed = tocomplex(uk.getPerturbedValue(spinx), + vk.getPerturbedValue(spinx)) * 0.5; + } + + if(evalCoherence || evalIQ) + { + resXX.setPerturbedParm(spinx, perturbedParm); + resXX.setPerturbedValue(spinx, + (iScaledPerturbed + qScaledPerturbed) * coherencePerturbed); + + resYY.setPerturbedParm(spinx, perturbedParm); + resYY.setPerturbedValue(spinx, + (iScaledPerturbed - qScaledPerturbed) * coherencePerturbed); + } + + if(evalCoherence || evalUV) + { + resXY.setPerturbedParm(spinx, perturbedParm); + resXY.setPerturbedValue(spinx, uvScaledPerturbed * coherencePerturbed); + + resYX.setPerturbedParm(spinx, perturbedParm); + resYX.setPerturbedValue(spinx, conj(uvScaledPerturbed) * coherencePerturbed); + } + } + + return result; +} + + +MeqMatrix MeqGaussianCoherency::computeCoherence(const MeqRequest &request, + const MeqMatrix &uBaseline, const MeqMatrix &vBaseline, + const MeqMatrix &major, const MeqMatrix &minor, const MeqMatrix &phi) +{ + // Compute dot product of rotated, scaled uv-vector with itself (1D in time). + MeqMatrix cosPhi(cos(phi)); + MeqMatrix sinPhi(sin(phi)); + MeqMatrix uvTransformedDotProduct = + sqr(major * (uBaseline * cosPhi - vBaseline * sinPhi)) + + sqr(minor * (uBaseline * sinPhi + vBaseline * cosPhi)); + + // Allocate the result matrix (2D). + MeqMatrix result; + double *valuep = result.setDoubleFormat(request.nx(), request.ny()); + + // Compute unnormalized spatial coherence (2D). + double cSqr = casa::C::c * casa::C::c; + for(int t = 0; t < request.ny(); ++t) + { + double freq = request.domain().startX() + request.stepX() / 2.0; + const double uv = casa::C::_2pi * uvTransformedDotProduct.getDouble(0, t); + + for(int f = 0; f < request.nx(); ++f) + { + *(valuep++) = exp(-(freq * freq * uv) / cSqr); + freq += request.stepX(); + } + } + + return result; +} + + +#ifdef EXPR_GRAPH +std::string MeqGaussianCoherency::getLabel() +{ + return std::string("MeqGaussianCoherency\\nSpatial coherence function of" + " an elliptical gaussian source\\n" + itsSource->getName() + " (" + + itsSource->getGroupName() + ")"); +} +#endif + +} // namespace BBS +} // namespace LOFAR diff --git a/CEP/BB/BBSKernel/src/MNS/MeqGaussSource.cc b/CEP/BB/BBSKernel/src/MNS/MeqGaussianSource.cc similarity index 59% rename from CEP/BB/BBSKernel/src/MNS/MeqGaussSource.cc rename to CEP/BB/BBSKernel/src/MNS/MeqGaussianSource.cc index 9543710f2ac..55d9b3ac3d9 100644 --- a/CEP/BB/BBSKernel/src/MNS/MeqGaussSource.cc +++ b/CEP/BB/BBSKernel/src/MNS/MeqGaussianSource.cc @@ -1,4 +1,4 @@ -//# MeqGaussSource.cc: Class holding the expressions defining a gauss source +//# MeqGaussianSource.cc: Class holding the expressions defining a gauss source //# //# Copyright (C) 2002 //# ASTRON (Netherlands Foundation for Research in Astronomy) @@ -21,33 +21,35 @@ //# $Id$ #include <lofar_config.h> -#include <BBSKernel/MNS/MeqGaussSource.h> - +#include <BBSKernel/MNS/MeqGaussianSource.h> namespace LOFAR { namespace BBS { -MeqGaussSource::MeqGaussSource (const string& name, - const string& groupName, - const MeqExpr& fluxI, const MeqExpr& fluxQ, - const MeqExpr& fluxU, const MeqExpr& fluxV, - const MeqExpr& ra, const MeqExpr& dec, - const MeqExpr& minor, const MeqExpr& major, - const MeqExpr& phi) -: MeqSource (name, groupName, ra, dec), - itsI (fluxI), - itsQ (fluxQ), - itsU (fluxU), - itsV (fluxV), - itsMinor (minor), - itsMajor (major), - itsPhi (phi) -{} +MeqGaussianSource::MeqGaussianSource(const string &name, + const string& groupName, + const MeqExpr& fluxI, const MeqExpr& fluxQ, + const MeqExpr& fluxU, const MeqExpr& fluxV, + const MeqExpr& ra, const MeqExpr& dec, + const MeqExpr& major, const MeqExpr& minor, + const MeqExpr& phi) + : MeqSource(name, groupName, ra, dec), + itsI(fluxI), + itsQ(fluxQ), + itsU(fluxU), + itsV(fluxV), + itsMajor(major), + itsMinor(minor), + itsPhi(phi) +{ +} + -MeqGaussSource::~MeqGaussSource() -{} +MeqGaussianSource::~MeqGaussianSource() +{ +} } // namespace BBS } // namespace LOFAR diff --git a/CEP/BB/BBSKernel/src/MNS/MeqJonesMul.cc b/CEP/BB/BBSKernel/src/MNS/MeqJonesMul.cc new file mode 100644 index 00000000000..2f9c5537d43 --- /dev/null +++ b/CEP/BB/BBSKernel/src/MNS/MeqJonesMul.cc @@ -0,0 +1,159 @@ +//# MeqJonesMul.h: Multiply each component of a MeqJonesExpr with a single +//# MeqExpr. +//# +//# Copyright (C) 2007 +//# ASTRON (Netherlands Foundation for Research in Astronomy) +//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands, seg@astron.nl +//# +//# This program 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 2 of the License, or +//# (at your option) any later version. +//# +//# This program 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 this program; if not, write to the Free Software +//# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +//# +//# $Id$ + +#include <lofar_config.h> + +#include <BBSKernel/MNS/MeqJonesMul.h> +#include <BBSKernel/MNS/MeqRequest.h> +#include <BBSKernel/MNS/MeqMatrix.h> +#include <BBSKernel/MNS/MeqMatrixTmp.h> + +namespace LOFAR +{ +namespace BBS +{ + +MeqJonesMul::MeqJonesMul(const MeqJonesExpr &left, const MeqExpr &right) + : itsLeft(left), + itsRight(right) +{ + addChild(left); + addChild(right); +} + + +MeqJonesMul::~MeqJonesMul() +{ +} + + +MeqJonesResult MeqJonesMul::getJResult(const MeqRequest &request) +{ + // Create the result object. + MeqJonesResult result(request.nspid()); + MeqResult& result11 = result.result11(); + MeqResult& result12 = result.result12(); + MeqResult& result21 = result.result21(); + MeqResult& result22 = result.result22(); + + // Evaluate the children. + MeqJonesResult leftResult; + MeqResult rightResult; + + const MeqJonesResult &left = itsLeft.getResultSynced(request, leftResult); + const MeqResult &right = itsRight.getResultSynced(request, rightResult); + + const MeqResult &l11 = left.getResult11(); + const MeqResult &l12 = left.getResult12(); + const MeqResult &l21 = left.getResult21(); + const MeqResult &l22 = left.getResult22(); + + const MeqMatrix &ml11 = l11.getValue(); + const MeqMatrix &ml12 = l12.getValue(); + const MeqMatrix &ml21 = l21.getValue(); + const MeqMatrix &ml22 = l22.getValue(); + + const MeqMatrix &mr = right.getValue(); + + // Compute result. + result11.setValue(ml11 * mr); + result12.setValue(ml12 * mr); + result21.setValue(ml21 * mr); + result22.setValue(ml22 * mr); + + // Compute perturbed values. + for (int spinx=0; spinx<request.nspid(); spinx++) + { + if(right.isDefined(spinx)) + { + const MeqMatrix &mr = right.getPerturbedValue(spinx); + const MeqMatrix &ml11 = l11.getPerturbedValue(spinx); + const MeqMatrix &ml12 = l12.getPerturbedValue(spinx); + const MeqMatrix &ml21 = l21.getPerturbedValue(spinx); + const MeqMatrix &ml22 = l22.getPerturbedValue(spinx); + + const MeqParmFunklet *perturbedParm = right.getPerturbedParm(spinx); + result11.setPerturbedParm(spinx, perturbedParm); + result11.setPerturbedValue(spinx, ml11 * mr); + result12.setPerturbedParm(spinx, perturbedParm); + result12.setPerturbedValue(spinx, ml12 * mr); + result21.setPerturbedParm(spinx, perturbedParm); + result21.setPerturbedValue(spinx, ml21 * mr); + result22.setPerturbedParm(spinx, perturbedParm); + result22.setPerturbedValue(spinx, ml22 * mr); + } + else + { + if(l11.isDefined(spinx)) + { + const MeqMatrix &mr = right.getPerturbedValue(spinx); + const MeqMatrix &ml11 = l11.getPerturbedValue(spinx); + + result11.setPerturbedParm(spinx, l11.getPerturbedParm(spinx)); + result11.setPerturbedValue(spinx, ml11 * mr); + } + + if(l12.isDefined(spinx)) + { + const MeqMatrix &mr = right.getPerturbedValue(spinx); + const MeqMatrix &ml12 = l12.getPerturbedValue(spinx); + + result12.setPerturbedParm(spinx, l12.getPerturbedParm(spinx)); + result12.setPerturbedValue(spinx, ml12 * mr); + } + + if(l21.isDefined(spinx)) + { + const MeqMatrix &mr = right.getPerturbedValue(spinx); + const MeqMatrix &ml21 = l21.getPerturbedValue(spinx); + + result21.setPerturbedParm(spinx, l21.getPerturbedParm(spinx)); + result21.setPerturbedValue(spinx, ml21 * mr); + } + + if(l22.isDefined(spinx)) + { + const MeqMatrix &mr = right.getPerturbedValue(spinx); + const MeqMatrix &ml22 = l22.getPerturbedValue(spinx); + + result22.setPerturbedParm(spinx, l22.getPerturbedParm(spinx)); + result22.setPerturbedValue(spinx, ml22 * mr); + } + } + } + + return result; +} + + +#ifdef EXPR_GRAPH +virtual std::string MeqJonesMul::getLabel() +{ + return std::string("MeqJonesMul\\n" + "Multiply MeqJonesExpr with a single MeqExpr"); +} +#endif + + +} // namespace BBS +} // namespace LOFAR diff --git a/CEP/BB/BBSKernel/src/MNS/MeqPhaseShift.cc b/CEP/BB/BBSKernel/src/MNS/MeqPhaseShift.cc new file mode 100644 index 00000000000..f8349986122 --- /dev/null +++ b/CEP/BB/BBSKernel/src/MNS/MeqPhaseShift.cc @@ -0,0 +1,219 @@ +//# MeqPhaseShift.cc: Phase delay due to baseline geometry with respect to +//# source direction. +//# +//# Copyright (C) 2005 +//# ASTRON (Netherlands Foundation for Research in Astronomy) +//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands, seg@astron.nl +//# +//# This program 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 2 of the License, or +//# (at your option) any later version. +//# +//# This program 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 this program; if not, write to the Free Software +//# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +//# +//# $Id$ + +#include <lofar_config.h> + +#include <BBSKernel/MNS/MeqPhaseShift.h> +#include <BBSKernel/MNS/MeqDFTPS.h> +#include <BBSKernel/MNS/MeqStatUVW.h> +#include <BBSKernel/MNS/MeqMatrixTmp.h> +#include <Common/LofarLogger.h> +#include <Common/lofar_iomanip.h> + +using namespace casa; + +namespace LOFAR +{ +namespace BBS +{ +using LOFAR::dcomplex; +using LOFAR::conj; + + +MeqPhaseShift::MeqPhaseShift (const MeqExpr& left, const MeqExpr& right) + : itsLeft (left), + itsRight (right) +{ + addChild (itsLeft); + addChild (itsRight); +} + + +MeqPhaseShift::~MeqPhaseShift() +{ +} + + +MeqResult MeqPhaseShift::getResult (const MeqRequest& request) +{ + MeqResult result(request.nspid()); + + // Get N (for the division). + // Assert it is a scalar value. +// MeqResultVec lmnRes; +// const MeqResultVec& lmn = itsLMN.getResultVecSynced (request, lmnRes); +// const MeqResult& nk = lmn[2]; +// DBGASSERT (nk.getValue().nelements() == request.ny() || +// nk.getValue().nelements() == 1); + + // Calculate the left and right station Jones matrix elements. + // A delta in the station source predict is only available if multiple + // frequency channels are used. + bool multFreq = request.nx() > 1; + MeqResultVec reslBuf, resrBuf; + const MeqResultVec& resl = itsLeft.getResultVecSynced (request, resrBuf); + const MeqResultVec& resr = itsRight.getResultVecSynced (request, reslBuf); + const MeqResult& left = resl[0]; + const MeqResult& right = resr[0]; + const MeqResult& leftDelta = resl[1]; + const MeqResult& rightDelta = resr[1]; + + // It is tried to compute the DFT as efficient as possible. + // Therefore the baseline contribution is split into its antenna parts. + // dft = exp(2i.pi(ul+vm+wn)) / n (with u,v,w in wavelengths) + // = (exp(2i.pi((u1.l+v1.m+w1.n) - (u2.l+v2.m+w2.n))/wvl))/n (u,v,w in m) + // = ((exp(i(u1.l+v1.m+w1.n)) / exp(i(u2.l+v2.m+w2.m))) ^ (2.pi/wvl)) / n + // So left and right return the exp values independent of wavelength. + // Thereafter they are scaled to the freq domain by raising the values + // for each time to the appropriate powers. + // Alas the rule + // x^(a*b) = (x^a)^b + // which is valid for real numbers, is only valid for complex numbers + // if b is an integer number. + // Therefore the station calculations (in MeqStatSources) are done as + // follows, where it should be noted that the frequencies are regularly + // strided. + // f = f0 + k.df (k = 0 ... nchan-1) + // s1 = (u1.l+v1.m+w1.n).2i.pi/c + // s2 = (u2.l+v2.m+w2.n).2i.pi/c + // dft = exp(s1(f0+k.df)) / exp(s2(f0+k.df)) / n + // = (exp(s1.f0)/exp(s2.f0)) . (exp(s1.k.df)/exp(s2.k.df)) / n + // = (exp(s1.f0)/exp(s2.f0)) . (exp(s1.df)/exp(s2.df))^k / n + // In principle the power is expensive, but because the frequencies are + // regularly strided, it is possible to use multiplication. + // So it gets + // dft(f0) = (exp(s1.f0)/exp(s2.f0)) / n + // dft(fj) = dft(fj-1) * (exp(s1.df)/exp(s2.df)) + // Using a python script (tuvw.py) is is checked that this way of + // calculation is accurate enough. + // Another optimization can be achieved in the division of the two + // complex numbers which can be turned into a cheaper multiplication. + // exp(x)/exp(y) = (cos(x) + i.sin(x)) / (cos(y) + i.sin(y)) + // = (cos(x) + i.sin(x)) * (cos(y) - i.sin(y)) + + /* + { + + MeqDFTPS *dftpsl = dynamic_cast<MeqDFTPS*>(itsLeft.rep()); + MeqDFTPS *dftpsr = dynamic_cast<MeqDFTPS*>(itsRight.rep()); + + const MeqResult& ul = dftpsl->uvw()->getU(request); + const MeqResult& vl = dftpsl->uvw()->getV(request); + const MeqResult& wl = dftpsl->uvw()->getW(request); + const MeqResult& ur = dftpsr->uvw()->getU(request); + const MeqResult& vr = dftpsr->uvw()->getV(request); + const MeqResult& wr = dftpsr->uvw()->getW(request); + cout << "U: " << setprecision(20) << ur.getValue() - ul.getValue() + << endl; + cout << "V: " << setprecision(20) << vr.getValue() - vl.getValue() + << endl; + cout << "W: " << setprecision(20) << wr.getValue() - wl.getValue() + << endl; + LOG_TRACE_FLOW ("U: " << ur.getValue() - ul.getValue()); + LOG_TRACE_FLOW ("V: " << vr.getValue() - vl.getValue()); + LOG_TRACE_FLOW ("W: " << wr.getValue() - wl.getValue()); + } + */ + + MeqMatrix res(makedcomplex(0,0), request.nx(), request.ny(), false); + for(int iy=0; iy<request.ny(); ++iy) + { + dcomplex tmpl = left.getValue().getDComplex(0,iy); + dcomplex tmpr = right.getValue().getDComplex(0,iy); + + // We have to divide by N. + // However, we divide by 2N to get the factor 0.5 needed in (I+Q)/2, etc. + // in MeqBaseLinPS. + // double tmpnk = 2. * nk.getValue().getDouble(0,iy); +// double tmpnk = 2.0; + dcomplex factor; + if(multFreq) + { + dcomplex deltal = leftDelta.getValue().getDComplex(0,iy); + dcomplex deltar = rightDelta.getValue().getDComplex(0,iy); + factor = deltar * conj(deltal); + } +// res.fillRowWithProducts (tmpr * conj(tmpl) / tmpnk, factor, iy); + res.fillRowWithProducts (tmpr * conj(tmpl), factor, iy); + } + result.setValue (res); + + // cout << "DFT:" << endl; + // cout << setprecision(20) << res << endl; + + // Evaluate (if needed) for the perturbed parameter values. + // Note that we do not have to test for perturbed values in nk, + // because the left and right value already depend on nk. + const MeqParmFunklet* perturbedParm; + + for(int spinx=0; spinx<request.nspid(); spinx++) + { + bool eval = false; + if(left.isDefined(spinx)) + { + eval = true; + perturbedParm = left.getPerturbedParm(spinx); + } + else if(right.isDefined(spinx)) + { + eval = true; + perturbedParm = right.getPerturbedParm(spinx); + } + + if(eval) + { + MeqMatrix pres(makedcomplex(0,0), request.nx(), request.ny(), false); + for(int iy=0; iy<request.ny(); ++iy) + { + dcomplex tmpl = left.getPerturbedValue(spinx).getDComplex(0,iy); + dcomplex tmpr = right.getPerturbedValue(spinx).getDComplex(0,iy); + + // double tmpnk = 2. * nk.getPerturbedValue(spinx).getDouble(0,iy); +// double tmpnk = 2.0; + dcomplex factor; + if(multFreq) + { + dcomplex deltal = leftDelta.getPerturbedValue(spinx).getDComplex(0,iy); + dcomplex deltar = rightDelta.getPerturbedValue(spinx).getDComplex(0,iy); + factor = deltar * conj(deltal); + } +// pres.fillRowWithProducts(tmpr * conj(tmpl) / tmpnk, factor, iy); + pres.fillRowWithProducts(tmpr * conj(tmpl), factor, iy); + result.setPerturbedValue (spinx, pres); + result.setPerturbedParm (spinx, perturbedParm); + } + } + } + + return result; +} + +#ifdef EXPR_GRAPH +std::string MeqPhaseShift::getLabel() +{ + return std::string("MeqPhaseShift\\nPhase delay due to baseline geometry."); +} +#endif + +} // namespace BBS +} // namespace LOFAR diff --git a/CEP/BB/BBSKernel/src/MNS/MeqPointCoherency.cc b/CEP/BB/BBSKernel/src/MNS/MeqPointCoherency.cc new file mode 100644 index 00000000000..30e8b82c0d9 --- /dev/null +++ b/CEP/BB/BBSKernel/src/MNS/MeqPointCoherency.cc @@ -0,0 +1,149 @@ +//# MeqPointCoherency.h: Spatial coherence function of a point source. +//# +//# Copyright (C) 2005 +//# ASTRON (Netherlands Foundation for Research in Astronomy) +//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands, seg@astron.nl +//# +//# This program 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 2 of the License, or +//# (at your option) any later version. +//# +//# This program 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 this program; if not, write to the Free Software +//# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +//# +//# $Id$ + +#include <lofar_config.h> +//#include <Common/Timer.h> + +#include <BBSKernel/MNS/MeqPointCoherency.h> +#include <BBSKernel/MNS/MeqMatrixTmp.h> +#include <Common/LofarLogger.h> + +using namespace casa; + +namespace LOFAR +{ +namespace BBS +{ + +MeqPointCoherency::MeqPointCoherency(const MeqPointSource *source) + : itsSource(source) +{ + addChild (itsSource->getI()); + addChild (itsSource->getQ()); + addChild (itsSource->getU()); + addChild (itsSource->getV()); +} + + +MeqPointCoherency::~MeqPointCoherency() +{ +} + + +MeqJonesResult MeqPointCoherency::getJResult(const MeqRequest &request) +{ + //static NSTimer timer("MeqPointCoherency::getResult", true); + //timer.start(); + + // Allocate the result. + MeqJonesResult result(request.nspid()); + MeqResult& resXX = result.result11(); + MeqResult& resXY = result.result12(); + MeqResult& resYX = result.result21(); + MeqResult& resYY = result.result22(); + + // Calculate the source fluxes. + MeqResult ikBuf, qkBuf, ukBuf, vkBuf; + const MeqResult& ik = itsSource->getI().getResultSynced (request, ikBuf); + const MeqResult& qk = itsSource->getQ().getResultSynced (request, qkBuf); + const MeqResult& uk = itsSource->getU().getResultSynced (request, ukBuf); + const MeqResult& vk = itsSource->getV().getResultSynced (request, vkBuf); + + // Calculate the XX values, etc. + MeqMatrix uvk_2 = tocomplex(uk.getValue(), vk.getValue()) * 0.5; + MeqMatrix ik_2 = ik.getValue() * 0.5; + MeqMatrix qk_2 = qk.getValue() * 0.5; + + resXX.setValue (ik_2 + qk_2); + resXY.setValue (uvk_2); + resYX.setValue (conj(uvk_2)); + resYY.setValue (ik_2 - qk_2); + + // Evaluate (if needed) for the perturbed parameter values. + const MeqParmFunklet* perturbedParm; + for(int spinx=0; spinx<request.nspid(); spinx++) + { + bool eval1 = false; + bool eval2 = false; + + if(ik.isDefined(spinx)) + { + eval1 = true; + perturbedParm = ik.getPerturbedParm (spinx); + } + else if(qk.isDefined(spinx)) + { + eval1 = true; + perturbedParm = qk.getPerturbedParm (spinx); + } + + if(uk.isDefined(spinx)) + { + eval2 = true; + perturbedParm = uk.getPerturbedParm (spinx); + } + else if(vk.isDefined(spinx)) + { + eval2 = true; + perturbedParm = vk.getPerturbedParm (spinx); + } + + if(eval1) + { + const MeqMatrix& ikp_2 = ik.getPerturbedValue(spinx) * 0.5; + const MeqMatrix& qkp_2 = qk.getPerturbedValue(spinx) * 0.5; + + resXX.setPerturbedParm (spinx, perturbedParm); + resXX.setPerturbedValue(spinx, (ikp_2 + qkp_2)); + + resYY.setPerturbedParm (spinx, perturbedParm); + resYY.setPerturbedValue(spinx, (ikp_2 - qkp_2)); + } + + if(eval2) + { + MeqMatrix uvkp_2 = tocomplex(uk.getPerturbedValue(spinx), + vk.getPerturbedValue(spinx)) * 0.5; + + resXY.setPerturbedParm(spinx, perturbedParm); + resXY.setPerturbedValue(spinx, uvkp_2); + + resYX.setPerturbedParm (spinx, perturbedParm); + resYX.setPerturbedValue(spinx, conj(uvkp_2)); + } + } + + //timer.stop(); + return result; +} + +#ifdef EXPR_GRAPH +std::string MeqPointCoherency::getLabel() +{ + return std::string("MeqPointCoherency\\nSpatial coherence function of a" + " point source\\n" + itsSource->getName() + " (" + + itsSource->getGroupName() + ")"); +} +#endif + +} // namespace BBS +} // namespace LOFAR diff --git a/CEP/BB/BBSKernel/src/MNS/MeqSourceList.cc b/CEP/BB/BBSKernel/src/MNS/MeqSourceList.cc index 2702765da80..a61d8714ec3 100644 --- a/CEP/BB/BBSKernel/src/MNS/MeqSourceList.cc +++ b/CEP/BB/BBSKernel/src/MNS/MeqSourceList.cc @@ -23,7 +23,7 @@ #include <lofar_config.h> #include <BBSKernel/MNS/MeqSourceList.h> #include <BBSKernel/MNS/MeqPointSource.h> -#include <BBSKernel/MNS/MeqGaussSource.h> +#include <BBSKernel/MNS/MeqGaussianSource.h> #include <BBSKernel/MNS/MeqParmFunklet.h> #include <Common/LofarLogger.h> @@ -106,14 +106,16 @@ MeqSourceList::MeqSourceList (LOFAR::ParmDB::ParmDB& parmTable, MeqParmGroup& gr if (std::find(gnams.begin(), gnams.end(), name) == gnams.end()) { add (new MeqPointSource(name, name, mi, mq, mu, mv, mr, md)); } else { + LOG_DEBUG_STR("Found gaussian source: " << name); + MeqExpr mmin = MeqParmFunklet::create ("Minor:"+name, group, &parmTable); MeqExpr mmaj = MeqParmFunklet::create ("Major:"+name, group, &parmTable); MeqExpr mphi = MeqParmFunklet::create ("Phi:"+name, group, &parmTable); - add (new MeqGaussSource(name, name, mi, mq, mu, mv, mr, md, - mmin, mmaj, mphi)); + add (new MeqGaussianSource(name, name, mi, mq, mu, mv, mr, md, + mmaj, mmin, mphi)); } // cout << "Found source " << name << " (srcnr=" << srcnr << ')' << endl; } diff --git a/CEP/BB/BBSKernel/src/Makefile.am b/CEP/BB/BBSKernel/src/Makefile.am index 90099c3602f..591dd78abb4 100644 --- a/CEP/BB/BBSKernel/src/Makefile.am +++ b/CEP/BB/BBSKernel/src/Makefile.am @@ -9,19 +9,19 @@ libbbskernel_la_SOURCES = Package__Version.cc \ VisData.cc \ VisSelection.cc \ MNS/MeqAzEl.cc \ - MNS/MeqBaseDFTPS.cc \ - MNS/MeqBaseLinPS.cc \ MNS/MeqDFTPS.cc \ MNS/MeqDiag.cc \ MNS/MeqDipoleBeam.cc \ MNS/MeqDomain.cc \ MNS/MeqExpr.cc \ MNS/MeqFunklet.cc \ - MNS/MeqGaussSource.cc \ + MNS/MeqGaussianCoherency.cc \ + MNS/MeqGaussianSource.cc \ MNS/MeqJonesCMul2.cc \ MNS/MeqJonesCMul3.cc \ MNS/MeqJonesExpr.cc \ MNS/MeqJonesInvert.cc \ + MNS/MeqJonesMul.cc \ MNS/MeqJonesMul2.cc \ MNS/MeqJonesMul3.cc \ MNS/MeqJonesNode.cc \ @@ -40,6 +40,8 @@ libbbskernel_la_SOURCES = Package__Version.cc \ MNS/MeqParmFunklet.cc \ MNS/MeqParmSingle.cc \ MNS/MeqPhaseRef.cc \ + MNS/MeqPhaseShift.cc \ + MNS/MeqPointCoherency.cc \ MNS/MeqPointSource.cc \ MNS/MeqPolc.cc \ MNS/MeqPolcLog.cc \ diff --git a/CEP/BB/BBSKernel/src/Model.cc b/CEP/BB/BBSKernel/src/Model.cc index b7ef623e280..719057951cc 100644 --- a/CEP/BB/BBSKernel/src/Model.cc +++ b/CEP/BB/BBSKernel/src/Model.cc @@ -30,8 +30,10 @@ #include <BBSKernel/MNS/MeqParmSingle.h> #include <BBSKernel/MNS/MeqDiag.h> #include <BBSKernel/MNS/MeqJonesInvert.h> -#include <BBSKernel/MNS/MeqBaseDFTPS.h> -#include <BBSKernel/MNS/MeqBaseLinPS.h> +#include <BBSKernel/MNS/MeqPhaseShift.h> +#include <BBSKernel/MNS/MeqPointCoherency.h> +#include <BBSKernel/MNS/MeqGaussianCoherency.h> +#include <BBSKernel/MNS/MeqJonesMul.h> #include <BBSKernel/MNS/MeqJonesCMul3.h> #include <BBSKernel/MNS/MeqJonesSum.h> #include <BBSKernel/MNS/MeqJonesVisData.h> @@ -362,14 +364,28 @@ void Model::makeEquations(EquationType type, const vector<string> &components, { // Phase shift (incorporates geometry and fringe stopping). MeqExpr shift - (new MeqBaseDFTPS(dft[baseline.first * nSources + j], - dft[baseline.second * nSources + j], itsLMNNodes[j])); + (new MeqPhaseShift(dft[baseline.first * nSources + j], + dft[baseline.second * nSources + j])); + + MeqJonesExpr sourceExpr; - // Point source. MeqPointSource *source = dynamic_cast<MeqPointSource*>(itsSourceNodes[j]); - MeqJonesExpr sourceExpr = new MeqBaseLinPS(shift, source); + if(source) + { + // Point source. + sourceExpr = new MeqPointCoherency(source); + } + else + { + MeqGaussianSource *source = dynamic_cast<MeqGaussianSource*>(itsSourceNodes[j]); + ASSERT(source); + sourceExpr = new MeqGaussianCoherency(source, itsUVWNodes[baseline.first].get(), itsUVWNodes[baseline.second].get()); + } + + // Phase shift the source coherency. + sourceExpr = new MeqJonesMul(sourceExpr, shift); // Apply image plane effects if required. if(mask[DIPOLE_BEAM]) -- GitLab