Task #2782: update ionospheric phase screen in BBS to conform implementation in awimager

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.

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:

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.

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();
 

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;
 }

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(),

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));
                 }
             }
         }

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

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));
 }
 

CEP/Calibration/ExpIon/CMakeLists.txt

-# $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)