diff --git a/CEP/BB/BBSKernel/src/MNS/MeqAzEl.cc b/CEP/BB/BBSKernel/src/MNS/MeqAzEl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e6359c41398531def5465294e3fadc4c6b667244
--- /dev/null
+++ b/CEP/BB/BBSKernel/src/MNS/MeqAzEl.cc
@@ -0,0 +1,180 @@
+//# MeqAzEl.cc: Azimuth and elevation for a direction (ra,dec) on the sky.
+//#
+//# 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/MeqAzEl.h>
+#include <BBSKernel/MNS/MeqSource.h>
+#include <BBSKernel/MNS/MeqStation.h>
+#include <BBSKernel/MNS/MeqRequest.h>
+#include <Common/LofarLogger.h>
+
+#include <measures/Measures/MDirection.h>
+#include <measures/Measures/MEpoch.h>
+#include <measures/Measures/MPosition.h>
+#include <measures/Measures/MeasFrame.h>
+#include <measures/Measures/MeasConvert.h>
+#include <casa/Quanta/Quantum.h>
+
+using namespace casa;
+
+namespace LOFAR
+{
+namespace BBS
+{
+
+MeqAzEl::MeqAzEl(MeqSource &source, MeqStation &station)
+{
+ addChild(source.getRa());
+ addChild(source.getDec());
+ addChild(station.getPosX());
+ addChild(station.getPosY());
+ addChild(station.getPosZ());
+}
+
+
+MeqResultVec MeqAzEl::getResultVec(const MeqRequest& request)
+{
+ // Check preconditions.
+ ASSERTSTR(request.ny() > 0, "Need time values.");
+
+ // Evaluate children.
+ MeqResult res_ra, res_dec, res_x, res_y, res_z;
+ const MeqResult &ra =
+ getChild(MeqAzEl::IN_RA).getResultSynced(request, res_ra);
+ const MeqResult &dec =
+ getChild(MeqAzEl::IN_DEC).getResultSynced(request, res_dec);
+ const MeqResult &x =
+ getChild(MeqAzEl::IN_X).getResultSynced(request, res_x);
+ const MeqResult &y =
+ getChild(MeqAzEl::IN_Y).getResultSynced(request, res_y);
+ const MeqResult &z =
+ getChild(MeqAzEl::IN_Z).getResultSynced(request, res_z);
+
+ // Check preconditions.
+ ASSERTSTR(!ra.getValue().isArray() && !dec.getValue().isArray(), "Variable"
+ " source positions are not supported yet.");
+ ASSERTSTR(!x.getValue().isArray() && !y.getValue().isArray()
+ && !z.getValue().isArray(), "Variable station positions are not"
+ " supported yet.");
+
+ // Create result.
+ MeqResultVec result(2, request.nspid());
+
+ // Evaluate main value.
+ evaluate(request, ra.getValue(), dec.getValue(), x.getValue(), y.getValue(),
+ z.getValue(), result[0].getValueRW(), result[1].getValueRW());
+
+ // Evaluate perturbed values.
+ const MeqParmFunklet *perturbedParm;
+ for(int i = 0; i < request.nspid(); ++i)
+ {
+ // Find out if this perturbed value needs to be computed.
+ if(ra.isDefined(i))
+ {
+ perturbedParm = ra.getPerturbedParm(i);
+ }
+ else if(dec.isDefined(i))
+ {
+ perturbedParm = dec.getPerturbedParm(i);
+ }
+ else if(x.isDefined(i))
+ {
+ perturbedParm = x.getPerturbedParm(i);
+ }
+ else if(y.isDefined(i))
+ {
+ perturbedParm = y.getPerturbedParm(i);
+ }
+ else if(z.isDefined(i))
+ {
+ perturbedParm = z.getPerturbedParm(i);
+ }
+ else
+ {
+ continue;
+ }
+
+ evaluate(request, ra.getPerturbedValue(i), dec.getPerturbedValue(i),
+ x.getPerturbedValue(i), y.getPerturbedValue(i),
+ z.getPerturbedValue(i), result[0].getPerturbedValueRW(i),
+ result[1].getPerturbedValueRW(i));
+
+ result[0].setPerturbedParm(i, perturbedParm);
+ result[1].setPerturbedParm(i, perturbedParm);
+ }
+
+ return result;
+}
+
+
+void MeqAzEl::evaluate(const MeqRequest& request, const MeqMatrix &in_ra,
+ const MeqMatrix &in_dec, const MeqMatrix &in_x, const MeqMatrix &in_y,
+ const MeqMatrix &in_z, MeqMatrix &out_az, MeqMatrix &out_el)
+{
+ MPosition position(MVPosition(in_x.getDouble(0, 0), in_y.getDouble(0, 0),
+ in_z.getDouble(0, 0)), MPosition::Ref(MPosition::ITRF));
+ Quantum<double> qepoch(0, "s");
+ MEpoch epoch(qepoch, MEpoch::UTC);
+
+ // Create and initialize a frame.
+ MeasFrame frame;
+ frame.set(position);
+ frame.set(epoch);
+
+ // Create conversion engine.
+ MDirection dir(MVDirection(in_ra.getDouble(0, 0), in_dec.getDouble(0, 0)),
+ MDirection::Ref(MDirection::J2000));
+ MDirection::Convert converter = MDirection::Convert(dir,
+ MDirection::Ref(MDirection::AZEL, frame));
+
+ // Result is only time variable.
+ double *az = out_az.setDoubleFormat(1, request.ny());
+ double *el = out_el.setDoubleFormat(1, request.ny());
+
+ for(int i = 0; i < request.ny(); ++i)
+ {
+ // Update reference frame.
+ qepoch.setValue(request.y(i));
+ epoch.set(qepoch);
+ frame.set(epoch);
+
+ // Compute azimuth and elevation.
+ MDirection azel(converter());
+ Vector<Double> vec_azel = azel.getValue().getAngle("rad").getValue();
+ *az = vec_azel(0);
+ *el = vec_azel(1);
+ ++az; ++el;
+ }
+}
+
+
+#ifdef EXPR_GRAPH
+std::string MeqAzEl::getLabel()
+{
+ return std::string("MeqAzEl\\nAzimuth and elevation of a source.");
+}
+#endif
+
+} //# namespace BBS
+} //# namespace LOFAR
+
diff --git a/CEP/BB/BBSKernel/src/MNS/MeqDipoleBeam.cc b/CEP/BB/BBSKernel/src/MNS/MeqDipoleBeam.cc
new file mode 100644
index 0000000000000000000000000000000000000000..b1b92fb77d3a4c788668be1ad2624ee9346bbaa4
--- /dev/null
+++ b/CEP/BB/BBSKernel/src/MNS/MeqDipoleBeam.cc
@@ -0,0 +1,228 @@
+//# MeqDipoleBeam.cc: Dipole voltage beam (analytic)
+//#
+//# 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/MeqDipoleBeam.h>
+#include <Common/lofar_complex.h>
+
+namespace LOFAR
+{
+namespace BBS
+{
+//using LOFAR::dcomplex;
+//using LOFAR::conj;
+
+
+MeqDipoleBeam::MeqDipoleBeam(MeqExpr azel, double height, double length,
+ double slant, double orientation)
+ : itsHeight(height),
+ itsLength(length),
+ itsSlant(slant),
+ itsOrientation(orientation)
+{
+ addChild(azel);
+}
+
+
+MeqJonesResult MeqDipoleBeam::getJResult(const MeqRequest &request)
+{
+ // Evaluate children.
+ MeqResultVec res_azel;
+ const MeqResultVec &azel =
+ getChild(MeqDipoleBeam::IN_AZEL).getResultVecSynced(request, res_azel);
+
+ // Create result.
+ MeqJonesResult result(request.nspid());
+ MeqResult& result11 = result.result11();
+ MeqResult& result12 = result.result12();
+ MeqResult& result21 = result.result21();
+ MeqResult& result22 = result.result22();
+
+ // Evaluate main value.
+ // Evaluate X polarization.
+ evaluate(request, azel[0].getValue(), azel[1].getValue(),
+ result11.getValueRW(), result12.getValueRW(), itsHeight, itsLength,
+ itsSlant, itsOrientation);
+ // Evaluate Y polarization.
+ evaluate(request, azel[0].getValue(), azel[1].getValue(),
+ result21.getValueRW(), result22.getValueRW(), itsHeight, itsLength,
+ itsSlant, itsOrientation - casa::C::pi_2);
+
+ // Evaluate perturbed values.
+ const MeqParmFunklet *perturbedParm;
+ for(int i = 0; i < request.nspid(); ++i)
+ {
+ // Find out if this perturbed value needs to be computed.
+ if(azel[0].isDefined(i))
+ {
+ perturbedParm = azel[0].getPerturbedParm(i);
+ }
+ else if(azel[1].isDefined(i))
+ {
+ perturbedParm = azel[1].getPerturbedParm(i);
+ }
+ else
+ {
+ continue;
+ }
+
+ // Evaluate X polarization.
+ evaluate(request, azel[0].getValue(), azel[1].getValue(),
+ result11.getPerturbedValueRW(i), result12.getPerturbedValueRW(i),
+ itsHeight, itsLength, itsSlant, itsOrientation);
+
+ // Evaluate Y polarization.
+ evaluate(request, azel[0].getValue(), azel[1].getValue(),
+ result21.getPerturbedValueRW(i), result22.getPerturbedValueRW(i),
+ itsHeight, itsLength, itsSlant, itsOrientation - casa::C::pi_2);
+
+ result11.setPerturbedParm(i, perturbedParm);
+ result12.setPerturbedParm(i, perturbedParm);
+ result21.setPerturbedParm(i, perturbedParm);
+ result22.setPerturbedParm(i, perturbedParm);
+ }
+
+ return result;
+}
+
+
+void MeqDipoleBeam::evaluate(const MeqRequest &request, const MeqMatrix &in_az,
+ const MeqMatrix &in_el, MeqMatrix &out_E_theta, MeqMatrix &out_E_phi,
+ double height, double length, double slant, double orientation)
+{
+ const double *az = in_az.doubleStorage();
+ const double *el = in_el.doubleStorage();
+
+ double *E_theta_re, *E_theta_im;
+ out_E_theta.setDCMat(request.nx(), request.ny());
+ out_E_theta.dcomplexStorage(E_theta_re, E_theta_im);
+
+ double *E_phi_re, *E_phi_im;
+ out_E_phi.setDCMat(request.nx(), request.ny());
+ out_E_phi.dcomplexStorage(E_phi_re, E_phi_im);
+
+ for(int t = 0; t < request.ny(); ++t)
+ {
+ double phi = orientation + (*az);
+ double theta = casa::C::pi_2 - (*el);
+
+ if(theta >= casa::C::_2pi)
+ {
+ // Below the horizon.
+ for(int f = 0; f < request.nx(); ++f)
+ {
+ *E_theta_re = 0.0;
+ *E_theta_im = 0.0;
+ E_theta_re++;
+ E_theta_im++;
+
+ *E_phi_re = 0.0;
+ *E_phi_im = 0.0;
+ E_phi_re++;
+ E_phi_im++;
+ }
+ }
+ else
+ {
+ // Common terms that depend on direction (and therefore on time).
+ double sin_theta = sin(theta);
+ double cos_theta = cos(theta);
+ double sin_phi = sin(phi);
+ double cos_phi = cos(phi);
+ double sin_alpha = sin(slant);
+ double cos_alpha = cos(slant);
+ double tan_alpha = tan(slant);
+ double inv_sin_alpha = 1.0 / sin_alpha;
+ double inv_sin_alpha_sqr = 1.0 / (sin_alpha * sin_alpha);
+
+ double A = sin_theta * cos_phi;
+ double B = cos_theta / tan_alpha;
+ double F = A + B;
+ double G = A - B;
+
+ double I = sin_alpha * sin_phi;
+ double J = cos_alpha * sin_theta;
+ double K = sin_alpha * cos_theta * cos_phi;
+ double M = J + K;
+ double N = J - K;
+
+ double freq = request.domain().startX() + request.stepX() / 2.0;
+
+ for(int f = 0; f < request.nx(); ++f)
+ {
+ // Common terms that depend on frequency.
+ double omega = casa::C::_2pi * freq;
+ double k = omega / casa::C::c;
+
+ double C = k * length;
+ double D = inv_sin_alpha * cos(C);
+ double E = k * height * cos_theta;
+
+ dcomplex term1 = makedcomplex(cos(E), sin(E));
+
+ double term2_F_phase = F * C * sin_alpha;
+ double term2_F_re = cos(term2_F_phase) * inv_sin_alpha - D;
+ double term2_F_im =
+ sin(term2_F_phase) * inv_sin_alpha - F * sin(C);
+ dcomplex term2_F = makedcomplex(term2_F_re, term2_F_im);
+ // k-term vanishes when computing E_phi/E_theta.
+ double term3_F = F * F - inv_sin_alpha_sqr;
+
+ double term2_G_phase = G * C * sin_alpha;
+ double term2_G_re = cos(term2_G_phase) * inv_sin_alpha - D;
+ double term2_G_im =
+ sin(term2_G_phase) * inv_sin_alpha - G * sin(C);
+ dcomplex term2_G = makedcomplex(term2_G_re, term2_G_im);
+ // k-term vanishes when computing E_phi/E_theta.
+ double term3_G = G * G - inv_sin_alpha_sqr;
+
+ dcomplex Gamma1 = (-term1 * term2_G) / term3_G;
+ dcomplex Gamma2 = (-term1 * term2_F) / term3_F;
+ dcomplex Gamma3 = (-conj(term1) * term2_G) / term3_G;
+ dcomplex Gamma4 = (-conj(term1) * term2_F) / term3_F;
+
+ // mu = 1e-7 * 4.0 * pi --> mu / (4 * pi) = 1e-7
+ double H = (1e-7 * inv_sin_alpha) * omega;
+ dcomplex E_phi =
+ H * (N * Gamma1 - M * Gamma2 - N * Gamma3 + M * Gamma4);
+ dcomplex E_theta = H * I * (Gamma1 + Gamma2 - Gamma3 - Gamma4);
+
+ *(E_theta_re) = real(E_theta);
+ *(E_theta_im) = imag(E_theta);
+ ++E_theta_re;
+ ++E_theta_im;
+
+ *(E_phi_re) = real(E_phi);
+ *(E_phi_im) = imag(E_phi);
+ ++E_phi_re;
+ ++E_phi_im;
+
+ freq += request.stepX();
+ }
+ }
+ ++az; ++el;
+ }
+}
+
+
+} //# namespace BBS
+} //# namespace LOFAR
diff --git a/CEP/BB/BBSKernel/src/MNS/MeqExpr.cc b/CEP/BB/BBSKernel/src/MNS/MeqExpr.cc
index adc4826e0c766d2cd4473193ce69da67db47526e..b50fe54528c789397b10c3642454f6e33931da3c 100644
--- a/CEP/BB/BBSKernel/src/MNS/MeqExpr.cc
+++ b/CEP/BB/BBSKernel/src/MNS/MeqExpr.cc
@@ -80,6 +80,12 @@ void MeqExprRep::removeChild(MeqExpr child)
itsChildren.erase(it);
}
+MeqExpr MeqExprRep::getChild(size_t index)
+{
+ ASSERT(index < itsChildren.size());
+ return itsChildren[index];
+}
+
int MeqExprRep::setLevel (int level)
{
if (level < itsMinLevel) itsMinLevel = level;
diff --git a/CEP/BB/BBSKernel/src/MNS/MeqRequest.cc b/CEP/BB/BBSKernel/src/MNS/MeqRequest.cc
index 96748f66cf7f8d37b1a1f0b7e92a9b771a682afe..de0a9fc44c83d473ca29ec397ca8f62bb7925b44 100644
--- a/CEP/BB/BBSKernel/src/MNS/MeqRequest.cc
+++ b/CEP/BB/BBSKernel/src/MNS/MeqRequest.cc
@@ -35,9 +35,9 @@ MeqRequestId MeqRequest::theirRequestId = 0;
MeqRequest::MeqRequest (const MeqDomain& domain, int nx, int ny, int nrSpid)
: itsRequestId (theirRequestId++),
- itsOffset (0, 0),
itsSourceNr (0),
- itsNspids (nrSpid)
+ itsNspids (nrSpid),
+ itsOffset (0, 0)
{
setDomain (domain, nx, ny);
}
@@ -46,9 +46,9 @@ MeqRequest::MeqRequest (const MeqDomain& domain, int nx,
const vector<double>& y,
int nrSpid)
: itsRequestId (theirRequestId++),
- itsOffset (0, 0),
itsSourceNr (0),
- itsNspids (nrSpid)
+ itsNspids (nrSpid),
+ itsOffset (0, 0)
{
setDomain (domain, nx, y);
}
@@ -57,9 +57,9 @@ MeqRequest::MeqRequest (const MeqRequest& req,
uint stx, uint nrx, uint sty, uint nry)
: itsRequestId (req.itsRequestId),
itsStepX (req.itsStepX),
- itsOffset (req.itsOffset),
itsSourceNr (req.itsSourceNr),
- itsNspids (req.itsNspids)
+ itsNspids (req.itsNspids),
+ itsOffset (req.itsOffset)
{
ASSERT (int(stx+nrx) <= req.itsNx && int(sty+nry) <= req.itsNy);
itsNx = nrx;
@@ -78,9 +78,9 @@ MeqRequest::MeqRequest (const MeqRequest& req)
itsStepX (req.itsStepX),
itsY (req.itsY),
itsYP (req.itsYP),
- itsOffset (req.itsOffset),
itsSourceNr (req.itsSourceNr),
- itsNspids (req.itsNspids)
+ itsNspids (req.itsNspids),
+ itsOffset (req.itsOffset)
{
if (itsY.size() > 0) {
itsYP = &itsY[0];
@@ -97,9 +97,9 @@ MeqRequest& MeqRequest::operator= (const MeqRequest& req)
itsStepX = req.itsStepX;
itsY = req.itsY;
itsYP = req.itsYP;
- itsOffset = req.itsOffset;
itsSourceNr = req.itsSourceNr;
itsNspids = req.itsNspids;
+ itsOffset = req.itsOffset;
if (itsY.size() > 0) {
itsYP = &itsY[0];
}
diff --git a/CEP/BB/BBSKernel/src/Makefile.am b/CEP/BB/BBSKernel/src/Makefile.am
index f7832b58d49f93d3cf66c27a1deb1ee844e0ebec..dbbc452f318cb9a3597d585ea8506c4c84823cc7 100644
--- a/CEP/BB/BBSKernel/src/Makefile.am
+++ b/CEP/BB/BBSKernel/src/Makefile.am
@@ -8,10 +8,12 @@ libbbskernel_la_SOURCES = \
Prediffer.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 \
diff --git a/CEP/BB/BBSKernel/src/Model.cc b/CEP/BB/BBSKernel/src/Model.cc
index 81388df57d58bc0a282f61a4d8c824d20c2c9c00..bd5fa25e4dd6b18289d6c8f35025fdf9a99f644f 100644
--- a/CEP/BB/BBSKernel/src/Model.cc
+++ b/CEP/BB/BBSKernel/src/Model.cc
@@ -35,6 +35,8 @@
#include <BBSKernel/MNS/MeqJonesCMul3.h>
#include <BBSKernel/MNS/MeqJonesSum.h>
#include <BBSKernel/MNS/MeqJonesVisData.h>
+#include <BBSKernel/MNS/MeqAzEl.h>
+#include <BBSKernel/MNS/MeqDipoleBeam.h>
#include <BBSKernel/MNS/MeqStation.h>
#include <BBSKernel/MNS/MeqStatUVW.h>
@@ -89,17 +91,14 @@ void Model::makeEquations(EquationType type, const vector<string> &components,
mask[GAIN] = true;
else if(*it == "DIRECTIONAL_GAIN")
mask[DIRECTIONAL_GAIN] = true;
+ else if(*it == "DIPOLE_BEAM")
+ mask[DIPOLE_BEAM] = true;
else if(*it == "BANDPASS")
mask[BANDPASS] = true;
else if(*it == "PHASORS")
mask[PHASORS] = true;
}
-/*
- if(mask{GAIN] && mask[DIRECTIONAL_GAIN])
- LOG_WARN("Model components GAIN and DIRECTIONAL_GAIN are mutually"
- " exclusive; using GAIN only.");
-*/
-
+
string part1("real:");
string part2("imag:");
if(mask[PHASORS])
@@ -121,19 +120,42 @@ void Model::makeEquations(EquationType type, const vector<string> &components,
size_t nStations = itsStationNodes.size();
size_t nSources = itsSourceNodes.size();
+ ASSERT(nStations > 0 && nSources > 0);
+
+ vector<MeqExpr> azel(nSources);
vector<MeqExpr> dft(nStations * nSources);
- vector<MeqJonesExpr> gain, inv_gain, bandpass;
+ vector<MeqJonesExpr> bandpass, gain, dir_gain, dipole_beam;
if(mask[BANDPASS])
+ {
bandpass.resize(nStations);
+ }
+
+ if(mask[GAIN])
+ {
+ gain.resize(nStations);
+ }
+
+ if(mask[DIRECTIONAL_GAIN])
+ {
+ dir_gain.resize(nStations * nSources);
+ }
- if(mask[GAIN] || mask[DIRECTIONAL_GAIN])
+ if(mask[DIPOLE_BEAM])
{
- inv_gain.resize(nStations);
- if(mask[GAIN])
- gain.resize(nStations);
- else
- gain.resize(nStations * nSources);
+ // Make an E-jones expression per source. For now, we use the instrument
+ // position as reference position on earth. In the future, we will need
+ // an E-jones expression per source, _per station_.
+ azel.resize(nSources);
+ dipole_beam.resize(nSources);
+ for(size_t i = 0; i < nSources; ++i)
+ {
+ azel[i] =
+ new MeqAzEl(*(itsSourceNodes[i]), *(itsStationNodes[i].get()));
+
+ dipole_beam[i] = new MeqDipoleBeam(azel[i], 1.706, 1.38,
+ casa::C::pi / 4.001, -casa::C::pi_4);
+ }
}
for(size_t i = 0; i < nStations; ++i)
@@ -153,13 +175,13 @@ void Model::makeEquations(EquationType type, const vector<string> &components,
MeqExpr B22(MeqParmFunklet::create("bandpass:22:" + suffix,
parmGroup, instrumentDBase));
- bandpass[i] = new MeqDiag(MeqExpr(B11), MeqExpr(B22));
+ bandpass[i] = new MeqDiag(B11, B22);
}
// Make a complex gain expression per station and possibly per source.
if(mask[GAIN])
{
- MeqExprRep *J11, *J12, *J21, *J22;
+ MeqExpr J11, J12, J21, J22;
string suffix = itsStationNodes[i]->getName();
// Make a J-jones expression per station
@@ -195,17 +217,17 @@ void Model::makeEquations(EquationType type, const vector<string> &components,
J22 = new MeqExprToComplex(J22_part1, J22_part2);
}
- gain[i] = new MeqJonesNode(MeqExpr(J11), MeqExpr(J12),
- MeqExpr(J21), MeqExpr(J22));
- inv_gain[i] = new MeqJonesInvert(gain[i]);
+ gain[i] = new MeqJonesNode(J11, J12, J21, J22);
+// inv_gain[i] = new MeqJonesInvert(gain[i]);
}
- else if(mask[DIRECTIONAL_GAIN])
+
+ if(mask[DIRECTIONAL_GAIN])
{
// Make a J-jones expression per station per source. Eventually,
// patches of several sources will be supported as well.
for(size_t j = 0; j < nSources; ++j)
{
- MeqExprRep *J11, *J12, *J21, *J22;
+ MeqExpr J11, J12, J21, J22;
string suffix = itsStationNodes[i]->getName() + ":"
+ itsSourceNodes[j]->getName();
@@ -241,32 +263,96 @@ void Model::makeEquations(EquationType type, const vector<string> &components,
J22 = new MeqExprToComplex(J22_part1, J22_part2);
}
- gain[i * nSources + j] = new MeqJonesNode(MeqExpr(J11),
- MeqExpr(J12), MeqExpr(J21), MeqExpr(J22));
+ dir_gain[i * nSources + j] =
+ new MeqJonesNode(J11, J12, J21, J22);
// Gain correction is always performed with respect to the
// direction of the first source (patch).
- if(j == 0)
- inv_gain[i] = new MeqJonesInvert(gain[i * nSources]);
+// if(j == 0)
+// inv_gain[i] = new MeqJonesInvert(gain[i * nSources]);
}
}
}
if(type == CORRECT)
{
+ ASSERTSTR(mask[GAIN] || mask[DIRECTIONAL_GAIN] || mask[DIPOLE_BEAM],
+ "Need at least one of GAIN, DIRECTIONAL_GAIN, or DIPOLE_BEAM to"
+ " correct for.");
+
+ vector<MeqJonesExpr> inv_bandpass, inv_gain, inv_dir_gain;
+ MeqJonesExpr inv_dipole_beam;
+
+ if(mask[BANDPASS])
+ {
+ inv_bandpass.resize(nStations);
+ for(size_t i = 0; i < nStations; ++i)
+ {
+ inv_bandpass[i] = new MeqJonesInvert(bandpass[i]);
+ }
+ }
+
+ if(mask[GAIN])
+ {
+ inv_gain.resize(nStations);
+ for(size_t i = 0; i < nStations; ++i)
+ {
+ inv_gain[i] = new MeqJonesInvert(gain[i]);
+ }
+ }
+
+ if(mask[DIRECTIONAL_GAIN])
+ {
+ inv_dir_gain.resize(nStations);
+ for(size_t i = 0; i < nStations; ++i)
+ {
+ // Always correct for the first source (direction).
+ inv_dir_gain[i] = new MeqJonesInvert(dir_gain[i * nSources]);
+ }
+ }
+
+ if(mask[DIPOLE_BEAM])
+ {
+ // Always correct for the first source (direction).
+ inv_dipole_beam = new MeqJonesInvert(dipole_beam[0]);
+ }
+
for(set<baseline_t>::const_iterator it = baselines.begin();
it != baselines.end();
++it)
{
const baseline_t &baseline = *it;
- itsEquations[baseline] = new MeqJonesCMul3(inv_gain[baseline.first],
- new MeqJonesVisData(buffer, baseline),
+
+ MeqJonesExpr vis = new MeqJonesVisData(buffer, baseline);
+
+ if(mask[BANDPASS])
+ {
+ vis = new MeqJonesCMul3(inv_bandpass[baseline.first], vis,
+ inv_bandpass[baseline.second]);
+ }
+
+ if(mask[GAIN])
+ {
+ vis = new MeqJonesCMul3(inv_gain[baseline.first], vis,
inv_gain[baseline.second]);
+ }
+
+ if(mask[DIRECTIONAL_GAIN])
+ {
+ vis = new MeqJonesCMul3(inv_dir_gain[baseline.first], vis,
+ inv_dir_gain[baseline.second]);
+ }
+
+ if(mask[DIPOLE_BEAM])
+ {
+ vis = new MeqJonesCMul3(inv_dipole_beam, vis, inv_dipole_beam);
+ }
+
+ itsEquations[baseline] = vis;
}
}
else
{
-
for(set<baseline_t>::const_iterator it = baselines.begin();
it != baselines.end();
++it)
@@ -276,23 +362,33 @@ void Model::makeEquations(EquationType type, const vector<string> &components,
vector<MeqJonesExpr> terms;
for(size_t j = 0; j < nSources; ++j)
{
- MeqExpr base
+ // Phase shift (incorporates geometry and fringe stopping).
+ MeqExpr shift
(new MeqBaseDFTPS(dft[baseline.first * nSources + j],
dft[baseline.second * nSources + j], itsLMNNodes[j]));
+ // Point source.
MeqPointSource *source =
dynamic_cast<MeqPointSource*>(itsSourceNodes[j]);
- MeqJonesExpr sourceTerm(new MeqBaseLinPS(base, source));
+ MeqJonesExpr sourceExpr = new MeqBaseLinPS(shift, source);
+
+ // Apply image plane effects if required.
+ if(mask[DIPOLE_BEAM])
+ {
+ sourceExpr = new MeqJonesCMul3(dipole_beam[j], sourceExpr,
+ dipole_beam[j]);
+ }
if(mask[DIRECTIONAL_GAIN])
{
- terms.push_back
- (new MeqJonesCMul3(gain[baseline.first * nSources + j],
- sourceTerm, gain[baseline.second * nSources + j]));
+ sourceExpr = new MeqJonesCMul3
+ (dir_gain[baseline.first * nSources + j],
+ sourceExpr,
+ dir_gain[baseline.second * nSources + j]);
}
- else
- terms.push_back(sourceTerm);
+
+ terms.push_back(sourceExpr);
}
MeqJonesExpr sum;
@@ -301,6 +397,7 @@ void Model::makeEquations(EquationType type, const vector<string> &components,
else
sum = MeqJonesExpr(new MeqJonesSum(terms));
+ // Apply UV-plane effects if required.
if(mask[GAIN])
{
sum = new MeqJonesCMul3(gain[baseline.first], sum,
@@ -312,6 +409,7 @@ void Model::makeEquations(EquationType type, const vector<string> &components,
sum = new MeqJonesCMul3(bandpass[baseline.first], sum,
bandpass[baseline.second]);
}
+
itsEquations[baseline] = sum;
}
}