DirectionCoordinate.h

Go to the documentation of this file.

// # DirectionCoordinate.h: Interconvert pixel positions and directions (e.g. RA/DEC)
// # Copyright (C) 1997,1998,1999,2000,2001,2002,2003,2004
// # Associated Universities, Inc. Washington DC, USA.
// #
// # This library is free software; you can redistribute it and/or modify it
// # under the terms of the GNU Library General Public License as published by
// # the Free Software Foundation; either version 2 of the License, or (at your
// # option) any later version.
// #
// # This library 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 Library General Public
// # License for more details.
// #
// # You should have received a copy of the GNU Library General Public License
// # along with this library; if not, write to the Free Software Foundation,
// # Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
// #
// # Correspondence concerning AIPS++ should be addressed as follows:
// #        Internet email: casa-feedback@nrao.edu.
// #        Postal address: AIPS++ Project Office
// #                        National Radio Astronomy Observatory
// #                        520 Edgemont Road
// #                        Charlottesville, VA 22903-2475 USA
 
#ifndef COORDINATES_DIRECTIONCOORDINATE_H
#define COORDINATES_DIRECTIONCOORDINATE_H
 
#include <casacore/casa/aips.h>
#include <casacore/coordinates/Coordinates/Coordinate.h>
#include <casacore/coordinates/Coordinates/Projection.h>
#include <casacore/casa/Arrays/Vector.h>
#include <casacore/measures/Measures/MDirection.h>
#include <casacore/measures/Measures/MeasConvert.h>
#include <casacore/casa/Quanta/RotMatrix.h>
#include <wcslib/wcs.h>
 
struct celprm;
struct prjprm;
struct wcsprm;
 
namespace casacore {  // # NAMESPACE CASACORE - BEGIN
 
class MVDirection;
class MVAngle;
class LogIO;
template <class T>
class Quantum;
 
// <summary>
// Interconvert pixel positions and directions (e.g. RA/DEC).
// </summary>
 
// <use visibility=export>
 
// <reviewed reviewer="Peter Barnes" date="1999/12/24" tests="tDirectionCoordinate">
// </reviewed>
 
// <prerequisite>
//   <li> Knowledge of astronomical coordinate conversions in general. Probably the
//        best documents are the papers by Mark Calabretta and Eric Greisen.
//        The initial draft from 1996 can be found at
//        http://www.atnf.csiro.au/~mcalabre.  It is this draft that the
//        Coordinate classes are based upon.  Since then, this paper has evolved
//        into three which can be found at the above address, and will be published in the
//        Astronomy and Astrophysics Supplement Series (probably in 2000).
//        The design has changed since the initial draft.  When these papers
//        are finalized, and the IAU has ratified the new standards, WCSLIB
//        (Mark Calabretta's implementation of these conventions) will be
//        revised for the new designs.  At that time, the Coordinate classes
//        may also be revised.
//   <li> <linkto class=Coordinate>Coordinate</linkto> defines the fundamental
//        interface to coordinate conversions.
//   <li> <linkto class=MDirection>MDirection</linkto> defines the types of
//        directions (J2000 etc.) which are defined. The measures machinery
//        also implements "astronomical" conversions which are outside the
//        scope of these coordinates (for example, <src>J2000</src> to
//        <src>B1950</src>).
//   <li> <linkto class=Projection>Projection</linkto> defines the types of
//        celestial projections which are available.
// </prerequisite>
//
// <synopsis>
// This class implements pixel to world coordinate conversions. This class
// implements geometric conversions (e.g. SIN projection) via the WCS library
// and also provides an interface to astronomical conversions (RA/DEC <--> l,b)
// via the <linkto module=Measures>Measures</linkto> module.
// </synopsis>
//
//
// <note role=caution>
// All absolute pixels coordinates are zero relative.
// </note>
//
// <example>
// Let's make a DirectionCoordinate --- used to represent a direction,
// usually an RA/DEC, but it could also be, e.g., an AZ/EL pair.
// <srcblock>
//    Matrix<Double> xform(2,2);                                    // 1
//    xform = 0.0; xform.diagonal() = 1.0;                          // 2
//    DirectionCoordinate radec(MDirection::J2000,                  // 3
//                            Projection(Projection::SIN),          // 4
//                            135*C::pi/180.0, 60*C::pi/180.0,      // 5
//                            -1*C::pi/180.0, 1*C::pi/180,          // 6
//                            xform,                                // 7
//                            128, 128);                            // 8
// </srcblock>
// <ul>
//    <li> <i>1-2:</i>Here we set up a diagonal transformation matrix.
//         Normally this matrix should be diagonal, however if you wanted
//         to introduce a rotation or skew, you would do it through this
//         matrix.
//    <li> <i>3:</i>This defines the astronomical type of the world
//         coordinate. Most of the time it will probably be J2000
//         or B1950, but many other possibilities are possible as listed
//         in the <linkto class=MDirection>MDirection</linkto> class
//         header.
//    <li> <i>4:</i>The <linkto class=Projection>Projection</linkto> class
//         defines the "geometry" that is used to map <src>xy<-->world</src>. SIN
//         is the most common projection for radio interferometers. Note that
//         SIN can optionally take parameters as defined in Calabretta and Greisen.
//         If not provided, they default to 0.0, which is the "old" SIN
//         convention.
//    <li> <i>5:</i>Set the reference position to RA=135, DEC=60 degrees.
//         Note that the native units of a Direction is radians.
//    <li> <i>6:</i> Set the increments to -1 degree in RA, and +1 degree
//         in DEC.
//    <li> <i>7:</i> Set the previously defined transformation matrix.
//    <li> <i>8:</i> Set the zero-relative reference pixel. Note that it does
//         not have to be incremental. At the reference pixel, the world
//         coordinate has the reference value.
// </ul>
//
// In this example it is more convenient to change the units to degrees. This can
// be accomplished as follows:
// <srcblock>
//    Vector<String> units(2); units = "deg";                       //  9
//    radec.setWorldAxisUnits(units);                               // 10
// </srcblock>
// The increment and reference value are updated appropriately.
//
// Set up a couple of vectors to use the world and pixel coordinate values.
// <srcblock>
//    Vector<Double> world(2), pixel(2);                            // 11
//    pixel = 138.0;                                                // 12
// </srcblock>
// We use 138 as an arbitrary pixel position which is near the reference pixel
// so we can tell if the answers look foolish or not.
// We can actually perform a transformation like this as follows. If
// it succeeds we print the value of the world coordinate.
// <srcblock>
//    Bool ok = radec.toWorld(world, pixel);                        // 13
//    if (!ok) {                                                    // 14
//      cout << "Error: " << radec.errorMessage() << endl;          // 15
//      return 1;                                                   // 16
//    }                                                             // 17
//    cout << world << " <--- " << pixel << endl;         // 18
// </srcblock>
// There is an overloaded "toWorld" function that produces an MDirection
// in case you want to, e.g., find out what the position in B1950 coordinates
// would be.
//
// The reverse transformation takes place similarly:
// <srcblock>
//    ok = radec.toPixel(pixel, world);                             // 19
// </srcblock>
// </example>
//
// <example>
// We could also have made the above DirectionCoordinate using the Quantum-based
// constructor, which is a little more elegant if you want to use degrees.
//
//    Matrix<Double> xform(2,2);
//    xform = 0.0; xform.diagonal() = 1.0;
//    Quantum<Double> refLon(135.0, "deg");
//    Quantum<Double> refLat(60.0, "deg");
//    Quantum<Double> incLon(-1.0, "deg");
//    Quantum<Double> incLat(1.0, "deg");
//    DirectionCoordinate radec(MDirection::J2000,
//                            Projection(Projection::SIN),
//                            refLon, refLat,
//                            incLon, incLat,
//                            xform,
//                            128, 128);
//
// But note that the constructor will have converted the native units
// of the DirectionCoordinate to radians.  So the Double-based toWorld and
// toPixel functions will be in terms of radians.   If you want the native
// units to be degrees, then again you can use
//
// <srcblock>
//    Vector<String> units(2); units = "deg";
//    radec.setWorldAxisUnits(units);
// </srcblock>
// and thereafter degrees are the native units.
// </example>
//
// <motivation>
// Directions in the sky are fundamental to astronomy.
// </motivation>
//
//
// <thrown>
//   <li>  AipsError
// </thrown>
//
// <todo asof="2000/01/01">
//   <li> Nothing
// </todo>
 
class DirectionCoordinate : public Coordinate {
 public:
  // The default constructor creates a J2000 DirectionCoordinate with a
  // CARtesion projection with longitude,latitude 0,0 at pixel 0,0 and an
  // increment of +1 radian per pixel on both axes.
  DirectionCoordinate();
 
  // Define the DirectionCoordinate transformation. <src>refLong</src> and
  // <src>refLat</src> will normally the the RA/DEC of the pixel described by
  // <src>refX/refY</src>. <src>incLat/incLong</src>
  // are the increments per pixel (RA is usually negative), and the <src>xform</src>
  // matrix is usually the unit diagonal matrix unless you have a rotation or
  // some other linear transformation between the pixel and world axes.
  //
  // Note that the units are radians initially. You can change it to degrees
  // or something else with the <src>setWorldAxisUnits</src> method later if you want.
  //
  // longPole and latPole are defined by Calabretta and Greisen (these
  // are reference points not at the native pole).  In general
  // you can leave these out and the default values will cause them
  // to be computed appropriately.  However, when reading from FITS
  // the LONPOLE and LATPOLE keywords are passed along here.
  DirectionCoordinate(MDirection::Types directionType, const Projection& projection, Double refLong,
                      Double refLat, Double incLong, Double incLat, const Matrix<Double>& xform,
                      Double refX, Double refY, Double longPole = 999.0, Double latPole = 999.0);
 
  // Create DirectionCoordinate with Quantum-based interface.
  // Parameters are the same as above.
  // Regardless of the units of the quanta, the initial units
  // of the DirectionCoordinate will be converted radians.
  // You can change it to degrees or something else with the
  // setWorldAxisUnits method later if you want.
  //
  // longPole and latPole are defined by Calabretta and Greisen (these
  // are reference points not at the native pole).  In general
  // you can leave these out and the default values will cause them
  // to be computed appropriately.  However, when reading from FITS
  // the LONPOLE and LATPOLE keywords are passed along here.
  // To get the default the 999.0 value should be used (units
  // are irrelevant in that case)
  DirectionCoordinate(MDirection::Types directionType, const Projection& projection,
                      const Quantum<Double>& refLong, const Quantum<Double>& refLat,
                      const Quantum<Double>& incLong, const Quantum<Double>& incLat,
                      const Matrix<Double>& xform, Double refX, Double refY,
                      const Quantum<Double>& longPole = Quantum<Double>(999.0, Unit("rad")),
                      const Quantum<Double>& latPole = Quantum<Double>(999.0, Unit("rad")));
 
  // Constructor from WCS structure; must hold ONLY a celestial wcs structure
  // Specify whether the absolute pixel coordinates in the wcs structure
  // are 0- or 1-relative.  The coordinate is always constructed with 0-relative
  // pixel coordinates
  DirectionCoordinate(MDirection::Types directionType, const ::wcsprm& wcs, Bool oneRel = True);
 
  // Copy constructor (copy semantics)
  DirectionCoordinate(const DirectionCoordinate& other);
 
  // Assignment (copy semantics).
  DirectionCoordinate& operator=(const DirectionCoordinate& other);
 
  // Destructor
  virtual ~DirectionCoordinate();
 
  // Return Coordinate::DIRECTION
  virtual Coordinate::Type type() const;
 
  // Always returns the String "Direction".
  virtual String showType() const;
 
  // Always returns 2.
  // <group>
  virtual uInt nPixelAxes() const;
  virtual uInt nWorldAxes() const;
  // </group>
 
  // Set extra conversion type.  Whenever a conversion from pixel to world is done,
  // the world value is then further converted to this MDirection::Types value.
  // For example, your DirectionCoordinate may be defined in J2000.
  // You can use this to get the world values out in say GALACTIC.
  // Similarly, whenever you convert from world to pixel, the world
  // value is assumed to be that appropriate to the conversionDirectionType.
  // It is first converted to the MDirection::Types with which the
  // DirectionCoordinate was constructed and from there to pixel.
  // If you don't call this function, or you set the same type
  // for which the DirectionCoordinate was constructed, no extra
  // conversions occur.   Some conversions will fail.  These are the
  // ones that require extra frame information (epoch, position) such
  // as to AZEL from J2000 etc.  This will be added later.
  //
  // In the mixed pixel/world conversion routine <src>toMix</src>
  // the implementation is only partial.  See the comments for this
  // function below.
  // <group>
  void setReferenceConversion(MDirection::Types type);
  void getReferenceConversion(MDirection::Types& type) const { type = conversionType_p; };
  // </group>
 
  // Convert a pixel position to a world position or vice versa. Returns True
  // if the conversion succeeds, otherwise it returns False and method
  // errorMessage returns its error message.
  // The output vectors are appropriately resized.
  // if <src>useConversionFrame</src>, if the coordinate has a conversion
  // layer frame, it is used. Else, the native frame is used for the conversion.
  // <group>
  virtual Bool toWorld(Vector<Double>& world, const Vector<Double>& pixel,
                       Bool useConversionFrame = True) const;
 
  // <src>world</src> values must have units equivalent to the world axis
  // units. If the coordinate has a conversion layer, the world coordinates
  // must be supplied in the conversion frame.
  virtual Bool toPixel(Vector<Double>& pixel, const Vector<Double>& world) const;
  // </group>
 
  // Mixed pixel/world coordinate conversion.
  // <src>worldIn</src> and <src>worldAxes</src> are of length
  // nWorldAxes.
  // <src>pixelIn</src> and <src>pixelAxes</src> are of length nPixelAxes.
  // <src>worldAxes(i)=True</src> specifies you have given a world
  // value in <src>worldIn(i)</src> to convert to pixel.
  // <src>pixelAxes(i)=True</src> specifies you have given a pixel
  // value in <src>pixelIn(i)</src> to convert to world.
  // You cannot specify the same axis via <src>worldAxes</src>
  // and <src>pixelAxes</src>.
  // Values in <src>pixelIn</src> are converted to world and
  // put into <src>worldOut</src> in the appropriate world axis
  // location.  Values in <src>worldIn</src> are copied to
  // <src>worldOut</src>.
  // Values in <src>worldIn</src> are converted to pixel and
  // put into <src>pixelOut</src> in the appropriate pixel axis
  // location.  Values in <src>pixelIn</src> are copied to
  // <src>pixelOut</src>.
  //
  // <src>worldMin</src> and <src>worldMax</src> specify the range of the world
  // coordinate (in the world axis units of that world axis
  // in the CoordinateSystem) being solved for in a mixed calculation
  // for each world axis.    Some mixed solutions can be degenerate, whereupon you
  // you must say which one you want.  Use functions <src>setWorldMixRanges</src>
  // and <src>worldMixMin, worldMixMax</src> to set these ranges,
  // If you don't know, use the defaults (function <src>setDefaultWorldMixRanges</src>.
  // Removed axes are handled (for example, a removed pixel
  // axis with remaining corresponding world axis will
  // correctly be converted to world using the replacement
  // value).
  // Returns True if the conversion succeeds, otherwise it returns False and
  // <src>errorMessage()</src> contains an error message. The output vectors
  // are resized.
  //
  // If you actually request a pure pixel to world or world to pixel
  // via <src>toMix</src>, then the functions <src>toWorld</src> or <src>toPixel</src>
  // will be invoked directly (see above) and the extra conversion layer
  // invoked through function <src>setReferenceConversion</src> will be active.
  // However, if you request a true mixed pixel/world conversion,
  // the extra conversion layer is not activated (because of the nature of mixed
  // conversions).  This situation may change in the future
  // with a partial implementation added.
  virtual Bool toMix(Vector<Double>& worldOut, Vector<Double>& pixelOut,
                     const Vector<Double>& worldIn, const Vector<Double>& pixelIn,
                     const Vector<Bool>& worldAxes, const Vector<Bool>& pixelAxes,
                     const Vector<Double>& worldMin, const Vector<Double>& worldMax) const;
 
  // Compute and retrieve the world min and max ranges, for use in function <src>toMix</src>,
  // for  a lattice of the given shape (for this coordinate).   Using these
  // ranges in <src>toMix</src> should speed it up and help avoid ambiguity.
  // If the shape is negative, that indicates that the shape is unknown
  // for that axis.  The default range is used for that axis.  This situation
  // arises in a CoordinateSystem for which a pixel, but not a world axis
  // has been removed.
  // The output vectors are resized.  Returns False if fails (and
  // then <src>setDefaultWorldMixRanges</src> generates the ranges)
  // with a reason in <src>errorMessage()</src>.
  // The <src>setDefaultWorldMixRanges</src> function
  // just gives you [-90->90], [-180,180] (in appropriate units)
  // <group>
  virtual Bool setWorldMixRanges(const IPosition& shape);
  virtual void setDefaultWorldMixRanges();
  // </group>
 
  // Non-virtual function.  When <src>which</src> is T, use the
  // world value as the center for the mix world range.
  void setWorldMixRanges(const Vector<Bool>& which, const Vector<Double>& world);
 
  // A convenient way to turn the world vector into an MDirection or MVDirection
  // for further processing in the Measures system.
  // <br>We could improve the performance of this if it would be useful. However it is
  // expected that normally one would just call this once to get a template
  // MDirection, and then call the vector versions.
  // <br>In case of a failure, the versions with a Bool return value will return
  // False. The other versions will throw an exception.
  // <group>
  Bool toWorld(MDirection& world, const Vector<Double>& pixel) const;
  Bool toPixel(Vector<Double>& pixel, const MDirection& world) const;
  Bool toWorld(MVDirection& world, const Vector<Double>& pixel) const;
  Bool toPixel(Vector<Double>& pixel, const MVDirection& world) const;
  MVDirection toWorld(const Vector<Double>& pixel) const;
  Vector<Double> toPixel(const MVDirection& world) const;
  Vector<Double> toPixel(const MDirection& world) const;
  //</group>
 
  // Batch up a lot of transformations. The first (most rapidly varying) axis
  // of the matrices contain the coordinates. Returns False if any conversion
  // failed  and  <src>errorMessage()</src> will hold a message.
  // The <src>failures</src> array is the length of the number of conversions
  // (True for failure, False for success)
  // <group>
  virtual Bool toWorldMany(Matrix<Double>& world, const Matrix<Double>& pixel,
                           Vector<Bool>& failures) const;
  virtual Bool toPixelMany(Matrix<Double>& pixel, const Matrix<Double>& world,
                           Vector<Bool>& failures) const;
  // </group>
 
  // Make absolute world coordinates relative and vice-versa (relative to
  // the reference value).  Note that these functions are independent
  // of the MDirection::Types  (set either at construction or by function
  // <src>setReferenceConversion</src>).  The vectors must be
  // of length <src>nWorldAxes</src> or memory access errors will occur
  //<group>
  virtual void makeWorldRelative(Vector<Double>& world) const;
  virtual void makeWorldRelative(MDirection& world) const;
  virtual void makeWorldAbsolute(Vector<Double>& world) const;
  virtual void makeWorldAbsolute(MDirection& world) const;
  //</group>
 
  // Make absolute coordinates relative and vice versa with respect
  // to the given reference value.  Add the other functions in this grouping
  // as needed.
  //<group>
  virtual void makeWorldAbsoluteRef(Vector<Double>& world, const Vector<Double>& refVal) const;
  //</group>
 
  // Recover the requested attribute.
  // <group>
  MDirection::Types directionType(Bool showConversion = False) const;
  Projection projection() const;
  virtual Vector<String> worldAxisNames() const;
  virtual Vector<String> worldAxisUnits() const;
  virtual Vector<Double> referenceValue() const;
  virtual Vector<Double> increment() const;
  virtual Matrix<Double> linearTransform() const;
  virtual Vector<Double> referencePixel() const;
  // </group>
 
  // Set the value of the requested attribute.  Note that these just
  // change the internal values, they do not cause any recomputation.
  // <group>
  virtual Bool setWorldAxisNames(const Vector<String>& names);
  virtual Bool setReferencePixel(const Vector<Double>& refPix);
  virtual Bool setLinearTransform(const Matrix<Double>& xform);
  virtual Bool setIncrement(const Vector<Double>& inc);
  virtual Bool setReferenceValue(const Vector<Double>& refval);
  // </group>
 
  // Change the world axis units.  Adjust the increment and
  // reference value by the ratio of the old and new units.
  // The units must be compatible with
  // angle. The units are initially "rad" (radians).
  virtual Bool setWorldAxisUnits(const Vector<String>& units);
 
  // Return canonical axis names for the given MDirection type,
  // giving FITS names if desired.
  // BEG think this should be in the MDirection class, but WNB
  // disagrees. Leave it here for now.
  static Vector<String> axisNames(MDirection::Types type, Bool FITSName = False);
 
  // Comparison function. Any private Double data members are compared
  // with the specified fractional tolerance.  Don't compare on the specified
  // axes in the Coordinate.  If the comparison returns False,  method
  // errorMessage returns a message about why.
  // <group>
  virtual Bool near(const Coordinate& other, Double tol = 1e-6) const;
  virtual Bool near(const Coordinate& other, const Vector<Int>& excludeAxes,
                    Double tol = 1e-6) const;
  // </group>
 
  // Format a DirectionCoordinate coordinate world value nicely through the
  // common format interface.  See <linkto class=Coordinate>Coordinate</linkto>
  // for basics.
  //
  // Formatting types that are allowed are SCIENTIFIC, FIXED, MIXED, and TIME
  // If you ask for format type Coordinate::DEFAULT then the
  // selected format depends upon what the value of the enum
  // MDirection::GlobalTypes is for this DirectionCoordinate.
  // For example, if it is GRADEC or GHADEC you would
  // get Coordinate::TIME style formatting (DD:MM:SS.SS), otherwise
  // you would get Coordinate::FIXED formatting by default.
  //
  // <src>axis</src> says which axis in this Coordinate we are formatting.
  // We have to know this because we may format Longitude and Latitude differently.
  // For Coordinate::TIME style formatting, precision
  // refers to the places after the decimal in the SS field.
  //
  // If you leave <src>units</src> empty, then it makes up a nice unit for you.
  //<group>
  virtual void getPrecision(Int& precision, Coordinate::formatType& format, Bool showAsAbsolute,
                            Int defPrecScientific, Int defPrecFixed, Int defPrecTime) const;
  virtual String format(String& units, Coordinate::formatType format, Double worldValue, uInt axis,
                        Bool isAbsolute, Bool showAsAbsolute, Int precision = -1,
                        Bool usePrecForMixed = False) const;
  //</group>
 
  // Fix cylindrical coordinates to put the longitude in [-180,180] range.
  // If False returned, it failed an an error is in <src>errorMessage</src>
  // This fix is not done automatically internally because of the dependence
  // on the image shape.  It should be called for any foreign image
  // (such as FITS) that is imported
  Bool cylindricalFix(Int shapeLong, Int shapeLat);
 
  // Find the Coordinate for when we Fourier Transform ourselves.  This pointer
  // must be deleted by the caller. Axes specifies which axes of the Coordinate
  // you wish to transform.   Shape specifies the shape of the image
  // associated with all the axes of the Coordinate.   Currently the
  // output reference pixel is always shape/2. If the pointer returned is 0,
  // it failed with a message in <src>errorMessage</src>
  virtual Coordinate* makeFourierCoordinate(const Vector<Bool>& axes,
                                            const Vector<Int>& shape) const;
 
  // Save the DirectionCoordinate into the supplied record using the supplied field name.
  // The field must not exist, otherwise <src>False</src> is returned.
  virtual Bool save(RecordInterface& container, const String& fieldName) const;
 
  // Recover the DirectionCoordinate from a record.
  // A null pointer means that the restoration did not succeed.
  static DirectionCoordinate* restore(const RecordInterface& container, const String& fieldName);
 
  // Make a copy of the DirectionCoordinate using new. The caller
  // is responsible for calling delete.
  virtual Coordinate* clone() const;
 
  // Fish out the ref and non-native poles (refLong, refLat, longPole, latPole)
  // Not for general use.  Units are degrees.
  Vector<Double> longLatPoles() const;
 
  // get the pixel area.
  Quantity getPixelArea() const;
 
  // Convert this coordinate to another reference frame by rotating it about
  // the reference pixel so the the axes of the new reference frame are
  // aligned along the current pixel axes. The reference pixel remains the
  // same and the conversion is exact for the reference pixel and in general
  // becomes less accurate as distance from reference pixel increases. The
  // latitude like and the longitude like pixel increments are preserved.
  // Conversions which require extra information such as epoch and position
  // are not supported. The <src>angle</src> parameter is the angle between
  // the new coordinate and the pixel coordinate, measured clockwise from the
  // positive y-axis of the new coordinate to the positive y-axis of the pixel
  // coordinate; ie, it is the clockwise angle through which the current world
  // coordinate would have to be rotated so that the new coordinate's axes
  // would be parallel to the pixel axes. The accuracy of the returned angle
  // is good to at least 7 digits.
  DirectionCoordinate convert(Quantity& angle, MDirection::Types directionType) const;
 
  // Set the projection.
  void setProjection(const Projection&);
 
  // Set the base (as opposed to conversion) reference frame.
  void setReferenceFrame(const MDirection::Types rf);
 
  // Are the pixels square?
  Bool hasSquarePixels() const;
 
  // Is the projection equivalent to NCP?
  Bool isNCP() const;
 
 private:
  // Direction type
  MDirection::Types type_p, conversionType_p;
 
  // Projection parameters
  Projection projection_p;
 
  // WCS structure.  This is mutable because the wcs functions
  // that do toPixel and toWorld (which have const signature)
  // require a non const wcs structure.  so either all of these
  // virtual functions lose their const or we use mutable...
  mutable ::wcsprm wcs_p;
 
  // WCS computes in degrees - use this to convert back and forth between
  // current DirectionCoordinate units and degrees or radians
  Vector<Double> to_degrees_p;  // From current units to degrees
  Vector<Double> to_radians_p;  // From current units to radians
 
  // Axis names.
  Vector<String> names_p;
 
  // Current units.
  Vector<String> units_p;
 
  // Rotation matrix used to handle relative coordinates
  RotMatrix rot_p;
 
  // Conversion machines.
  // "To"   handles type_p -> conversionType_p
  // "From" handles conversionType_p -> type_p;
  mutable MDirection::Convert* pConversionMachineTo_p;
  mutable MDirection::Convert* pConversionMachineFrom_p;
 
  // Interconvert between the current units and wcs units (degrees)
  // <group>
  void toCurrent(Vector<Double>& degrees) const;
  void fromCurrent(Vector<Double>& current) const;
  // </group>
 
  // Check formatting types.
  void checkFormat(Coordinate::formatType& format, Bool absolute) const;
 
  // Format a latitude.
  String formatLatitude(String& units, MVAngle& mVA, Bool absolute, Coordinate::formatType form,
                        Int prec) const;
  // Format a longitude.
  String formatLongitude(String& units, MVAngle& mVA, MDirection::GlobalTypes gtype, Bool absolute,
                         Coordinate::formatType form, Int prec) const;
 
  // Mixed pixel/world coordinate conversion.  Vector in must
  // be length nWorldAxes (2).  Specify whether longitude
  // (in(0)) or latitude (in(1)) is the world coordinate . It is
  // assumed that the other value is the pixel coordinate.
  Bool toMix2(Vector<Double>& out, const Vector<Double>& in, const Vector<Double>& minWorld,
              const Vector<Double>& maxWorld, Bool longIsWorld) const;
 
  // Initialize unit conversion vectors and units
  void initializeFactors();
 
  // Helper functions interfacing to WCS.
  // <group>
  void makeDirectionCoordinate(MDirection::Types directionType, const Projection& proj,
                               Double refLong, Double refLat, Double incLong, Double incLat,
                               const Matrix<Double>& xform, Double refX, Double refY,
                               Double longPole, Double latPole);
  //
  void makeWCS(::wcsprm& wcs, const Matrix<Double>& xform, const Projection& proj,
               MDirection::Types directionType, Double refPixLong, Double refPixLat, Double refLong,
               Double refLat, Double incLong, Double incLat, Double longPole, Double latPole);
  // </group>
 
  // Normalize each row of the PC matrix such that increment() will return the actual
  // angular increment and any scale factors are removed from the PC matrix
  // (modifies wcs_p.pc _and_ wcs_p.cdelt _and_ wcs_p.altlin,
  // executes set_wcs() and hence wcsset() on the struct)
  // See Greisen & Calabretta, A&A 395, 1061-1075 (2002), equation (4)
  void normalizePCMatrix();
 
  Double putLongInPiRange(Double lon, const String& unit) const;
 
  // Set up conversion machine
  void makeConversionMachines();
 
  // Convert from type_p -> conversionType_p
  // <group>
  virtual void convertTo(Vector<Double>& world) const;
  virtual void convertFrom(Vector<Double>& world) const;
  // </group>
 
  // Copy private data
  void copy(const DirectionCoordinate& other);
 
  // Set up the offset coordinate rotation matrix.  Units
  // of long and lat are current world units
  // <group>
  void setRotationMatrix();
  void setRotationMatrix(RotMatrix& rot, Double lon, Double lat) const;
  // </group>
 
  // Return unit conversion vector for converting to current units
  const Vector<Double> toCurrentFactors() const;
};
 
}  // namespace casacore
 
#endif