casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SpectralCoordinate.h
Go to the documentation of this file.
00001 //# SpectralCoordinate.h: Interconvert between pixel and frequency.
00002 //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003,2004
00003 //# Associated Universities, Inc. Washington DC, USA.
00004 //#
00005 //# This library is free software; you can redistribute it and/or modify it
00006 //# under the terms of the GNU Library General Public License as published by
00007 //# the Free Software Foundation; either version 2 of the License, or (at your
00008 //# option) any later version.
00009 //#
00010 //# This library is distributed in the hope that it will be useful, but WITHOUT
00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00012 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00013 //# License for more details.
00014 //#
00015 //# You should have received a copy of the GNU Library General Public License
00016 //# along with this library; if not, write to the Free Software Foundation,
00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00018 //#
00019 //# Correspondence concerning AIPS++ should be addressed as follows:
00020 //#        Internet email: aips2-request@nrao.edu.
00021 //#        Postal address: AIPS++ Project Office
00022 //#                        National Radio Astronomy Observatory
00023 //#                        520 Edgemont Road
00024 //#                        Charlottesville, VA 22903-2475 USA
00025 //#
00026 //#
00027 //# $Id: SpectralCoordinate.h 20491 2009-01-16 08:33:56Z gervandiepen $
00028 
00029 
00030 #ifndef COORDINATES_SPECTRALCOORDINATE_H
00031 #define COORDINATES_SPECTRALCOORDINATE_H
00032 
00033 #include <casa/aips.h>
00034 #include <casa/Arrays/Vector.h>
00035 #include <coordinates/Coordinates/Coordinate.h>
00036 #include <coordinates/Coordinates/ObsInfo.h>
00037 #include <measures/Measures/MFrequency.h>
00038 #include <measures/Measures/MDoppler.h>
00039 #include <measures/Measures/MDirection.h>
00040 #include <measures/Measures/MPosition.h>
00041 #include <measures/Measures/MEpoch.h>
00042 #include <casa/Quanta/Quantum.h>
00043 
00044 #include <wcslib/wcs.h>
00045 
00046 #include <memory>
00047 
00048 namespace casa { //# NAMESPACE CASA - BEGIN
00049 
00050 
00051 class TabularCoordinate;
00052 class LogIO;
00053 class MVFrequency;
00054 class VelocityMachine;
00055 template<class T> class Quantum;
00056 
00057 
00058 // <summary>
00059 // Interconvert pixel and frequency values.
00060 // </summary>
00061 
00062 // <use visibility=export>
00063 
00064 // <reviewed reviewer="Peter Barnes" date="1999/12/24" tests="tSpectralCoordinate">
00065 // </reviewed>
00066 //
00067 // <prerequisite>
00068 //   <li> <linkto class=Coordinate>Coordinate</linkto>
00069 //   <li> <linkto class=MFrequency>MFrequency</linkto>,
00070 //        <linkto class=MDoppler>MDoppler</linkto> and
00071 //        <linkto class=VelocityMachine>VelocityMachine</linkto> 
00072 //        classes if you want radial velocities.
00073 // </prerequisite>
00074 //
00075 // <synopsis>
00076 // This class performs the mapping from pixel to frequency. 
00077 // This can be done via a Tabular lookup or via an algorithmic
00078 // implementation which may be linear or non-linear.  The latter
00079 // is implemented via the WCS library.
00080 // 
00081 // </synopsis>
00082 //
00083 
00084 // <note role=caution>
00085 // All pixels coordinates are zero relative.
00086 // </note>
00087 //
00088 // <example>
00089 // Let us make a linear SpectralCoordinate first
00090 // <srcblock>
00091 //   Double restfreq = 1.420405752E9;
00092 //   Double crpix = 10.0;
00093 //   Double crval = 1.4e9;
00094 //   Double cdelt = 1.0e6;
00095 //   SpectralCoordinate sc(MFrequency::TOPO, crval, cdelt, crpix, restfreq);
00096 //
00097 //   Double world, pixel;
00098 //   pixel = 12.1;
00099 //   if (!sc.toWorld(world, pixel)) {
00100 //      cerr << "Error : " << sc.errorMessage() << endl;
00101 //   } else {
00102 //      cerr << "pixel, world = " << pixel << ", " << world << endl;
00103 //   }
00104 //
00105 // </srcblock>
00106 // </example>
00107 //
00108 // <example>
00109 // Now we make a non-linear SpectralCoordinate 
00110 // <srcblock>
00111 //   Vector<Double> freqs(5);
00112 //   freqs(0) = 1.4e9; freqs(1) = 1.41e9;
00113 //   freqs(2) = 1.43e9; freqs(3) = 1.44e9;
00114 //   freqs(4) = 1.47e9;
00115 //   SpectralCoordinate sc(MFrequency::LSRK, freqs, restfreq);
00116 //
00117 //   Double world, pixel;
00118 //   world = 1.42e9;
00119 //   if (!sc.toPixel(pixel, world)) {
00120 //      cerr << "Error : " << sc.errorMessage() << endl;
00121 //   } else {
00122 //      cerr << "world, pixel = " << world << ", " << pixel << endl;
00123 //   }
00124 //
00125 // </srcblock>
00126 // </example>
00127 //
00128 // <motivation>
00129 //  Spectral-line astronomy requires a specialized SpectralCoordinate.
00130 // </motivation>
00131 
00132 // <todo asof="2000/01/01">
00133 //   <li> Allow other than linear interpolations for frequency lookup.
00134 // </todo>
00135 // 
00136 
00137 class SpectralCoordinate : public Coordinate
00138 {
00139 public:
00140          enum SpecType { // taken from the FITS spectral coordinate type codes
00141                  FREQ,
00142                  VRAD,
00143                  VOPT,
00144                  BETA,
00145                  WAVE,
00146                  AWAV
00147          };
00148 
00149          // Default constructor.    It is equivalent to doing
00150     // SpectralCoordinate(MFrequency::TOPO, 0.0, 1.0, 0.0)
00151     SpectralCoordinate();
00152 
00153     // Create a linear frequency axis SpectralCoordinate
00154     // <src>f0</src> is the frequency of the reference pixel, <src>inc</src> is the pixel increment,
00155     // <src>refPix</src> is the reference pixel. You can
00156     // optionally store the rest frequency for later use in calculating radial
00157     // velocities.  Use 0 for restFrequency if continuum.
00158     //
00159     // Frequencies and increments initially in Hz.
00160     SpectralCoordinate(MFrequency::Types type, Double f0, Double inc, 
00161                        Double refPix, Double restFrequency = 0.0);
00162 
00163     // Create linear frequency axis SpectralCoordinate with Quantum-based interface. 
00164     // Parameters are the same as above.
00165     // Regardless of the units of the Quanta, the initial units
00166     // of the SpectralCoordinate will be Hz.  You can change it to 
00167     // something else with the setWorldAxisUnits method later if you want.
00168     // Use 0 for restFrequency if continuum.
00169     SpectralCoordinate(MFrequency::Types type, const Quantum<Double>& f0, 
00170                        const Quantum<Double>& inc, Double refPix, 
00171                        const Quantum<Double>& restFrequency = Quantum<Double>(0.0,"Hz"));
00172 
00173     // Construct a SpectralCoordinate with the specified frequencies (in Hz).
00174     // This axis can be nonlinear; the increments and related 
00175     // functions return the <src>average</src> values
00176     // (calculated from the first and last pixel's frequencies).
00177     //
00178     // A linear interpolation/extrapolation is used for pixels which are
00179     // not supplied. The reference pixel is chosen to be 0.
00180     // The frequencies must increase or decrease monotonically (otherwise
00181     // the toPixel lookup would not be possible).
00182     // Use 0 for restFrequency if continuum.
00183     SpectralCoordinate(MFrequency::Types type, const Vector<Double> &freqs,
00184                        Double restFrequency = 0.0);
00185     
00186     // Construct a SpectralCoordinate with the specified frequencies
00187     // with Quantum-based interface. 
00188     // Parameters are the same as above.
00189     // Regardless of the units of the Quanta, the initial units
00190     // of the SpectralCoordinate will be Hz.
00191     // Use 0 for restFrequency if continuum.
00192     SpectralCoordinate(MFrequency::Types type, const Quantum<Vector<Double> >& freqs,
00193                        const Quantum<Double>& restFrequency = Quantum<Double>(0.0,"Hz"));
00194 
00195     // Construct a SpectralCoordinate with the specified velocities (in km/s).
00196     // They will be converted to Hz and the SpectralCoordinate constructed.
00197     // This axis can be nonlinear; the increments and related 
00198     // functions return the <src>average</src> values
00199     // (calculated from the first and last pixel's frequencies).
00200     //
00201     // A linear interpolation/extrapolation is used for pixels which are
00202     // not supplied. The reference pixel is chosen to be 0.
00203     // The velocities must increase or decrease monotonically (otherwise
00204     // the toPixel lookup would not be possible).
00205     SpectralCoordinate(MFrequency::Types freqType, MDoppler::Types velType, 
00206                        const Vector<Double>& velocities,  const String& velUnit,
00207                        Double restFrequency = 0.0);
00208 
00209     // Construct a SpectralCoordinate with the specified wavelengths (in mm).
00210     // They will be converted to Hz and the SpectralCoordinate constructed.
00211     // This axis can be nonlinear; the increments and related 
00212     // functions return the <src>average</src> values
00213     // (calculated from the first and last pixel's frequencies).
00214     // If inAir is True, the input wavelengths are assumed to be Air Wavelengths.
00215     // They are converted to vacuum frequency using the refractive index
00216     // which is calculated based on the mean input air wavelength.
00217     //
00218     // A linear interpolation/extrapolation is used for pixels which are
00219     // not supplied. The reference pixel is chosen to be 0.
00220     // The wavelengths must increase or decrease monotonically (otherwise
00221     // the toPixel lookup would not be possible).
00222     SpectralCoordinate(MFrequency::Types freqType,     
00223                        const Vector<Double>& wavelengths,  const String& waveUnit,
00224                        Double restFrequency = 0.0, Bool inAir = False);
00225 
00226     // Construct from wcs structure.  Must hold only a spectral wcs structure
00227     // Specify whether the absolute pixel coordinates in the wcs structure
00228     // are 0- or 1-relative.  The coordinate is always constructed with 0-relative
00229     // pixel coordinates
00230     SpectralCoordinate(MFrequency::Types freqType, const ::wcsprm& wcs, Bool oneRel=True);
00231     
00232     // Copy constructor (copy semantics).
00233     SpectralCoordinate(const SpectralCoordinate &other);
00234 
00235     // Assignment (copy semantics).
00236     SpectralCoordinate &operator=(const SpectralCoordinate &other);
00237 
00238     // Destructor.
00239     virtual ~SpectralCoordinate();
00240 
00241     // Always returns Coordinate::SPECTRAL.
00242     virtual Coordinate::Type type() const;
00243 
00244     // Always returns the String "Spectral".
00245     virtual String showType() const;
00246 
00247     // Always returns 1.
00248     // <group>
00249     virtual uInt nPixelAxes() const;
00250     virtual uInt nWorldAxes() const;
00251     // </group>
00252 
00253     // Set extra conversion layer.  Whenever a conversion from pixel to world is done,
00254     // the world value is then further converted to this MFrequency::Types value.
00255     // For example, your SpectralCoordinate may be defined in LSRK.
00256     // You can use this to get the world values out in say BARY. You must
00257     // specify the position on earth, the epoch and the direction for the conversions
00258     // and it is your responsibility to ensure they are viable.
00259     // Similarly, whenever you convert from world to pixel, the world
00260     // value is assumed to be that appropriate to the setReferenceConversion type.
00261     // It is first converted to the MFrequency::Types with which the
00262     // SpectralCoordinate was constructed and from there to pixel.
00263     // If you don't call this function, or you set the same type
00264     // for which the SpectralCoordinate was constructed, no extra
00265     // conversions occur.   Some conversions will fail.  These are the
00266     // ones that require extra frame information (radial velocity) such
00267     // as to REST. This will be added later.  In this case this function
00268     // returns False (and the conversion parameters are all left as they were),
00269    //  else it returns True.
00270     // <group>
00271     Bool setReferenceConversion (MFrequency::Types type,
00272                                  const MEpoch& epoch, const MPosition& position,
00273                                  const MDirection& direction);
00274     void getReferenceConversion (MFrequency::Types& type,
00275                                  MEpoch& epoch, MPosition& position,
00276                                  MDirection& direction) const
00277        {type = conversionType_p; epoch=epoch_p; 
00278         position=position_p; direction=direction_p;};
00279     // </group>
00280 
00281     // Convert a pixel to a world coordinate or vice versa. Returns True
00282     // if the conversion succeeds, otherwise it returns False and
00283     // <src>errorMessage()</src> contains an error message.  The input vectors
00284     // must be of length one and the output vectors are resized if they are not
00285     // already of length one.
00286     // <group>
00287     virtual Bool toWorld(Vector<Double> &world, 
00288                          const Vector<Double> &pixel) const;
00289     virtual Bool toPixel(Vector<Double> &pixel, 
00290                          const Vector<Double> &world) const;
00291     Bool toWorld(Double& world, const Double& pixel) const;
00292     Bool toPixel(Double& pixel, const Double& world) const;
00293     // </group>
00294 
00295     // Convert a pixel (channel number) into an MFrequency or MVFrequency and vice 
00296     // versa. Usually you will do
00297     // this for calculating velocities or converting frequencies from one frame
00298     // to another.
00299     // <group>
00300     Bool toWorld(MFrequency &world,
00301                  Double pixel) const;
00302     Bool toPixel(Double& pixel, const MFrequency &world) const;
00303     Bool toWorld(MVFrequency &world,
00304                  Double pixel) const;
00305     Bool toPixel(Double& pixel, const MVFrequency &world) const;
00306     // </group>
00307 
00308     // Batch up a lot of transformations. The first (most rapidly varying) axis
00309     // of the matrices contain the coordinates. Returns False if any conversion
00310     // failed  and  <src>errorMessage()</src> will hold a message.
00311     // The <src>failures</src> array (True for fail, False for success)
00312     // is the length of the number of conversions and
00313     // holds an error status for each conversion.  
00314     // <group>
00315     virtual Bool toWorldMany(Matrix<Double>& world,
00316                               const Matrix<Double>& pixel,
00317                               Vector<Bool>& failures) const;
00318     virtual Bool toPixelMany(Matrix<Double>& pixel,
00319                              const Matrix<Double>& world,   
00320                              Vector<Bool>& failures) const;
00321     // </group>
00322 
00323     // Set the state that is used for conversions from pixel and frequency to velocity
00324     // or wavelength. The SpectralCoordinate is constructed  with 
00325     // <src>MDoppler::RADIO</src> and <src>km/s</src> as the velocity conversion state
00326     // and <src>mm</src> as the wavelength conversion state.
00327     // The functions in this class which use this state are those that convert
00328     // to or from velocity.  Also, function <src>format</src> uses the Doppler
00329     // state set here.  If the function returns False it means the unit was 
00330     // not valid.  There will be an error message in function <src>errorMessage</src>
00331     // <group>
00332     Bool setVelocity (const String& velUnit=String("km/s"),
00333                       MDoppler::Types velType=MDoppler::RADIO);
00334 
00335     MDoppler::Types velocityDoppler () const {return velType_p;};
00336     String velocityUnit () const {return velUnit_p;};
00337     //
00338     Bool setWavelengthUnit (const String& waveUnit=String("mm"));
00339     String wavelengthUnit () const {return waveUnit_p;};
00340     //
00341     Bool setNativeType (const SpectralCoordinate::SpecType spcType);
00342     SpectralCoordinate::SpecType nativeType() const {return nativeType_p;}
00343 
00344     // </group>
00345     // Functions to convert to velocity (uses the current active
00346     // rest frequency) or wavelength.  There is no reference frame
00347     // change but you can specify the velocity Doppler and the output
00348     // units of the velocity with function <src>setVelocity</src>
00349     // or <src>setWavelength</src> respectively. When the input is a frequency stored 
00350     // as a Double it must be in the current units of the SpectralCoordinate.  
00351     // 
00352     // Note that the extra conversion layer (see function <src>setReferenceConversion</src>)
00353     // is active in the <src>pixelToVelocity</src> functions (because internally
00354     // the use <src>toWorld</src>) but not in the <src>frequencyToVelocity</src> 
00355     // or <src>frequencyToWavelength</src> functions.
00356     // <group>  
00357     Bool pixelToVelocity (Quantum<Double>& velocity, Double pixel) const;
00358     Bool pixelToVelocity (Double& velocity, Double pixel) const;
00359     Bool pixelToVelocity (Vector<Double>& velocity, const Vector<Double>& pixel) const;
00360     //
00361     Bool frequencyToVelocity (Quantum<Double>& velocity, Double frequency) const;
00362     Bool frequencyToVelocity (Quantum<Double>& velocity, const MFrequency& frequency) const;
00363     Bool frequencyToVelocity (Quantum<Double>& velocity, const MVFrequency& frequency) const;
00364     Bool frequencyToVelocity (Double& velocity, Double frequency) const;
00365     Bool frequencyToVelocity (Vector<Double>& velocity, const Vector<Double>& frequency) const;
00366     //
00367     Bool frequencyToWavelength (Vector<Double>& wavelength, const Vector<Double>& frequency) const;
00368     Bool frequencyToAirWavelength (Vector<Double>& wavelength, const Vector<Double>& frequency) const;
00369     // The refractive index of air (argument can be wavelength or airwavelength)
00370     // according to Greisen et al., 2006, A&A, 464, 746.
00371     // If airwavelength is used there is an error of the order of 1E-9.
00372     // Argument must be in micrometers!  
00373     //static Double refractiveIndex(const Double& lambda_um);
00374     // </group>
00375 
00376     // Functions to convert from velocity (uses the current active
00377     // rest frequency) or wavelength.  There is no reference frame
00378     // change but you can specify the velocity Doppler and the output
00379     // units of the velocity with function <src>setVelocity</src>
00380     // and those of the wavelength with <src>setWavelength</src>. 
00381     // When the input is a frequency stored 
00382     // as a Double it must be  in the current units of the SpectralCoordinate.  
00383     //
00384     // Note that the extra conversion layer (see function <src>setReferenceConversion</src>)
00385     // is active in the <src>pixelToVelocity</src> functions (because internally
00386     // the use <src>toPixel</src>) but not in the <src>frequencyToVelocity</src> functions.
00387     // <group>  
00388     Bool velocityToPixel (Double& pixel, Double velocity) const;
00389     Bool velocityToPixel (Vector<Double>& pixel, const Vector<Double>& velocity) const;
00390     // 
00391     Bool velocityToFrequency (Double& frequency, Double velocity) const;
00392     Bool velocityToFrequency (Vector<Double>& frequency, const Vector<Double>& velocity) const;
00393     //
00394     Bool wavelengthToFrequency (Vector<Double>& frequency, const Vector<Double>& wavelength) const;
00395     Bool airWavelengthToFrequency (Vector<Double>& frequency, const Vector<Double>& wavelength) const;
00396     // </group>
00397 
00398     // The SpectralCoordinate can maintain a list of rest frequencies
00399     // (e.g. multiple lines within a band).  However, only
00400     // one of them is active (e.g. for velocity conversions) at any 
00401     // one time.  Function <src>restFrequency</src> returns that
00402     // frequency.    Function <src>restFrequencies</src> returns   
00403     // all of the possible restfrequencies.
00404     //
00405     // When you construct the SpectralCoordinate, you give it one rest frequency
00406     // and it is the active one.  Thereafter you can add a new restfrequency
00407     // with function <src>setRestFrequency</src> (<src>append=True</src>) and 
00408     // that frequency will become the active one.    With this function
00409     // and <src>append=False</src>, the current active restfrequency will
00410     // be replaced by the one you give.
00411     //
00412     // You can change the list of
00413     // restfrequencies with function <src>setRestFrequencies</src>. When
00414     // you do so, you can either replace the list of rest frequencies or append to it.
00415     // You specify which frequency of the new (appended) list
00416     // becomes active.  
00417     //
00418     // You can also select the active rest frequency either by an index into
00419     // the current list (exception if out of range) given by
00420     // <src>restFrequencies()</src> or by the value in the list
00421     // nearest to the frequency you give.
00422     //
00423     // Whenever you change the active rest frequency, the class internals
00424     // are adjusted (e.g. the velocity machine is updated).
00425     // <group>
00426     Double restFrequency() const;
00427     const Vector<Double>& restFrequencies() const;
00428     Bool setRestFrequency(Double newFrequency, Bool append=False);
00429     void setRestFrequencies(const Vector<Double>& newFrequencies, uInt which=0,
00430                             Bool append=False);
00431     void selectRestFrequency(uInt which);
00432     void selectRestFrequency(Double frequency);
00433     String formatRestFrequencies () const;
00434     // </group>
00435   
00436     // Retrieve/set the frequency system.  Note that setting the
00437     // frequency system just changes the internal value of the
00438     // frequency system.  In addition, it will reset the internal
00439     // conversion frequency system to the new type and delete any
00440     // conversion machines.  
00441     // <group>
00442     MFrequency::Types frequencySystem(Bool showConversion=False) const;
00443     void  setFrequencySystem(MFrequency::Types type);
00444     // </group>
00445 
00446     // Report the value of the requested attribute.
00447     // <group>
00448     virtual Vector<String> worldAxisNames() const;
00449     virtual Vector<Double> referencePixel() const;
00450     virtual Matrix<Double> linearTransform() const;
00451     virtual Vector<Double> increment() const;
00452     virtual Vector<Double> referenceValue() const;
00453     // </group>
00454 
00455     // Set the value of the requested attribute. Note that these just
00456     // change the internal values, they do not cause any recomputation.
00457     // <group>
00458     virtual Bool setWorldAxisNames(const Vector<String> &names);
00459     virtual Bool setReferencePixel(const Vector<Double> &refPix);
00460     virtual Bool setLinearTransform(const Matrix<Double> &xform);
00461     virtual Bool setIncrement(const Vector<Double> &inc) ;
00462     virtual Bool setReferenceValue(const Vector<Double> &refval);
00463     // </group>
00464 
00465     // Get the table, i.e. the pixel and world values. The length of these
00466     // Vectors will be zero if this axis is pure linear (i.e. if the
00467     // channel and frequencies are related through an increment and offset).
00468     // <group>
00469     Vector<Double> pixelValues() const;
00470     Vector<Double> worldValues() const;
00471     // </group>
00472 
00473     // Set/get the unit. Adjust the increment and
00474     // reference value by the ratio of the old and new units.
00475     // The unit must be compatible with  frequency.
00476     //<group>
00477     virtual Bool setWorldAxisUnits(const Vector<String> &units);
00478     virtual Vector<String> worldAxisUnits() const;
00479     //</group>
00480 
00481     // Comparison function. Any private Double data members are compared
00482     // with the specified fractional tolerance.  Don't compare on the specified 
00483     // axes in the Coordinate.  If the comparison returns False, 
00484     // <src>errorMessage()</src> contains a message about why.
00485     // <group>
00486     virtual Bool near(const Coordinate& other, 
00487                       Double tol=1e-6) const;
00488     virtual Bool near(const Coordinate& other, 
00489                       const Vector<Int>& excludeAxes,
00490                       Double tol=1e-6) const;
00491     // </group>
00492 
00493 
00494     // Find the Coordinate for when we Fourier Transform ourselves.  This pointer
00495     // must be deleted by the caller. Axes specifies which axes of the Coordinate
00496     // you wish to transform.   Shape specifies the shape of the image
00497     // associated with all the axes of the Coordinate.   Currently the
00498     // output reference pixel is always shape/2.  Cannot transform tabular
00499     // coordinates.  If the pointer returned is 0, it failed with a message
00500     // in <src>errorMessage</src>
00501     virtual Coordinate* makeFourierCoordinate (const Vector<Bool>& axes,
00502                                                const Vector<Int>& shape) const;
00503 
00504 
00505     // Format a SpectralCoordinate coordinate world value nicely through the
00506     // common format interface.  See <linkto class=Coordinate>Coordinate</linkto>
00507     // for basics.
00508     //
00509     // Format types SCIENTIFIC, FIXED, MIXED and DEFAULT are supported.
00510     // DEFAULT will use MIXED.
00511     //
00512     // The world value must always be given in native frequency units.
00513     // Use argument <src>unit</src> to determine what it will be 
00514     // converted to for formatting. If <src>unit</src> is given, it 
00515     // must be dimensionally consistent with Hz, m, or m/s.   
00516     // If you give a unit consistent with m/s then the
00517     // appropriate velocity Doppler type is taken from that set by
00518     // function <src>setVelocity</src>.  There is no frame conversion.
00519     // If <src>unit</src> is empty, the unit given by <src>setFormatUnit</src> 
00520     // is used.  If this is turn empty, then native units are used.
00521     virtual String format(String& unit,
00522                           Coordinate::formatType format,
00523                           Double worldValue,  
00524                           uInt worldAxis,
00525                           Bool isAbsolute=True,
00526                           Bool showAsAbsolute=True,
00527                           Int precision=-1, Bool usePrecForFixed=False) const;
00528 
00529     // Set the default formatter unit (which is initialized to empty).  Must 
00530     // be consistent with Hz or km/s.  
00531     // If the given unit is illegal, False is returned and the internal state unchanged.
00532     // This unit is used by the function <src>format</src> when the given
00533     // unit is empty.  
00534     // <group>
00535     String formatUnit () const {return formatUnit_p;}
00536     Bool setFormatUnit (const String& unit);
00537     // </group>
00538 
00539     // Convert to FITS header record.  When writing the FITS record,
00540     // the fields "ctype, crval, crpix", and "cdelt" must already be created. Other header
00541     // words are created as needed.  Use <src>oneRelative=True</src> to
00542     // convert zero-relative SpectralCoordinate pixel coordinates to 
00543     // one-relative FITS coordinates, and vice-versa.  If <src>preferVelocity=True</src>
00544     // the primary axis type will be velocity, if <src>preferWavelength=True</src> it will
00545     // be wavelength, else frequency. For a velocity axis, if <src>opticalVelDef=False</src>,
00546     // the radio velocity definition will be used, else optical definition. Similarly for a
00547     // wavelength axis, if <src>airWaveDef=True</src> air wavelength will be used, the
00548     // default is vacuum wavelength.
00549     //<group>
00550     void toFITS(RecordInterface &header, uInt whichAxis, 
00551                 LogIO &logger, Bool oneRelative=True,
00552                 Bool preferVelocity=True, Bool opticalVelDef=True,
00553                 Bool preferWavelength=False, Bool airWaveDef=False) const;
00554 
00555 // Old interface.  Handled by wcs in new interface in FITSCoordinateUtil.cc
00556 //    static Bool fromFITSOld(SpectralCoordinate &out, String &error,
00557 //                       const RecordInterface &header, 
00558 //                       uInt whichAxis,
00559 //                       LogIO &logger, Bool oneRelative=True);
00560     //</group>
00561 
00562     // Save the SpectralCoordinate into the supplied record using the supplied field name.
00563     // The field must not exist, otherwise <src>False</src> is returned.
00564     virtual Bool save(RecordInterface &container,
00565                     const String &fieldName) const;
00566 
00567     // Recover the SpectralCoordinate from a record.
00568     // A null pointer means that the restoration did not succeed.
00569     static SpectralCoordinate* restore(const RecordInterface &container,
00570                                        const String &fieldName);
00571 
00572     // Convert from String to spectral type and vice versa.
00573     static Bool specTypetoString(String &stypeString, const SpecType &specType);
00574     static Bool stringtoSpecType(SpecType &specType, const String &stypeString);
00575 
00576     // Make a copy of the SpectralCoordinate using new. The caller is responsible for calling
00577     // delete.
00578     virtual Coordinate* clone() const;
00579 
00580     ostream& print(ostream& os) const;
00581 
00582 private:
00583 
00584     std::auto_ptr<TabularCoordinate> _tabular;                     // Tabular coordinate OR
00585     mutable ::wcsprm wcs_p;                              // wcs structure is used 
00586     Double to_hz_p;                                    // Convert from current world units to Hz
00587     Double to_m_p;                                     // Convert from current wavelength units to m
00588 //
00589     MFrequency::Types type_p, conversionType_p;        // Frequency system and conversion system
00590     Vector<Double> restfreqs_p;                        // List of possible rest frequencies
00591     uInt restfreqIdx_p;                                // Current active rest frequency index
00592 
00593                                                                // Conversion machines; for pixel<->world conversions only.
00594     mutable MFrequency::Convert* pConversionMachineTo_p;       // For type_p -> conversionType_p
00595     mutable MFrequency::Convert* pConversionMachineFrom_p;     // For conversionType_p -> type_p
00596 
00597     VelocityMachine* pVelocityMachine_p;           // The velocity machine does all conversions between world & velocity.
00598     MDoppler::Types velType_p;                     // Velocity Doppler
00599     String velUnit_p;                              // Velocity unit
00600 //
00601     String waveUnit_p;                             // Wavelength unit for conversions between world & wavelength
00602     SpectralCoordinate::SpecType nativeType_p;     // The native spectral type
00603 //
00604     Unit unit_p;                                   // World axis unit
00605     String axisName_p;                             // The axis name
00606     String formatUnit_p;                           // The default unit for the format function
00607 // 
00608     MDirection direction_p;                // These are a part of the frame set for
00609     MPosition position_p;                  // the reference conversions machines
00610     MEpoch epoch_p;                        // They are only private so we can save their state
00611 
00612 // Format checker
00613     void checkFormat(Coordinate::formatType& format,
00614                      const Bool ) const;
00615 
00616 // Copy private data
00617    void copy (const SpectralCoordinate& other);
00618 
00619 // Convert to and from conversion reference type
00620     virtual void convertTo (Vector<Double>& world) const;
00621     virtual void convertFrom (Vector<Double>& world) const;
00622 
00623 // Deletes and sets pointer to 0
00624     void deleteVelocityMachine ();
00625 
00626 // Deletes and sets pointers to 0
00627     void deleteConversionMachines ();
00628 
00629 // Set up pixel<->world conversion machines
00630 // Returns: 3 (machines were noOPs, machines deleted)
00631 //          2 (types the same, machines deleted), 
00632 //          1 (machines created and functioning)
00633 //         -1 (machines could not make trial conversion, machines deleted)
00634     Int makeConversionMachines (MFrequency::Types type,  MFrequency::Types conversionType,
00635                                  const MEpoch& epoch, 
00636                                  const MPosition& position, 
00637                                  const MDirection& direction);
00638 
00639 // Create velocity<->frequency machine 
00640     void makeVelocityMachine (const String& velUnit,                  
00641                               MDoppler::Types velType,
00642                               const Unit& freqUnit,
00643                               MFrequency::Types freqType,
00644                               Double restFreq);
00645 
00646 // Make spectral wcs structure (items in Hz)
00647    static void makeWCS(wcsprm& wcs, const String& ctype, Double refPix, Double refVal, 
00648                        Double inc, Double pc, Double restFreq);
00649 
00650 // Record restoration handling
00651 // <group>
00652    static SpectralCoordinate* restoreVersion1 (const RecordInterface& container);
00653    static SpectralCoordinate* restoreVersion2 (const RecordInterface& container);
00654    static void restoreVelocity (SpectralCoordinate*& pSpectral,
00655                                 const RecordInterface& subrec);
00656    static void restoreRestFrequencies (SpectralCoordinate*& pSpectral,
00657                                        const RecordInterface& subrec,
00658                                        Double restFreq);
00659    static void restoreConversion (SpectralCoordinate*& pSpectral,
00660                                   const RecordInterface& subrec);
00661 
00662 // </group>
00663 
00664 // Interconvert between the current units and wcs units (Hz)
00665 // <group>
00666     void toCurrent(Vector<Double>& value) const;
00667     void fromCurrent(Vector<Double>& value) const;
00668 // </group>
00669 
00670 // Return unit conversion vector for converting to current units
00671    const Vector<Double> toCurrentFactors () const;
00672 
00673 // Update Velocity Machine
00674    void updateVelocityMachine (const String& velUnit, 
00675                                MDoppler::Types velType);
00676 // Restore wcs stuff from Record
00677    static Bool wcsRestore (Double& crval, Double& crpix, Double& cdelt,
00678                            Double& pc, String& ctype,
00679                            const RecordInterface& rec);
00680 
00681 // Save wcs stuff into Record
00682    Bool wcsSave (RecordInterface& rec, const wcsprm& wcs,
00683                  const String& fieldName) const;
00684 
00685    void _setTabulatedFrequencies(const Vector<Double>& freqs);
00686 
00687 };
00688 
00689 ostream &operator<<(ostream &os, const SpectralCoordinate& spcoord);
00690 
00691 } //# NAMESPACE CASA - END
00692 
00693 
00694 #endif