SpectralCoordinate.h

Go to the documentation of this file.

// # SpectralCoordinate.h: Interconvert between pixel and frequency.
// # 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_SPECTRALCOORDINATE_H
#define COORDINATES_SPECTRALCOORDINATE_H
 
#include <casacore/casa/aips.h>
#include <casacore/casa/Arrays/Vector.h>
#include <casacore/coordinates/Coordinates/Coordinate.h>
#include <casacore/coordinates/Coordinates/ObsInfo.h>
#include <casacore/measures/Measures/MFrequency.h>
#include <casacore/measures/Measures/MDoppler.h>
#include <casacore/measures/Measures/MDirection.h>
#include <casacore/measures/Measures/MPosition.h>
#include <casacore/measures/Measures/MEpoch.h>
#include <casacore/casa/Quanta/Quantum.h>
#include <memory>
 
#include <wcslib/wcs.h>
 
namespace casacore {  // # NAMESPACE CASACORE - BEGIN
 
class TabularCoordinate;
class LogIO;
class MVFrequency;
class VelocityMachine;
template <class T>
class Quantum;
 
// <summary>
// Interconvert pixel and frequency values.
// </summary>
 
// <use visibility=export>
 
// <reviewed reviewer="Peter Barnes" date="1999/12/24" tests="tSpectralCoordinate">
// </reviewed>
//
// <prerequisite>
//   <li> <linkto class=Coordinate>Coordinate</linkto>
//   <li> <linkto class=MFrequency>MFrequency</linkto>,
//        <linkto class=MDoppler>MDoppler</linkto> and
//        <linkto class=VelocityMachine>VelocityMachine</linkto>
//        classes if you want radial velocities.
// </prerequisite>
//
// <synopsis>
// This class performs the mapping from pixel to frequency.
// This can be done via a Tabular lookup or via an algorithmic
// implementation which may be linear or non-linear.  The latter
// is implemented via the WCS library.
//
// </synopsis>
//
 
// <note role=caution>
// All pixels coordinates are zero relative.
// </note>
//
// <example>
// Let us make a linear SpectralCoordinate first
// <srcblock>
//   Double restfreq = 1.420405752E9;
//   Double crpix = 10.0;
//   Double crval = 1.4e9;
//   Double cdelt = 1.0e6;
//   SpectralCoordinate sc(MFrequency::TOPO, crval, cdelt, crpix, restfreq);
//
//   Double world, pixel;
//   pixel = 12.1;
//   if (!sc.toWorld(world, pixel)) {
//      cerr << "Error : " << sc.errorMessage() << endl;
//   } else {
//      cerr << "pixel, world = " << pixel << ", " << world << endl;
//   }
//
// </srcblock>
// </example>
//
// <example>
// Now we make a non-linear SpectralCoordinate
// <srcblock>
//   Vector<Double> freqs(5);
//   freqs(0) = 1.4e9; freqs(1) = 1.41e9;
//   freqs(2) = 1.43e9; freqs(3) = 1.44e9;
//   freqs(4) = 1.47e9;
//   SpectralCoordinate sc(MFrequency::LSRK, freqs, restfreq);
//
//   Double world, pixel;
//   world = 1.42e9;
//   if (!sc.toPixel(pixel, world)) {
//      cerr << "Error : " << sc.errorMessage() << endl;
//   } else {
//      cerr << "world, pixel = " << world << ", " << pixel << endl;
//   }
//
// </srcblock>
// </example>
//
// <motivation>
//  Spectral-line astronomy requires a specialized SpectralCoordinate.
// </motivation>
 
// <todo asof="2000/01/01">
//   <li> Allow other than linear interpolations for frequency lookup.
// </todo>
//
 
class SpectralCoordinate : public Coordinate {
 public:
  enum SpecType {  // taken from the FITS spectral coordinate type codes
    FREQ,
    VRAD,
    VOPT,
    BETA,
    WAVE,
    AWAV
  };
 
  // Default constructor.    It is equivalent to doing
  // SpectralCoordinate(MFrequency::TOPO, 0.0, 1.0, 0.0)
  SpectralCoordinate();
 
  // Create a linear frequency axis SpectralCoordinate
  // <src>f0</src> is the frequency of the reference pixel, <src>inc</src> is the pixel increment,
  // <src>refPix</src> is the reference pixel. You can
  // optionally store the rest frequency for later use in calculating radial
  // velocities.  Use 0 for restFrequency if continuum.
  //
  // Frequencies and increments initially in Hz.
  SpectralCoordinate(MFrequency::Types type, Double f0, Double inc, Double refPix,
                     Double restFrequency = 0.0);
 
  // Create linear frequency axis SpectralCoordinate with Quantum-based interface.
  // Parameters are the same as above.
  // Regardless of the units of the Quanta, the initial units
  // of the SpectralCoordinate will be Hz.  You can change it to
  // something else with the setWorldAxisUnits method later if you want.
  // Use 0 for restFrequency if continuum.
  SpectralCoordinate(MFrequency::Types type, const Quantum<Double>& f0, const Quantum<Double>& inc,
                     Double refPix,
                     const Quantum<Double>& restFrequency = Quantum<Double>(0.0, "Hz"));
 
  // Construct a SpectralCoordinate with the specified frequencies (in Hz).
  // This axis can be nonlinear; the increments and related
  // functions return the <src>average</src> values
  // (calculated from the first and last pixel's frequencies).
  //
  // A linear interpolation/extrapolation is used for pixels which are
  // not supplied. The reference pixel is chosen to be 0.
  // The frequencies must increase or decrease monotonically (otherwise
  // the toPixel lookup would not be possible).
  // Use 0 for restFrequency if continuum.
  SpectralCoordinate(MFrequency::Types type, const Vector<Double>& freqs,
                     Double restFrequency = 0.0);
 
  // Construct a SpectralCoordinate with the specified frequencies
  // with Quantum-based interface.
  // Parameters are the same as above.
  // Regardless of the units of the Quanta, the initial units
  // of the SpectralCoordinate will be Hz.
  // Use 0 for restFrequency if continuum.
  SpectralCoordinate(MFrequency::Types type, const Quantum<Vector<Double>>& freqs,
                     const Quantum<Double>& restFrequency = Quantum<Double>(0.0, "Hz"));
 
  // Construct a SpectralCoordinate with the specified velocities (in km/s).
  // They will be converted to Hz and the SpectralCoordinate constructed.
  // This axis can be nonlinear; the increments and related
  // functions return the <src>average</src> values
  // (calculated from the first and last pixel's frequencies).
  //
  // A linear interpolation/extrapolation is used for pixels which are
  // not supplied. The reference pixel is chosen to be 0.
  // The velocities must increase or decrease monotonically (otherwise
  // the toPixel lookup would not be possible).
  SpectralCoordinate(MFrequency::Types freqType, MDoppler::Types velType,
                     const Vector<Double>& velocities, const String& velUnit,
                     Double restFrequency = 0.0);
 
  // Construct a SpectralCoordinate with the specified wavelengths (in mm).
  // They will be converted to Hz and the SpectralCoordinate constructed.
  // This axis can be nonlinear; the increments and related
  // functions return the <src>average</src> values
  // (calculated from the first and last pixel's frequencies).
  // If inAir is True, the input wavelengths are assumed to be Air Wavelengths.
  // They are converted to vacuum frequency using the refractive index
  // which is calculated based on the mean input air wavelength.
  //
  // A linear interpolation/extrapolation is used for pixels which are
  // not supplied. The reference pixel is chosen to be 0.
  // The wavelengths must increase or decrease monotonically (otherwise
  // the toPixel lookup would not be possible).
  SpectralCoordinate(MFrequency::Types freqType, const Vector<Double>& wavelengths,
                     const String& waveUnit, Double restFrequency = 0.0, Bool inAir = False);
 
  // Construct from wcs structure.  Must hold only a spectral 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
  SpectralCoordinate(MFrequency::Types freqType, const ::wcsprm& wcs, Bool oneRel = True);
 
  // Copy constructor (copy semantics).
  SpectralCoordinate(const SpectralCoordinate& other);
 
  // Assignment (copy semantics).
  SpectralCoordinate& operator=(const SpectralCoordinate& other);
 
  // Destructor.
  virtual ~SpectralCoordinate();
 
  // Always returns Coordinate::SPECTRAL.
  virtual Coordinate::Type type() const;
 
  // Always returns the String "Spectral".
  virtual String showType() const;
 
  // Always returns 1.
  // <group>
  virtual uInt nPixelAxes() const;
  virtual uInt nWorldAxes() const;
  // </group>
 
  // Set extra conversion layer.  Whenever a conversion from pixel to world is done,
  // the world value is then further converted to this MFrequency::Types value.
  // For example, your SpectralCoordinate may be defined in LSRK.
  // You can use this to get the world values out in say BARY. You must
  // specify the position on earth, the epoch and the direction for the conversions
  // and it is your responsibility to ensure they are viable.
  // Similarly, whenever you convert from world to pixel, the world
  // value is assumed to be that appropriate to the setReferenceConversion type.
  // It is first converted to the MFrequency::Types with which the
  // SpectralCoordinate was constructed and from there to pixel.
  // If you don't call this function, or you set the same type
  // for which the SpectralCoordinate was constructed, no extra
  // conversions occur.   Some conversions will fail.  These are the
  // ones that require extra frame information (radial velocity) such
  // as to REST. This will be added later.  In this case this function
  // returns False (and the conversion parameters are all left as they were),
  //  else it returns True.
  // <group>
  Bool setReferenceConversion(MFrequency::Types type, const MEpoch& epoch,
                              const MPosition& position, const MDirection& direction);
  void getReferenceConversion(MFrequency::Types& type, MEpoch& epoch, MPosition& position,
                              MDirection& direction) const {
    type = conversionType_p;
    epoch = epoch_p;
    position = position_p;
    direction = direction_p;
  };
  // </group>
 
  // Convert a pixel to a world coordinate or vice versa. Returns True
  // if the conversion succeeds, otherwise it returns False and
  // <src>errorMessage()</src> contains an error message.  The input vectors
  // must be of length one and the output vectors are resized if they are not
  // already of length one.
  // 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;
  virtual Bool toPixel(Vector<Double>& pixel, const Vector<Double>& world) const;
  Bool toWorld(Double& world, const Double& pixel) const;
  Bool toPixel(Double& pixel, const Double& world) const;
  // </group>
 
  // Convert a pixel (channel number) into an MFrequency or MVFrequency and vice
  // versa. Usually you will do
  // this for calculating velocities or converting frequencies from one frame
  // to another.
  // <group>
  Bool toWorld(MFrequency& world, Double pixel) const;
  Bool toPixel(Double& pixel, const MFrequency& world) const;
  Bool toWorld(MVFrequency& world, Double pixel) const;
  Bool toPixel(Double& pixel, const MVFrequency& 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 (True for fail, False for success)
  // is the length of the number of conversions and
  // holds an error status for each conversion.
  // <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>
 
  // Set the state that is used for conversions from pixel and frequency to velocity
  // or wavelength. The SpectralCoordinate is constructed  with
  // <src>MDoppler::RADIO</src> and <src>km/s</src> as the velocity conversion state
  // and <src>mm</src> as the wavelength conversion state.
  // The functions in this class which use this state are those that convert
  // to or from velocity.  Also, function <src>format</src> uses the Doppler
  // state set here.  If the function returns False it means the unit was
  // not valid.  There will be an error message in function <src>errorMessage</src>
  // <group>
  Bool setVelocity(const String& velUnit = String("km/s"),
                   MDoppler::Types velType = MDoppler::RADIO);
 
  MDoppler::Types velocityDoppler() const { return velType_p; };
  String velocityUnit() const { return velUnit_p; };
  //
  Bool setWavelengthUnit(const String& waveUnit = String("mm"));
  String wavelengthUnit() const { return waveUnit_p; };
  //
  Bool setNativeType(const SpectralCoordinate::SpecType spcType);
  SpectralCoordinate::SpecType nativeType() const { return nativeType_p; }
 
  // </group>
  // Functions to convert to velocity (uses the current active
  // rest frequency) or wavelength.  There is no reference frame
  // change but you can specify the velocity Doppler and the output
  // units of the velocity with function <src>setVelocity</src>
  // or <src>setWavelength</src> respectively. When the input is a frequency stored
  // as a Double it must be in the current units of the SpectralCoordinate.
  //
  // Note that the extra conversion layer (see function <src>setReferenceConversion</src>)
  // is active in the <src>pixelToVelocity</src> functions (because internally
  // the use <src>toWorld</src>) but not in the <src>frequencyToVelocity</src>
  // or <src>frequencyToWavelength</src> functions.
  // <group>
  Bool pixelToVelocity(Quantum<Double>& velocity, Double pixel) const;
  Bool pixelToVelocity(Double& velocity, Double pixel) const;
  Bool pixelToVelocity(Vector<Double>& velocity, const Vector<Double>& pixel) const;
  //
  Bool frequencyToVelocity(Quantum<Double>& velocity, Double frequency) const;
  Bool frequencyToVelocity(Quantum<Double>& velocity, const MFrequency& frequency) const;
  Bool frequencyToVelocity(Quantum<Double>& velocity, const MVFrequency& frequency) const;
  Bool frequencyToVelocity(Double& velocity, Double frequency) const;
  Bool frequencyToVelocity(Vector<Double>& velocity, const Vector<Double>& frequency) const;
  //
  Bool frequencyToWavelength(Vector<Double>& wavelength, const Vector<Double>& frequency) const;
  Bool frequencyToAirWavelength(Vector<Double>& wavelength, const Vector<Double>& frequency) const;
  // The refractive index of air (argument can be wavelength or airwavelength)
  // according to Greisen et al., 2006, A&A, 464, 746.
  // If airwavelength is used there is an error of the order of 1E-9.
  // Argument must be in micrometers!
  // static Double refractiveIndex(const Double& lambda_um);
  // </group>
 
  // Functions to convert from velocity (uses the current active
  // rest frequency) or wavelength.  There is no reference frame
  // change but you can specify the velocity Doppler and the output
  // units of the velocity with function <src>setVelocity</src>
  // and those of the wavelength with <src>setWavelength</src>.
  // When the input is a frequency stored
  // as a Double it must be  in the current units of the SpectralCoordinate.
  //
  // Note that the extra conversion layer (see function <src>setReferenceConversion</src>)
  // is active in the <src>pixelToVelocity</src> functions (because internally
  // the use <src>toPixel</src>) but not in the <src>frequencyToVelocity</src> functions.
  // <group>
  Bool velocityToPixel(Double& pixel, Double velocity) const;
  Bool velocityToPixel(Vector<Double>& pixel, const Vector<Double>& velocity) const;
  //
  Bool velocityToFrequency(Double& frequency, Double velocity) const;
  Bool velocityToFrequency(Vector<Double>& frequency, const Vector<Double>& velocity) const;
  //
  Bool wavelengthToFrequency(Vector<Double>& frequency, const Vector<Double>& wavelength) const;
  Bool airWavelengthToFrequency(Vector<Double>& frequency, const Vector<Double>& wavelength) const;
  // </group>
 
  // The SpectralCoordinate can maintain a list of rest frequencies
  // (e.g. multiple lines within a band).  However, only
  // one of them is active (e.g. for velocity conversions) at any
  // one time.  Function <src>restFrequency</src> returns that
  // frequency.    Function <src>restFrequencies</src> returns
  // all of the possible restfrequencies.
  //
  // When you construct the SpectralCoordinate, you give it one rest frequency
  // and it is the active one.  Thereafter you can add a new restfrequency
  // with function <src>setRestFrequency</src> (<src>append=True</src>) and
  // that frequency will become the active one.    With this function
  // and <src>append=False</src>, the current active restfrequency will
  // be replaced by the one you give.
  //
  // You can change the list of
  // restfrequencies with function <src>setRestFrequencies</src>. When
  // you do so, you can either replace the list of rest frequencies or append to it.
  // You specify which frequency of the new (appended) list
  // becomes active.
  //
  // You can also select the active rest frequency either by an index into
  // the current list (exception if out of range) given by
  // <src>restFrequencies()</src> or by the value in the list
  // nearest to the frequency you give.
  //
  // Whenever you change the active rest frequency, the class internals
  // are adjusted (e.g. the velocity machine is updated).
  // <group>
  Double restFrequency() const;
  const Vector<Double>& restFrequencies() const;
  Bool setRestFrequency(Double newFrequency, Bool append = False);
  void setRestFrequencies(const Vector<Double>& newFrequencies, uInt which = 0,
                          Bool append = False);
  void selectRestFrequency(uInt which);
  void selectRestFrequency(Double frequency);
  String formatRestFrequencies() const;
  // </group>
 
  // Retrieve/set the frequency system.  Note that setting the
  // frequency system just changes the internal value of the
  // frequency system.  In addition, it will reset the internal
  // conversion frequency system to the new type and delete any
  // conversion machines.
  // <group>
  MFrequency::Types frequencySystem(Bool showConversion = False) const;
  void setFrequencySystem(MFrequency::Types type, Bool verbose = True);
 
  // Transform the SpectralCoordinate to a different native reference frame
  // keeping the conversion layer as is
  Bool transformFrequencySystem(MFrequency::Types type, const MEpoch& epoch,
                                const MPosition& position, const MDirection& direction);
  // </group>
 
  // Report the value of the requested attribute.
  // <group>
  virtual Vector<String> worldAxisNames() const;
  virtual Vector<Double> referencePixel() const;
  virtual Matrix<Double> linearTransform() const;
  virtual Vector<Double> increment() const;
  virtual Vector<Double> referenceValue() 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>
 
  // Get the table, i.e. the pixel and world values. The length of these
  // Vectors will be zero if this axis is pure linear (i.e. if the
  // channel and frequencies are related through an increment and offset).
  // <group>
  Vector<Double> pixelValues() const;
  Vector<Double> worldValues() const;
  // </group>
 
  // Set/get the unit. Adjust the increment and
  // reference value by the ratio of the old and new units.
  // The unit must be compatible with  frequency.
  //<group>
  virtual Bool setWorldAxisUnits(const Vector<String>& units);
  virtual Vector<String> worldAxisUnits() const;
  //</group>
 
  // 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,
  // <src>errorMessage()</src> contains 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>
 
  // 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.  Cannot transform tabular
  // coordinates.  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;
 
  // Format a SpectralCoordinate coordinate world value nicely through the
  // common format interface.  See <linkto class=Coordinate>Coordinate</linkto>
  // for basics.
  //
  // Format types SCIENTIFIC, FIXED, MIXED and DEFAULT are supported.
  // DEFAULT will use MIXED.
  //
  // The world value must always be given in native frequency units.
  // Use argument <src>unit</src> to determine what it will be
  // converted to for formatting. If <src>unit</src> is given, it
  // must be dimensionally consistent with Hz, m, or m/s.
  // If you give a unit consistent with m/s then the
  // appropriate velocity Doppler type is taken from that set by
  // function <src>setVelocity</src>.  There is no frame conversion.
  // If <src>unit</src> is empty, the unit given by <src>setFormatUnit</src>
  // is used.  If this is turn empty, then native units are used.
  virtual String format(String& unit, Coordinate::formatType format, Double worldValue,
                        uInt worldAxis, Bool isAbsolute = True, Bool showAsAbsolute = True,
                        Int precision = -1, Bool usePrecForFixed = False) const;
 
  // Set the default formatter unit (which is initialized to empty).  Must
  // be consistent with Hz or km/s.
  // If the given unit is illegal, False is returned and the internal state unchanged.
  // This unit is used by the function <src>format</src> when the given
  // unit is empty.
  // <group>
  String formatUnit() const { return formatUnit_p; }
  Bool setFormatUnit(const String& unit);
  // </group>
 
  // Convert to FITS header record.  When writing the FITS record,
  // the fields "ctype, crval, crpix", and "cdelt" must already be created. Other header
  // words are created as needed.  Use <src>oneRelative=True</src> to
  // convert zero-relative SpectralCoordinate pixel coordinates to
  // one-relative FITS coordinates, and vice-versa.  If <src>preferVelocity=True</src>
  // the primary axis type will be velocity, if <src>preferWavelength=True</src> it will
  // be wavelength, else frequency. For a velocity axis, if <src>opticalVelDef=False</src>,
  // the radio velocity definition will be used, else optical definition. Similarly for a
  // wavelength axis, if <src>airWaveDef=True</src> air wavelength will be used, the
  // default is vacuum wavelength.
  //<group>
  void toFITS(RecordInterface& header, uInt whichAxis, LogIO& logger, Bool oneRelative = True,
              Bool preferVelocity = True, Bool opticalVelDef = True, Bool preferWavelength = False,
              Bool airWaveDef = False) const;
 
  // Old interface.  Handled by wcs in new interface in FITSCoordinateUtil.cc
  //    static Bool fromFITSOld(SpectralCoordinate &out, String &error,
  //             const RecordInterface &header,
  //             uInt whichAxis,
  //             LogIO &logger, Bool oneRelative=True);
  //</group>
 
  // Save the SpectralCoordinate 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 SpectralCoordinate from a record.
  // A null pointer means that the restoration did not succeed.
  static SpectralCoordinate* restore(const RecordInterface& container, const String& fieldName);
 
  // Convert from String to spectral type and vice versa.
  static Bool specTypetoString(String& stypeString, const SpecType& specType);
  static Bool stringtoSpecType(SpecType& specType, const String& stypeString);
 
  // Make a copy of the SpectralCoordinate using new. The caller is responsible for calling
  // delete.
  virtual Coordinate* clone() const;
 
  ostream& print(ostream& os) const;
 
  // is this a tabular coordinate?
  Bool isTabular() const;
 
 private:
  std::unique_ptr<TabularCoordinate> _tabular;  // Tabular coordinate OR
  mutable ::wcsprm wcs_p;                       // wcs structure is used
  Double to_hz_p;                               // Convert from current world units to Hz
  Double to_m_p;                                // Convert from current wavelength units to m
                                                //
  MFrequency::Types type_p, conversionType_p;   // Frequency system and conversion system
  Vector<Double> restfreqs_p;                   // List of possible rest frequencies
  uInt restfreqIdx_p;                           // Current active rest frequency index
 
  // Conversion machines; for pixel<->world conversions only.
  mutable MFrequency::Convert* pConversionMachineTo_p;    // For type_p -> conversionType_p
  mutable MFrequency::Convert* pConversionMachineFrom_p;  // For conversionType_p -> type_p
 
  VelocityMachine*
      pVelocityMachine_p;     // The velocity machine does all conversions between world & velocity.
  MDoppler::Types velType_p;  // Velocity Doppler
  String velUnit_p;           // Velocity unit
                              //
  String waveUnit_p;          // Wavelength unit for conversions between world & wavelength
  SpectralCoordinate::SpecType nativeType_p;  // The native spectral type
                                              //
  Unit unit_p;                                // World axis unit
  String axisName_p;                          // The axis name
  String formatUnit_p;                        // The default unit for the format function
                                              //
  MDirection direction_p;                     // These are a part of the frame set for
  MPosition position_p;                       // the reference conversions machines
  MEpoch epoch_p;                             // They are only private so we can save their state
 
  // Format checker
  void checkFormat(Coordinate::formatType& format, const Bool) const;
 
  // Copy private data
  void copy(const SpectralCoordinate& other);
 
  // Convert to and from conversion reference type
  virtual void convertTo(Vector<Double>& world) const;
  virtual void convertFrom(Vector<Double>& world) const;
 
  // Deletes and sets pointer to 0
  void deleteVelocityMachine();
 
  // Deletes and sets pointers to 0
  void deleteConversionMachines();
 
  // Set up pixel<->world conversion machines
  // Returns: 3 (machines were noOPs, machines deleted)
  //          2 (types the same, machines deleted),
  //          1 (machines created and functioning)
  //         -1 (machines could not make trial conversion, machines deleted)
  Int makeConversionMachines(MFrequency::Types type, MFrequency::Types conversionType,
                             const MEpoch& epoch, const MPosition& position,
                             const MDirection& direction);
 
  // Create velocity<->frequency machine
  void makeVelocityMachine(const String& velUnit, MDoppler::Types velType, const Unit& freqUnit,
                           MFrequency::Types freqType, Double restFreq);
 
  // Make spectral wcs structure (items in Hz)
  static void makeWCS(wcsprm& wcs, const String& ctype, Double refPix, Double refVal, Double inc,
                      Double pc, Double restFreq);
 
  // Record restoration handling
  // <group>
  static SpectralCoordinate* restoreVersion1(const RecordInterface& container);
  static SpectralCoordinate* restoreVersion2(const RecordInterface& container);
  static void restoreVelocity(SpectralCoordinate*& pSpectral, const RecordInterface& subrec);
  static void restoreRestFrequencies(SpectralCoordinate*& pSpectral, const RecordInterface& subrec,
                                     Double restFreq);
  static void restoreConversion(SpectralCoordinate*& pSpectral, const RecordInterface& subrec);
 
  // </group>
 
  // Interconvert between the current units and wcs units (Hz)
  // <group>
  void toCurrent(Vector<Double>& value) const;
  void fromCurrent(Vector<Double>& value) const;
  // </group>
 
  // Return unit conversion vector for converting to current units
  const Vector<Double> toCurrentFactors() const;
 
  // Update Velocity Machine
  void updateVelocityMachine(const String& velUnit, MDoppler::Types velType);
  // Restore wcs stuff from Record
  static Bool wcsRestore(Double& crval, Double& crpix, Double& cdelt, Double& pc, String& ctype,
                         const RecordInterface& rec);
 
  // Save wcs stuff into Record
  Bool wcsSave(RecordInterface& rec, const wcsprm& wcs, const String& fieldName) const;
 
  void _setTabulatedFrequencies(const Vector<Double>& freqs);
};
 
ostream& operator<<(ostream& os, const SpectralCoordinate& spcoord);
 
}  // namespace casacore
 
#endif