From 591526fc71ce1b9ef260d948a97c7b84d07b74ae Mon Sep 17 00:00:00 2001
From: Bas van der Tol <vdtol@strw.leidenuniv.nl>
Date: Wed, 25 Sep 2013 21:24:32 +0000
Subject: [PATCH] Task #2782: update ionospheric phase screen in BBS to conform
implementation in awimager
---
.../include/BBSKernel/Expr/PiercePoint.h | 6 +-
.../include/BBSKernel/IonosphereExpr.h | 6 +-
.../BBSKernel/MeasurementExprLOFARUtil.h | 2 +-
.../BBSKernel/src/Expr/IonPhaseShift.cc | 67 ++---
.../BBSKernel/src/Expr/PiercePoint.cc | 234 +++++++-----------
.../BBSKernel/src/IonosphereExpr.cc | 8 +-
.../BBSKernel/src/MeasurementExprLOFAR.cc | 25 +-
.../BBSKernel/src/MeasurementExprLOFARUtil.cc | 4 +-
.../BBSKernel/src/StationExprLOFAR.cc | 21 +-
CEP/Calibration/ExpIon/CMakeLists.txt | 10 +-
10 files changed, 152 insertions(+), 231 deletions(-)
diff --git a/CEP/Calibration/BBSKernel/include/BBSKernel/Expr/PiercePoint.h b/CEP/Calibration/BBSKernel/include/BBSKernel/Expr/PiercePoint.h
index 1e239b484dd..41886dc2988 100644
--- a/CEP/Calibration/BBSKernel/include/BBSKernel/Expr/PiercePoint.h
+++ b/CEP/Calibration/BBSKernel/include/BBSKernel/Expr/PiercePoint.h
@@ -45,19 +45,19 @@ namespace BBS
// \addtogroup Expr
// @{
-class PiercePoint: public BasicBinaryExpr<Vector<2>, Scalar, Vector<4> >
+class PiercePoint: public BasicBinaryExpr<Vector<3>, Scalar, Vector<4> >
{
public:
typedef shared_ptr<PiercePoint> Ptr;
typedef shared_ptr<const PiercePoint> ConstPtr;
PiercePoint(const casa::MPosition &position,
- const Expr<Vector<2> >::ConstPtr &azel,
+ const Expr<Vector<3> >::ConstPtr &direction,
const Expr<Scalar>::ConstPtr &height);
protected:
virtual const Vector<4>::View evaluateImpl(const Grid &grid,
- const Vector<2>::View &azel, const Scalar::View &height) const;
+ const Vector<3>::View &direction, const Scalar::View &height) const;
private:
// Station position in ITRF coordinates.
diff --git a/CEP/Calibration/BBSKernel/include/BBSKernel/IonosphereExpr.h b/CEP/Calibration/BBSKernel/include/BBSKernel/IonosphereExpr.h
index f891fcea684..8a6a6090731 100644
--- a/CEP/Calibration/BBSKernel/include/BBSKernel/IonosphereExpr.h
+++ b/CEP/Calibration/BBSKernel/include/BBSKernel/IonosphereExpr.h
@@ -62,7 +62,7 @@ public:
virtual ~IonosphereExpr();
virtual Expr<JonesMatrix>::Ptr construct(const casa::MPosition &refPosition,
- const casa::MPosition &station, const Expr<Vector<2> >::ConstPtr &azel)
+ const casa::MPosition &station, const Expr<Vector<3> >::ConstPtr &direction)
const = 0;
};
@@ -75,7 +75,7 @@ public:
MIMExpr(const IonosphereConfig &config, Scope &scope);
virtual Expr<JonesMatrix>::Ptr construct(const casa::MPosition &refPosition,
- const casa::MPosition &station, const Expr<Vector<2> >::ConstPtr &azel)
+ const casa::MPosition &station, const Expr<Vector<3> >::ConstPtr &direction)
const;
private:
@@ -92,7 +92,7 @@ public:
ExpIonExpr(const IonosphereConfig&, Scope &scope);
virtual Expr<JonesMatrix>::Ptr construct(const casa::MPosition &refPosition,
- const casa::MPosition &station, const Expr<Vector<2> >::ConstPtr &azel)
+ const casa::MPosition &station, const Expr<Vector<3> >::ConstPtr &direction)
const;
private:
diff --git a/CEP/Calibration/BBSKernel/include/BBSKernel/MeasurementExprLOFARUtil.h b/CEP/Calibration/BBSKernel/include/BBSKernel/MeasurementExprLOFARUtil.h
index ff52a9201a7..46bc71833e0 100644
--- a/CEP/Calibration/BBSKernel/include/BBSKernel/MeasurementExprLOFARUtil.h
+++ b/CEP/Calibration/BBSKernel/include/BBSKernel/MeasurementExprLOFARUtil.h
@@ -138,7 +138,7 @@ makeScalarPhaseExpr(Scope &scope,
Expr<JonesMatrix>::Ptr
makeIonosphereExpr(const Station::ConstPtr &station,
const casa::MPosition &refPosition,
- const Expr<Vector<2> >::Ptr &exprAzEl,
+ const Expr<Vector<3> >::Ptr &exprDirection,
const IonosphereExpr::Ptr &exprIonosphere);
// Right multiply \p lhs by \p rhs. Return \p rhs if \p lhs is uninitialized.
diff --git a/CEP/Calibration/BBSKernel/src/Expr/IonPhaseShift.cc b/CEP/Calibration/BBSKernel/src/Expr/IonPhaseShift.cc
index af0ae825c2b..5ccee50ae4b 100644
--- a/CEP/Calibration/BBSKernel/src/Expr/IonPhaseShift.cc
+++ b/CEP/Calibration/BBSKernel/src/Expr/IonPhaseShift.cc
@@ -160,13 +160,14 @@ const JonesMatrix IonPhaseShift::evaluateExpr(const Request &request,
// Compute (differential) total electron content (TEC) at the piercepoint
// (see memo by Bas van der Tol).
+ ASSERT(!r0.value().isArray());
Matrix r0sqr = r0.value() * r0.value();
+ ASSERT(!beta.value().isArray());
Matrix beta_2 = beta.value() * 0.5;
// LOG_DEBUG_STR("r0sqr: " << r0sqr);
// LOG_DEBUG_STR("beta: " << beta_2);
- Matrix tecX(0.0, 1, nTime);
- Matrix tecY(0.0, 1, nTime);
+ Matrix tec(0.0, nFreq, nTime);
for(size_t i = 0; i < calPiercePoint.size(); ++i)
{
Matrix dx = calPiercePoint[i].value(0) - piercePoint.value(0);
@@ -174,43 +175,34 @@ const JonesMatrix IonPhaseShift::evaluateExpr(const Request &request,
Matrix dz = calPiercePoint[i].value(2) - piercePoint.value(2);
Matrix weight = pow((dx * dx + dy * dy + dz * dz) / r0sqr, beta_2);
-// LOG_DEBUG_STR("dx: " << dx);
-// LOG_DEBUG_STR("dy: " << dy);
-// LOG_DEBUG_STR("dz: " << dz);
+// LOG_DEBUG_STR("dx: " << dx.getDouble());
+// LOG_DEBUG_STR("dy: " << dy.getDouble());
+// LOG_DEBUG_STR("dz: " << dz.getDouble());
// LOG_DEBUG_STR("weight: " << weight);
- tecX += weight * tecWhite[i].value(0);
- tecY += weight * tecWhite[i].value(1);
+ tec += weight * tecWhite[i].value(0);
}
- tecX *= -0.5;
- tecY *= -0.5;
+ tec *= -0.5;
// Correct TEC value for the slant (angle of the line of sight w.r.t. the
// normal to the ionospheric thin layer at the pierce point position). A
// large slant implies a longer path through the ionosphere, and thus a
// higher TEC value.
- tecX /= cos(piercePoint.value(3));
- tecY /= cos(piercePoint.value(3));
-
+ tec *= piercePoint.value(3);
+
// Convert from slanted TEC to phase shift.
//
// TODO: Because MeqMatrix cannot handle the elementwise product of an N x 1
// times a 1 x M array, we have to write it out ourselves. Maybe this could
// be made possible in MeqMatrix instead?
- Matrix shiftX, shiftY;
- shiftX.setDCMat(nFreq, nTime);
- shiftY.setDCMat(nFreq, nTime);
- double *shiftX_re, *shiftX_im;
- double *shiftY_re, *shiftY_im;
- shiftX.dcomplexStorage(shiftX_re, shiftX_im);
- shiftY.dcomplexStorage(shiftY_re, shiftY_im);
-
- ASSERT(!tecX.isComplex() && tecX.isArray()
- && static_cast<size_t>(tecX.nelements()) == nTime);
- const double *tecItX = tecX.doubleStorage();
- ASSERT(!tecY.isComplex() && tecY.isArray()
- && static_cast<size_t>(tecY.nelements()) == nTime);
- const double *tecItY = tecY.doubleStorage();
+ Matrix shift;
+ shift.setDCMat(nFreq, nTime);
+ double *shift_re, *shift_im;
+ shift.dcomplexStorage(shift_re, shift_im);
+
+ ASSERT(!tec.isComplex() && tec.isArray()
+ && static_cast<size_t>(tec.nelements()) == nFreq*nTime);
+ const double *tecIt = tec.doubleStorage();
for(size_t t = 0; t < nTime; ++t)
{
@@ -230,30 +222,25 @@ const JonesMatrix IonPhaseShift::evaluateExpr(const Request &request,
double k_nu = 8.44797245e9 / reqGrid[FREQ]->center(f);
- double phaseX = k_nu * (*tecItX);
- double phaseY = k_nu * (*tecItY);
+ double phase = k_nu * (*tecIt);
- *shiftX_re++ = std::cos(phaseX);
- *shiftX_im++ = std::sin(phaseX);
- *shiftY_re++ = std::cos(phaseY);
- *shiftY_im++ = std::sin(phaseY);
+ *shift_re++ = std::cos(phase);
+ *shift_im++ = std::sin(phase);
+
+ ++tecIt;
}
- ++tecItX;
- ++tecItY;
}
-// LOG_DEBUG_STR("TEC X: " << tecX);
-// LOG_DEBUG_STR("TEC Y: " << tecY);
-// LOG_DEBUG_STR("SHIFT X: " << shiftX);
-// LOG_DEBUG_STR("SHIFT Y: " << shiftY);
+// LOG_DEBUG_STR("TEC: " << tec);
+// LOG_DEBUG_STR("SHIFT: " << shift);
JonesMatrix result;
result.setFlags(mergeFlags(flags.begin(), flags.end()));
- result.assign(0, 0, shiftX);
+ result.assign(0, 0, shift);
result.assign(0, 1, Matrix(makedcomplex(0.0, 0.0)));
result.assign(1, 0, Matrix(makedcomplex(0.0, 0.0)));
- result.assign(1, 1, shiftY);
+ result.assign(1, 1, shift);
EXPR_TIMER_STOP();
diff --git a/CEP/Calibration/BBSKernel/src/Expr/PiercePoint.cc b/CEP/Calibration/BBSKernel/src/Expr/PiercePoint.cc
index 47ff1d1b3ce..3e884ec3f2e 100644
--- a/CEP/Calibration/BBSKernel/src/Expr/PiercePoint.cc
+++ b/CEP/Calibration/BBSKernel/src/Expr/PiercePoint.cc
@@ -44,66 +44,41 @@ namespace LOFAR
namespace BBS
{
+const double earth_ellipsoid_a = 6378137.0;
+const double earth_ellipsoid_a2 = earth_ellipsoid_a*earth_ellipsoid_a;
+const double earth_ellipsoid_b = 6356752.3142;
+const double earth_ellipsoid_b2 = earth_ellipsoid_b*earth_ellipsoid_b;
+const double earth_ellipsoid_e2 = (earth_ellipsoid_a2 - earth_ellipsoid_b2) / earth_ellipsoid_a2;
+
+
PiercePoint::PiercePoint(const casa::MPosition &position,
- const Expr<Vector<2> >::ConstPtr &direction,
+ const Expr<Vector<3> >::ConstPtr &direction,
const Expr<Scalar>::ConstPtr &height)
- : BasicBinaryExpr<Vector<2>, Scalar, Vector<4> >(direction, height),
+ : BasicBinaryExpr<Vector<3>, Scalar, Vector<4> >(direction, height),
itsPosition(casa::MPosition::Convert(position,
casa::MPosition::ITRF)())
{
- // Get geodetic longitude, geodetic lattitude, and height above the
- // ellipsoid (WGS84).
- casa::MPosition mPositionWGS84 =
- casa::MPosition::Convert(itsPosition, casa::MPosition::WGS84)();
- const casa::MVPosition &mvPositionWGS84 = mPositionWGS84.getValue();
- itsLon = mvPositionWGS84.getLong();
- itsLat = mvPositionWGS84.getLat();
-
- // Compute the distance from the center of gravity of the earth to the
- // station position.
- itsPositionRadius = itsPosition.getValue().getLength().getValue();
-
- // Use height above the ellipsoid to compute the earth radius at the
- // position of the station.
- //
- // TODO: You cannot just subtract the height above the ellipsoid from the
- // length of the ITRF position vector and expect to get the earth radius at
- // the station position. Assuming "earth radius" means distance from the
- // center of mass of the earth to the position on the ellipsoid specified by
- // itsLat, itsLon, then still this is incorrect because the position vector
- // is generally not parallel to the ellipsoidal normal vector at the
- // position on the ellipsoid at itsLat, itsLon (along which the height above
- // the ellipsoid is measured).
- //
- // TODO: This earth radius is specific to the station position, yet it is
- // also used when computing the length of the pierce point vector (which
- // intersects the earth's surface at a different position, where the earth
- // radius may be different!).
- itsEarthRadius = itsPositionRadius - mvPositionWGS84.getLength().getValue();
-
-// LOG_DEBUG_STR("Antenna: Longitude: " << (itsLon * 180.0) / casa::C::pi
-// << " deg, Lattitude: " << (itsLat * 180.0) / casa::C::pi
-// << " deg, Height: " << itsHeight << " m, Radius: " << radius
-// << " m, Earth radius: " << itsEarthRadius << " m");
-// LOG_DEBUG_STR("Antenna: Longitude: " << itsLon
-// << " rad, Lattitude: " << itsLat
-// << " rad, Height: " << itsHeight << " m, Radius: " << radius
-// << " m, Earth radius: " << itsEarthRadius << " m");
}
const Vector<4>::View PiercePoint::evaluateImpl(const Grid &grid,
- const Vector<2>::View &azel, const Scalar::View &height) const
+ const Vector<3>::View &direction, const Scalar::View &height) const
{
const size_t nTime = grid[TIME]->size();
+ const size_t nFreq = grid[FREQ]->size();
// Allocate space for the result.
// TODO: This is a hack! The Matrix class does not support 1xN or Nx1
// "matrices".
- Matrix out_x, out_y, out_z, out_alpha;
- double *x = out_x.setDoubleFormat(1, nTime);
- double *y = out_y.setDoubleFormat(1, nTime);
- double *z = out_z.setDoubleFormat(1, nTime);
- double *alpha = out_alpha.setDoubleFormat(1, nTime);
+ Matrix out_x, out_y, out_z, out_airmass;
+ double *out_x_ptr = out_x.setDoubleFormat(nFreq, nTime);
+ double *out_y_ptr = out_y.setDoubleFormat(nFreq, nTime);
+ double *out_z_ptr = out_z.setDoubleFormat(nFreq, nTime);
+ double *out_airmass_ptr = out_airmass.setDoubleFormat(nFreq, nTime);
+ double pp_x, pp_y, pp_z, pp_airmass;
+
+
+ ASSERT(!height().isArray());
+ double h = height().getDouble();
// Get station position in ITRF coordinates.
const casa::MVPosition &mPosition = itsPosition.getValue();
@@ -111,109 +86,90 @@ const Vector<4>::View PiercePoint::evaluateImpl(const Grid &grid,
double stationY = mPosition(1);
double stationZ = mPosition(2);
-// LOG_DEBUG_STR("Geocentric lattitude: " << std::atan2(stationZ, std::sqrt(stationX * stationX + stationY * stationY)));
-
- // Precompute rotation matrix to transform from local ENU (East, North, Up)
- // coordinates to ITRF coordinates (see e.g.
- // http://en.wikipedia.org/wiki/Geodetic_system).
- //
- // TODO: Use measures to compute a direction vector in ITRF directly,
- // instead of going from (RA, Dec) to (azimuth, elevation), to ENU to ITRF.
- double sinlon = std::sin(itsLon);
- double coslon = std::cos(itsLon);
- double sinlat = std::sin(itsLat);
- double coslat = std::cos(itsLat);
- double R[3][3] = {{-sinlon, -sinlat * coslon, coslat * coslon},
- { coslon, -sinlat * sinlon, coslat * sinlon},
- { 0.0, coslat, sinlat}};
-
-// {
-// // A first stab at using the measures to compute a direction vector
-// // in ITRF coordinates directly. Needs more work.
-// casa::Quantity qEpoch(grid[TIME]->center(0), "s");
-// casa::MEpoch mEpoch(qEpoch, casa::MEpoch::UTC);
-
-// // Create and initialize a frame.
-// casa::MeasFrame frame;
-// frame.set(itsPosition);
-// frame.set(mEpoch);
-
-// // Create conversion engine.
-// casa::MDirection mDirection(casa::MVDirection(4.33961, 1.09537),
-// casa::MDirection::Ref(casa::MDirection::J2000));
-
-// casa::MDirection::Convert converter =
-// casa::MDirection::Convert(mDirection,
-// casa::MDirection::Ref(casa::MDirection::ITRF, frame));
-
-// LOG_DEBUG_STR("Convertor: " << converter);
-
-// // Compute XYZ direction vector.
-// casa::MVDirection mvXYZ(converter().getValue());
-// const casa::Vector<casa::Double> xyz = mvXYZ.getValue();
-// LOG_DEBUG_STR("XYZ direction (from measures): " << xyz(0) << " "
-// << xyz(1) << " " << xyz(2));
-// }
-
- ASSERT(azel(0).isArray() && azel(0).nx() == 1);
- ASSERT(azel(1).isArray() && azel(1).nx() == 1);
- ASSERT(!height().isArray());
+ const double ion_ellipsoid_a = earth_ellipsoid_a + h;
+ const double ion_ellipsoid_a2_inv = 1.0 / (ion_ellipsoid_a * ion_ellipsoid_a);
+ const double ion_ellipsoid_b = earth_ellipsoid_b + h;
+ const double ion_ellipsoid_b2_inv = 1.0 / (ion_ellipsoid_b * ion_ellipsoid_b);
- // TODO: itsEarthRadius is valid only for the station position, which is
- // generally not equal to the position where the pierce point vector
- // intersects the ellipsoid. The question is how the height of the
- // ionospheric layer should be defined. A layer with a fixed distance from
- // the earth center of mass is probably easiest.
- double ionosphereRadius = itsEarthRadius + height().getDouble();
+ double x = stationX/ion_ellipsoid_a;
+ double y = stationY/ion_ellipsoid_a;
+ double z = stationZ/ion_ellipsoid_b;
+ double c = x*x + y*y + z*z - 1.0;
for(size_t i = 0; i < nTime; ++i)
{
- double az = azel(0).getDouble(0, i);
- double el = azel(1).getDouble(0, i);
-
- // Calculate alpha, this is the angle of the line of sight with the
- // normal of the ionospheric layer at the pierce point location (i.e.
- // alpha' (alpha prime) in the memo).
- *alpha =
- std::asin(std::cos(el) * (itsPositionRadius / ionosphereRadius));
-
- // Compute direction vector in local ENU (East, North, Up) coordinates.
- double cosel = std::cos(el);
- double dx = std::sin(az) * cosel;
- double dy = std::cos(az) * cosel;
- double dz = std::sin(el);
-
- // Transform the direction vector from the local ENU frame to the ITRF
- // frame.
- double dxr = R[0][0] * dx + R[0][1] * dy + R[0][2] * dz;
- double dyr = R[1][0] * dx + R[1][1] * dy + R[1][2] * dz;
- double dzr = R[2][0] * dx + R[2][1] * dy + R[2][2] * dz;
-
-// if(i == 0)
-// {
-// LOG_DEBUG_STR("XYZ direction (from az, el): " << dxr << " " << dyr
-// << " " << dzr);
-// }
-
- // Compute the distance from the station to the pierce point, i.e.
- // equation 6 in the memo. The equation below is equivalent to that in
- // the memo, but it is expressed in elevation instead of zenith angle.
- double radius = ionosphereRadius * std::cos(el + *alpha) / cosel;
-
- // Compute the pierce point location in ITRF coordinates.
- *x = stationX + radius * dxr;
- *y = stationY + radius * dyr;
- *z = stationZ + radius * dzr;
-
- ++x; ++y; ++z; ++alpha;
+ double directionX = direction(0).getDouble(0, i);
+ double directionY = direction(1).getDouble(0, i);
+ double directionZ = direction(2).getDouble(0, i);
+
+ double dx = directionX / ion_ellipsoid_a;
+ double dy = directionY / ion_ellipsoid_a;
+ double dz = directionZ / ion_ellipsoid_b;
+
+ double a = dx*dx + dy*dy + dz*dz;
+ double b = x*dx + y*dy + z*dz;
+ double alpha = (-b + std::sqrt(b*b - a*c))/a;
+
+ pp_x = stationX + alpha*directionX;
+ pp_y = stationY + alpha*directionY;
+ pp_z = stationZ + alpha*directionZ;
+ double normal_x = pp_x * ion_ellipsoid_a2_inv;
+ double normal_y = pp_y * ion_ellipsoid_a2_inv;
+ double normal_z = pp_z * ion_ellipsoid_b2_inv;
+ double norm_normal2 = normal_x*normal_x + normal_y*normal_y + normal_z*normal_z;
+ double norm_normal = std::sqrt(norm_normal2);
+ double sin_lat2 = normal_z*normal_z / norm_normal2;
+
+ double g = 1.0 - earth_ellipsoid_e2*sin_lat2;
+ double sqrt_g = std::sqrt(g);
+
+ double M = earth_ellipsoid_b2 / ( earth_ellipsoid_a * g * sqrt_g );
+ double N = earth_ellipsoid_a / sqrt_g;
+
+ double local_ion_ellipsoid_e2 = (M-N) / ((M+h)*sin_lat2 - N - h);
+ double local_ion_ellipsoid_a = (N+h) * std::sqrt(1.0 - local_ion_ellipsoid_e2*sin_lat2);
+ double local_ion_ellipsoid_b = local_ion_ellipsoid_a*std::sqrt(1.0 - local_ion_ellipsoid_e2);
+
+ double z_offset = ((1.0-earth_ellipsoid_e2)*N + h - (1.0-local_ion_ellipsoid_e2)*(N+h)) * std::sqrt(sin_lat2);
+
+ double x1 = stationX/local_ion_ellipsoid_a;
+ double y1 = stationY/local_ion_ellipsoid_a;
+ double z1 = (stationZ-z_offset)/local_ion_ellipsoid_b;
+ double c1 = x1*x1 + y1*y1 + z1*z1 - 1.0;
+
+ dx = directionX / local_ion_ellipsoid_a;
+ dy = directionY / local_ion_ellipsoid_a;
+ dz = directionZ / local_ion_ellipsoid_b;
+ a = dx*dx + dy*dy + dz*dz;
+ b = x1*dx + y1*dy + z1*dz;
+ alpha = (-b + std::sqrt(b*b - a*c1))/a;
+
+ pp_x = stationX + alpha*directionX;
+ pp_y = stationY + alpha*directionY;
+ pp_z = stationZ + alpha*directionZ;
+
+ normal_x = pp_x / (local_ion_ellipsoid_a * local_ion_ellipsoid_a);
+ normal_y = pp_y / (local_ion_ellipsoid_a * local_ion_ellipsoid_a);
+ normal_z = (pp_z-z_offset) / (local_ion_ellipsoid_b * local_ion_ellipsoid_b);
+
+ norm_normal2 = normal_x*normal_x + normal_y*normal_y + normal_z*normal_z;
+ norm_normal = std::sqrt(norm_normal2);
+
+ pp_airmass = norm_normal / (directionX*normal_x + directionY*normal_y + directionZ*normal_z);
+
+ for(size_t j = 0; j < nFreq; ++j) {
+ *(out_x_ptr++) = pp_x;
+ *(out_y_ptr++) = pp_y;
+ *(out_z_ptr++) = pp_z;
+ *(out_airmass_ptr++) = pp_airmass;
+ }
}
-
// Create result.
Vector<4>::View result;
result.assign(0, out_x);
result.assign(1, out_y);
result.assign(2, out_z);
- result.assign(3, out_alpha);
+ result.assign(3, out_airmass);
return result;
}
diff --git a/CEP/Calibration/BBSKernel/src/IonosphereExpr.cc b/CEP/Calibration/BBSKernel/src/IonosphereExpr.cc
index 97d665e136f..3e1ed42e6f4 100644
--- a/CEP/Calibration/BBSKernel/src/IonosphereExpr.cc
+++ b/CEP/Calibration/BBSKernel/src/IonosphereExpr.cc
@@ -111,10 +111,10 @@ MIMExpr::MIMExpr(const IonosphereConfig &config, Scope &scope)
}
Expr<JonesMatrix>::Ptr MIMExpr::construct(const casa::MPosition &refPosition,
- const casa::MPosition &station, const Expr<Vector<2> >::ConstPtr &azel)
+ const casa::MPosition &station, const Expr<Vector<3> >::ConstPtr &direction)
const
{
- PiercePoint::Ptr piercePoint(new PiercePoint(station, azel, itsHeight));
+ PiercePoint::Ptr piercePoint(new PiercePoint(station, direction, itsHeight));
PolynomialLayer::Ptr shift(new PolynomialLayer(refPosition, piercePoint,
itsCoeff.begin(), itsCoeff.end()));
@@ -166,10 +166,10 @@ ExpIonExpr::ExpIonExpr(const IonosphereConfig&, Scope &scope)
}
Expr<JonesMatrix>::Ptr ExpIonExpr::construct(const casa::MPosition&,
- const casa::MPosition &station, const Expr<Vector<2> >::ConstPtr &azel)
+ const casa::MPosition &station, const Expr<Vector<3> >::ConstPtr &direction)
const
{
- PiercePoint::Ptr piercePoint(new PiercePoint(station, azel, itsHeight));
+ PiercePoint::Ptr piercePoint(new PiercePoint(station, direction, itsHeight));
IonPhaseShift::Ptr shift(new IonPhaseShift(piercePoint, itsR0, itsBeta));
shift->setCalibratorPiercePoints(itsCalPiercePoint.begin(),
diff --git a/CEP/Calibration/BBSKernel/src/MeasurementExprLOFAR.cc b/CEP/Calibration/BBSKernel/src/MeasurementExprLOFAR.cc
index 02b613bf7eb..eb2ae9af1ba 100644
--- a/CEP/Calibration/BBSKernel/src/MeasurementExprLOFAR.cc
+++ b/CEP/Calibration/BBSKernel/src/MeasurementExprLOFAR.cc
@@ -243,13 +243,13 @@ void MeasurementExprLOFAR::makeForwardExpr(SourceDB &sourceDB,
{
// Create an AZ, EL expression for the centroid direction of the
// patch.
- Expr<Vector<2> >::Ptr exprAzEl =
- makeAzElExpr(instrument->station(j)->position(),
- exprPatch->position());
+// Expr<Vector<2> >::Ptr exprAzEl =
+// makeAzElExpr(instrument->station(j)->position(),
+// exprPatch->position());
exprDDE[j] = compose(exprDDE[j],
makeIonosphereExpr(instrument->station(j),
- instrument->position(), exprAzEl, exprIonosphere));
+ instrument->position(), exprPatchPositionITRF, exprIonosphere));
}
}
@@ -549,15 +549,9 @@ void MeasurementExprLOFAR::makeInverseExpr(SourceDB &sourceDB,
// Ionosphere.
if(config.useIonosphere())
{
- // Create an AZ, EL expression for the phase reference
- // direction.
- Expr<Vector<2> >::Ptr exprAzEl =
- makeAzElExpr(instrument->station(i)->position(),
- exprRefPhase);
-
stationExpr[i] = compose(stationExpr[i],
makeIonosphereExpr(instrument->station(i),
- instrument->position(), exprAzEl, exprIonosphere));
+ instrument->position(), exprRefPhaseITRF, exprIonosphere));
}
}
}
@@ -583,7 +577,6 @@ void MeasurementExprLOFAR::makeInverseExpr(SourceDB &sourceDB,
makeDirectionalGainExpr(itsScope,
instrument->station(i), patch, config.usePhasors()));
}
-
// Beam.
if(config.useBeam())
{
@@ -629,15 +622,9 @@ void MeasurementExprLOFAR::makeInverseExpr(SourceDB &sourceDB,
// Ionosphere.
if(config.useIonosphere())
{
- // Create an AZ, EL expression for the centroid direction of
- // the patch.
- Expr<Vector<2> >::Ptr exprAzEl =
- makeAzElExpr(instrument->station(i)->position(),
- exprPatch->position());
-
stationExpr[i] = compose(stationExpr[i],
makeIonosphereExpr(instrument->station(i),
- instrument->position(), exprAzEl, exprIonosphere));
+ instrument->position(), exprPatchPositionITRF, exprIonosphere));
}
}
}
diff --git a/CEP/Calibration/BBSKernel/src/MeasurementExprLOFARUtil.cc b/CEP/Calibration/BBSKernel/src/MeasurementExprLOFARUtil.cc
index bf7ab794df1..4020dd42f30 100644
--- a/CEP/Calibration/BBSKernel/src/MeasurementExprLOFARUtil.cc
+++ b/CEP/Calibration/BBSKernel/src/MeasurementExprLOFARUtil.cc
@@ -397,11 +397,11 @@ makeScalarPhaseExpr(Scope &scope,
Expr<JonesMatrix>::Ptr
makeIonosphereExpr(const Station::ConstPtr &station,
const casa::MPosition &refPosition,
- const Expr<Vector<2> >::Ptr &exprAzEl,
+ const Expr<Vector<3> >::Ptr &exprDirection,
const IonosphereExpr::Ptr &exprIonosphere)
{
return exprIonosphere->construct(refPosition, station->position(),
- exprAzEl);
+ exprDirection);
}
Expr<JonesMatrix>::Ptr
diff --git a/CEP/Calibration/BBSKernel/src/StationExprLOFAR.cc b/CEP/Calibration/BBSKernel/src/StationExprLOFAR.cc
index e12e0b57f78..249e2cee4bd 100644
--- a/CEP/Calibration/BBSKernel/src/StationExprLOFAR.cc
+++ b/CEP/Calibration/BBSKernel/src/StationExprLOFAR.cc
@@ -134,7 +134,7 @@ void StationExprLOFAR::initialize(SourceDB &sourceDB, const BufferMap &buffers,
" correction can only be applied for a single direction on the"
" sky.");
}
-
+
// Beam reference position on the sky.
Expr<Vector<2> >::Ptr exprRefDelay = makeDirectionExpr(refDelay);
Expr<Vector<3> >::Ptr exprRefDelayITRF =
@@ -187,15 +187,9 @@ void StationExprLOFAR::initialize(SourceDB &sourceDB, const BufferMap &buffers,
// Ionosphere.
if(config.useIonosphere())
{
- // Create an AZ, EL expression for the phase reference
- // direction.
- Expr<Vector<2> >::Ptr exprAzEl =
- makeAzElExpr(instrument->station(i)->position(),
- exprRefPhase);
-
itsExpr[i] = compose(itsExpr[i],
makeIonosphereExpr(instrument->station(i),
- instrument->position(), exprAzEl, exprIonosphere));
+ instrument->position(), exprRefPhaseITRF, exprIonosphere));
}
}
}
@@ -221,7 +215,6 @@ void StationExprLOFAR::initialize(SourceDB &sourceDB, const BufferMap &buffers,
makeDirectionalGainExpr(itsScope,
instrument->station(i), patch, config.usePhasors()));
}
-
// Beam.
if(config.useBeam())
{
@@ -266,15 +259,9 @@ void StationExprLOFAR::initialize(SourceDB &sourceDB, const BufferMap &buffers,
// Ionosphere.
if(config.useIonosphere())
{
- // Create an AZ, EL expression for the centroid direction of
- // the patch.
- Expr<Vector<2> >::Ptr exprAzEl =
- makeAzElExpr(instrument->station(i)->position(),
- exprPatch->position());
-
itsExpr[i] = compose(itsExpr[i],
makeIonosphereExpr(instrument->station(i),
- instrument->position(), exprAzEl, exprIonosphere));
+ instrument->position(), exprPatchPositionITRF, exprIonosphere));
}
}
}
@@ -477,7 +464,7 @@ PatchExprBase::Ptr StationExprLOFAR::makePatchExpr(const string &name,
it->second));
}
- return PatchExprBase::Ptr(new PatchExpr(itsScope, sourceDB, name,
+ return PatchExprBase::Ptr(new PatchExpr(itsScope, sourceDB, name,
refPhase));
}
diff --git a/CEP/Calibration/ExpIon/CMakeLists.txt b/CEP/Calibration/ExpIon/CMakeLists.txt
index 4c3e74c5daf..d76bc8a1b6f 100644
--- a/CEP/Calibration/ExpIon/CMakeLists.txt
+++ b/CEP/Calibration/ExpIon/CMakeLists.txt
@@ -1,10 +1,14 @@
-# $Id: CMakeLists.txt 14609 2009-12-07 09:10:48Z loose $
+# $Id$
+
+include(LofarFindPackage)
+include(LofarSphinx)
-# Do not split the following line, otherwise makeversion will fail!
lofar_package(ExpIon 1.0 DEPENDS pyparameterset pyparmdb)
-include(LofarFindPackage)
lofar_find_package(Pyrap REQUIRED)
lofar_find_package(Boost REQUIRED COMPONENTS python thread)
lofar_find_package(Casacore REQUIRED COMPONENTS scimath)
+
add_subdirectory(src)
+
+lofar_add_sphinx_doc(sphinx-doc)
--
GitLab