diff options
Diffstat (limited to 'ast/specmap.c')
| -rw-r--r-- | ast/specmap.c | 4671 |
1 files changed, 4671 insertions, 0 deletions
diff --git a/ast/specmap.c b/ast/specmap.c new file mode 100644 index 0000000..164179a --- /dev/null +++ b/ast/specmap.c @@ -0,0 +1,4671 @@ +/* +*class++ +* Name: +* SpecMap + +* Purpose: +* Sequence of spectral coordinate conversions. + +* Constructor Function: +c astSpecMap (also see astSpecAdd) +f AST_SPECMAP (also see AST_SPECADD) + +* Description: +* A SpecMap is a specialised form of Mapping which can be used to +* represent a sequence of conversions between standard spectral +* coordinate systems. +* +* When an SpecMap is first created, it simply performs a unit +c (null) Mapping. Using the astSpecAdd +f (null) Mapping. Using the AST_SPECADD +c function, a series of coordinate conversion steps may then be +f routine, a series of coordinate conversion steps may then be +* added. This allows multi-step conversions between a variety of +* spectral coordinate systems to be assembled out of a set of building +* blocks. +* +* Conversions are available to transform between standards of rest. +* Such conversions need to know the source position as an RA and DEC. +* This information can be supplied in the form of parameters for +* the relevant conversions, in which case the SpecMap is 1-dimensional, +* simply transforming the spectral axis values. This means that the +* same source position will always be used by the SpecMap. However, this +* may not be appropriate for an accurate description of a 3-D spectral +* cube, where changes of spatial position can produce significant +* changes in the Doppler shift introduced when transforming between +* standards of rest. For this situation, a 3-dimensional SpecMap can +* be created in which axes 2 and 3 correspond to the source RA and DEC +* The SpecMap simply copies values for axes 2 and 3 from input to +* output), but modifies axis 1 values (the spectral axis) appropriately. +* +* For details of the individual coordinate conversions available, +c see the description of the astSpecAdd function. +f see the description of the AST_SPECADD routine. + +* Inheritance: +* The SpecMap class inherits from the Mapping class. + +* Attributes: +* The SpecMap class does not define any new attributes beyond those +* which are applicable to all Mappings. + +* Functions: +c In addition to those functions applicable to all Mappings, the +c following function may also be applied to all SpecMaps: +f In addition to those routines applicable to all Mappings, the +f following routine may also be applied to all SpecMaps: +* +c - astSpecAdd: Add a spectral coordinate conversion to an SpecMap +f - AST_SPECADD: Add a spectral coordinate conversion to an SpecMap + +* Copyright: +* Copyright (C) 1997-2006 Council for the Central Laboratory of the +* Research Councils +* Copyright (C) 2009 Science & Technology Facilities Council. +* All Rights Reserved. + +* Licence: +* This program is free software: you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation, either +* version 3 of the License, or (at your option) any later +* version. +* +* This program 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 Lesser General Public License for more details. +* +* You should have received a copy of the GNU Lesser General +* License along with this program. If not, see +* <http://www.gnu.org/licenses/>. + +* Authors: +* DSB: David S. Berry (Starlink) + +* History: +* 6-NOV-2002 (DSB): +* Original version. +* 14-JUL-2003 (DSB): +* Added checks for NAN values produced by transformation functions. +* 17-SEP-2003 (DSB): +* - Improve FRTOAW accuracy by iterating. +* - Changed Refrac to use algorithm given in FITS-WCS paper 3 +* version dated 21/9/03. +* 14-FEB-2006 (DSB): +* Override astGetObjSize. +* 10-MAY-2006 (DSB): +* Override astEqual. +* 15-NOV-2006 (DSB): +* Guard against division by zero when converting freq to wave in +* SystemChange. +* 18-JUN-2009 (DSB): +* Add OBSALT argument to TPF2HL and HLF2TP conversions. +* Change GEOLON/LAT to OBSLON/LAT for consistency with other +* classes. +* 2-OCT-2012 (DSB): +* Check for Infs as well as NaNs. +*class-- +*/ + +/* Module Macros. */ +/* ============== */ +/* Set the name of the class we are implementing. This indicates to + the header files that define class interfaces that they should make + "protected" symbols available. */ +#define astCLASS SpecMap + +/* Codes to identify spectral coordinate conversions. */ +#define AST__SPEC_NULL 0 /* Null value */ +#define AST__FRTOVL 1 /* Frequency to relativistic velocity */ +#define AST__VLTOFR 2 /* Relativistic velocity to Frequency */ +#define AST__ENTOFR 3 /* Energy to frequency */ +#define AST__FRTOEN 4 /* Frequency to energy */ +#define AST__WNTOFR 5 /* Wave number to frequency */ +#define AST__FRTOWN 6 /* Frequency to wave number */ +#define AST__WVTOFR 7 /* Wavelength (vacuum) to frequency */ +#define AST__FRTOWV 8 /* Frequency to wavelength (vacuum) */ +#define AST__AWTOFR 9 /* Wavelength (air) to frequency */ +#define AST__FRTOAW 10 /* Frequency to wavelength (air) */ +#define AST__VRTOVL 11 /* Radio to relativistic velocity */ +#define AST__VLTOVR 12 /* Relativistic to radio velocity */ +#define AST__VOTOVL 13 /* Optical to relativistic velocity */ +#define AST__VLTOVO 14 /* Relativistic to optical velocity */ +#define AST__ZOTOVL 15 /* Redshift to relativistic velocity */ +#define AST__VLTOZO 16 /* Relativistic velocity to redshift */ +#define AST__BTTOVL 17 /* Beta factor to relativistic velocity */ +#define AST__VLTOBT 18 /* Relativistic velocity to beta factor */ +#define AST__USF2HL 19 /* User-defined to heliocentric frequency */ +#define AST__HLF2US 20 /* Heliocentric to user-defined frequency */ +#define AST__TPF2HL 21 /* Topocentric to heliocentric frequency */ +#define AST__HLF2TP 22 /* Heliocentric to topocentric frequency */ +#define AST__GEF2HL 23 /* Geocentric to heliocentric frequency */ +#define AST__HLF2GE 24 /* Heliocentric to geocentric frequency */ +#define AST__BYF2HL 25 /* Barycentric to heliocentric frequency */ +#define AST__HLF2BY 26 /* Heliocentric to barycentric frequency */ +#define AST__LKF2HL 27 /* LSRK to heliocentric frequency */ +#define AST__HLF2LK 28 /* Heliocentric to LSRK frequency */ +#define AST__LDF2HL 29 /* LSRD to heliocentric frequency */ +#define AST__HLF2LD 30 /* Heliocentric to LSRD frequency */ +#define AST__LGF2HL 31 /* Local group to heliocentric frequency */ +#define AST__HLF2LG 32 /* Heliocentric to local group frequency */ +#define AST__GLF2HL 33 /* Galactic to heliocentric frequency */ +#define AST__HLF2GL 34 /* Heliocentric to galactic frequency */ + +/* Maximum number of arguments required by a conversion. */ +#define MAX_ARGS 7 + +/* The alphabet (used for generating keywords for arguments). */ +#define ALPHABET "abcdefghijklmnopqrstuvwxyz" + +/* Angle conversion */ +#define PI 3.141592653589793238462643 +#define PIBY2 (PI/2.0) +#define D2R (PI/180.0) +#define R2D (180.0/PI) + +/* Include files. */ +/* ============== */ +/* Interface definitions. */ +/* ---------------------- */ +#include "pal.h" /* SLALIB interface */ + +#include "globals.h" /* Thread-safe global data access */ +#include "error.h" /* Error reporting facilities */ +#include "memory.h" /* Memory allocation facilities */ +#include "object.h" /* Base Object class */ +#include "pointset.h" /* Sets of points/coordinates */ +#include "mapping.h" /* Coordinate Mappings (parent class) */ +#include "unitmap.h" /* Unit (null) Mappings */ +#include "specmap.h" /* Interface definition for this class */ + +/* Error code definitions. */ +/* ----------------------- */ +#include "ast_err.h" /* AST error codes */ + +/* C header files. */ +/* --------------- */ +#include <ctype.h> +#include <stddef.h> +#include <stdio.h> +#include <string.h> +#include <math.h> + +/* Module Variables. */ +/* ================= */ + +/* Address of this static variable is used as a unique identifier for + member of this class. */ +static int class_check; + +/* Pointers to parent class methods which are extended by this class. */ +static int (* parent_getobjsize)( AstObject *, int * ); +static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * ); +static double (* parent_rate)( AstMapping *, double *, int, int, int * ); + + +#ifdef THREAD_SAFE +/* Define how to initialise thread-specific globals. */ +#define GLOBAL_inits \ + globals->Class_Init = 0; + +/* Create the function that initialises global data for this module. */ +astMAKE_INITGLOBALS(SpecMap) + +/* Define macros for accessing each item of thread specific global data. */ +#define class_init astGLOBAL(SpecMap,Class_Init) +#define class_vtab astGLOBAL(SpecMap,Class_Vtab) + + +#include <pthread.h> + + +#else + + +/* Define the class virtual function table and its initialisation flag + as static variables. */ +static AstSpecMapVtab class_vtab; /* Virtual function table */ +static int class_init = 0; /* Virtual function table initialised? */ + +#endif + +/* Structure to hold parameters and intermediate values describing a + reference frame */ +typedef struct FrameDef { + double obsalt; /* Observers geodetic altitude (m) */ + double obslat; /* Observers geodetic latitude (rads) */ + double obslon; /* Observers geodetic longitude (rads, +ve east) */ + double epoch; /* Julian epoch of observation */ + double refdec; /* RA of reference point (FK5 J2000) */ + double refra; /* DEC of reference point (FK5 J2000) */ + double veluser; /* Heliocentric velocity of user-defined system (m/s) */ + double last; /* Local apparent sideral time */ + double amprms[21]; /* Mean to apparent parameters */ + double vuser[3]; /* Used-defined velocity as a FK5 J2000 vector */ + double dvh[3]; /* Earth-sun velocity */ + double dvb[3]; /* Barycentre-sun velocity */ +} FrameDef; + +/* External Interface Function Prototypes. */ +/* ======================================= */ +/* The following functions have public prototypes only (i.e. no + protected prototypes), so we must provide local prototypes for use + within this module. */ +AstSpecMap *astSpecMapId_( int, int, const char *, ... ); + +/* Prototypes for Private Member Functions. */ +/* ======================================== */ +static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * ); +static const char *CvtString( int, const char **, int *, int *, int *, int *, const char *[ MAX_ARGS ], int * ); +static double BaryVel( double, double, FrameDef *, int * ); +static double GalVel( double, double, FrameDef *, int * ); +static double GeoVel( double, double, FrameDef *, int * ); +static double LgVel( double, double, FrameDef *, int * ); +static double LsrdVel( double, double, FrameDef *, int * ); +static double LsrkVel( double, double, FrameDef *, int * ); +static double Rate( AstMapping *, double *, int, int, int * ); +static double Refrac( double, int * ); +static double Rverot( double, double, double, double, double, int * ); +static double TopoVel( double, double, FrameDef *, int * ); +static double UserVel( double, double, FrameDef *, int * ); +static int CvtCode( const char *, int * ); +static int Equal( AstObject *, AstObject *, int * ); +static int FrameChange( int, int, double *, double *, double *, double *, int, int * ); +static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * ); +static int SystemChange( int, int, double *, double *, int, int * ); +static void AddSpecCvt( AstSpecMap *, int, const double *, int * ); +static void Copy( const AstObject *, AstObject *, int * ); +static void Delete( AstObject *, int * ); +static void Dump( AstObject *, AstChannel *, int * ); +static void SpecAdd( AstSpecMap *, const char *, const double[], int * ); + +static int GetObjSize( AstObject *, int * ); +/* Member functions. */ +/* ================= */ +static int Equal( AstObject *this_object, AstObject *that_object, int *status ) { +/* +* Name: +* Equal + +* Purpose: +* Test if two SpecMaps are equivalent. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* int Equal( AstObject *this, AstObject *that, int *status ) + +* Class Membership: +* SpecMap member function (over-rides the astEqual protected +* method inherited from the astMapping class). + +* Description: +* This function returns a boolean result (0 or 1) to indicate whether +* two SpecMaps are equivalent. + +* Parameters: +* this +* Pointer to the first Object (a SpecMap). +* that +* Pointer to the second Object. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* One if the SpecMaps are equivalent, zero otherwise. + +* Notes: +* - A value of zero will be returned if this function is invoked +* with the global status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + AstSpecMap *that; + AstSpecMap *this; + const char *argdesc[ MAX_ARGS ]; + const char *comment; + int argdec; + int argra; + int i, j; + int nargs; + int nin; + int nout; + int result; + int szargs; + +/* Initialise. */ + result = 0; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Obtain pointers to the two SpecMap structures. */ + this = (AstSpecMap *) this_object; + that = (AstSpecMap *) that_object; + +/* Check the second object is a SpecMap. We know the first is a + SpecMap since we have arrived at this implementation of the virtual + function. */ + if( astIsASpecMap( that ) ) { + +/* Get the number of inputs and outputs and check they are the same for both. */ + nin = astGetNin( this ); + nout = astGetNout( this ); + if( astGetNin( that ) == nin && astGetNout( that ) == nout ) { + +/* If the Invert flags for the two SpecMaps differ, it may still be possible + for them to be equivalent. First compare the SpecMaps if their Invert + flags are the same. In this case all the attributes of the two SpecMaps + must be identical. */ + if( astGetInvert( this ) == astGetInvert( that ) ) { + if( this->ncvt == that->ncvt ) { + result = 1; + for( i = 0; i < this->ncvt && result; i++ ) { + if( this->cvttype[ i ] != that->cvttype[ i ] ) { + result = 0; + } else { + CvtString( this->cvttype[ i ], &comment, &argra, + &argdec, &nargs, &szargs, argdesc, status ); + for( j = 0; j < nargs; j++ ) { + if( !astEQUAL( this->cvtargs[ i ][ j ], + that->cvtargs[ i ][ j ] ) ){ + result = 0; + break; + } + } + } + } + } + +/* If the Invert flags for the two SpecMaps differ, the attributes of the two + SpecMaps must be inversely related to each other. */ + } else { + +/* In the specific case of a SpecMap, Invert flags must be equal. */ + result = 0; + + } + } + } + +/* If an error occurred, clear the result value. */ + if ( !astOK ) result = 0; + +/* Return the result, */ + return result; +} + +static int GetObjSize( AstObject *this_object, int *status ) { +/* +* Name: +* GetObjSize + +* Purpose: +* Return the in-memory size of an Object. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* int GetObjSize( AstObject *this, int *status ) + +* Class Membership: +* SpecMap member function (over-rides the astGetObjSize protected +* method inherited from the parent class). + +* Description: +* This function returns the in-memory size of the supplied SpecMap, +* in bytes. + +* Parameters: +* this +* Pointer to the SpecMap. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* The Object size, in bytes. + +* Notes: +* - A value of zero will be returned if this function is invoked +* with the global status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + AstSpecMap *this; /* Pointer to SpecMap structure */ + int result; /* Result value to return */ + int cvt; /* Loop counter for coordinate conversions */ + +/* Initialise. */ + result = 0; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Obtain a pointers to the SpecMap structure. */ + this = (AstSpecMap *) this_object; + +/* Invoke the GetObjSize method inherited from the parent class, and then + add on any components of the class structure defined by thsi class + which are stored in dynamically allocated memory. */ + result = (*parent_getobjsize)( this_object, status ); + + for ( cvt = 0; cvt < this->ncvt; cvt++ ) { + result += astTSizeOf( this->cvtargs[ cvt ] ); + } + + result += astTSizeOf( this->cvtargs ); + result += astTSizeOf( this->cvttype ); + +/* If an error occurred, clear the result value. */ + if ( !astOK ) result = 0; + +/* Return the result, */ + return result; +} + +static void AddSpecCvt( AstSpecMap *this, int cvttype, const double *args, int *status ) { +/* +* Name: +* AddSpecCvt + +* Purpose: +* Add a coordinate conversion step to an SpecMap. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* void AddSpecCvt( AstSpecMap *this, int cvttype, const double *args ) + +* Class Membership: +* SpecMap member function. + +* Description: +* This function allows one of the supported spectral coordinate +* conversions to be appended to a SpecMap. When a SpecMap is first +* created (using astSpecMap), it simply performs a unit mapping. By +* using AddSpecCvt repeatedly, a series of coordinate conversions may +* then be specified which the SpecMap will subsequently perform in +* sequence. This allows a complex coordinate conversion to be +* assembled out of the basic building blocks. The SpecMap will also +* perform the inverse coordinate conversion (applying the individual +* conversion steps in reverse) if required. + +* Parameters: +* this +* Pointer to the SpecMap. +* cvttype +* A code to identify which spectral coordinate conversion is to be +* appended. See the "Coordinate Conversions" section for details +* of those available. +* args +* Pointer to an array of double containing the argument values +* required to fully specify the required coordinate +* conversion. The number of arguments depends on the conversion +* (see the "Coordinate Conversions" section for details). This +* value is ignored and may be NULL if no arguments are required. + +* Returned Value: +* void. + +* Coordinate Conversions: +* The following values may be supplied for the "cvttype" parameter +* in order to specify the coordinate conversion to be performed. +* The argument(s) required to fully specify each conversion are +* indicated in parentheses after each value, and described at the end +* of the list. Values for these should be given in the array pointed +* at by "args". +* +* AST__FRTOVL( RF ) +* Convert frequency to relativistic velocity. +* AST__VLTOFR( RF ) +* Convert relativistic velocity to Frequency. +* AST__ENTOFR +* Convert energy to frequency. +* AST__FRTOEN +* Convert frequency to energy. +* AST__WNTOFR +* Convert wave number to frequency. +* AST__FRTOWN +* Convert frequency to wave number. +* AST__WVTOFR +* Convert wavelength (vacuum) to frequency. +* AST__FRTOWV +* Convert frequency to wavelength (vacuum). +* AST__AWTOFR +* Convert wavelength (air) to frequency. +* AST__FRTOAW +* Convert frequency to wavelength (air). +* AST__VRTOVL +* Convert radio to relativistic velocity. +* AST__VLTOVR +* Convert relativistic to radio velocity. +* AST__VOTOVL +* Convert optical to relativistic velocity. +* AST__VLTOVO +* Convert relativistic to optical velocity. +* AST__ZOTOVL +* Convert redshift to relativistic velocity. +* AST__VLTOZO +* Convert relativistic velocity to redshift. +* AST__BTTOVL +* Convert beta factor to relativistic velocity. +* AST__VLTOBT +* Convert relativistic velocity to beta factor. +* AST_USF2HL( VOFF, RA, DEC ) +* Convert frequency from a user-defined reference frame to +* heliocentric. +* AST__HLF2US( VOFF, RA, DEC ) +* Convert frequency from heliocentric reference frame to +* user-defined. +* AST__TPF2HL( OBSLON, OBSLAT, OBSALT, EPOCH, RA, DEC ) +* Convert from Topocentric to heliocentric frequency +* AST__HLF2TP( OBSLON, OBSLAT, OBSALT, EPOCH, RA, DEC ) +* Convert from Heliocentric to topocentric frequency. +* AST__GEF2HL( EPOCH, RA, DEC ) +* Convert from Geocentric to heliocentric frequency. +* AST__HLF2GE( EPOCH, RA, DEC ) +* Convert from Heliocentric to geocentric frequency. +* AST__BYF2HL( EPOCH, RA, DEC ) +* Convert from Barycentric to heliocentric frequency. +* AST__HLF2BY( EPOCH, RA, DEC ) +* Convert from Heliocentric to barycentric frequency. +* AST__LKF2HL( RA, DEC ) +* Convert from LSRK to heliocentric frequency. +* AST__HLF2LK( RA, DEC ) +* Convert from Heliocentric to LSRK frequency. +* AST__LDF2HL( RA, DEC ) +* Convert from LSRD to heliocentric frequency. +* AST__HLF2LD( RA, DEC ) +* Convert from Heliocentric to LSRD frequency. +* AST__LGF2HL( RA, DEC ) +* Convert from Local group to heliocentric frequency. +* AST__HLF2LG( RA, DEC ) +* Convert from Heliocentric to local group frequency. +* AST__GLF2HL( RA, DEC ) +* Convert from Galactic to heliocentric frequency. +* AST__HLF2GL( RA, DEC ) +* Convert from Heliocentric to galactic frequency. +* +* The units for the values processed by the above conversions are as +* follows: +* +* - all velocities: metres per second. +* - frequency: Hertz. +* - all wavelengths: metres. +* - energy: Joules. +* - wave number: cycles per metre. +* +* The arguments used in the above conversions are as follows: +* +* - RF: Rest frequency (Hz). +* - OBSALT: Geodetic altitude of observer (IAU 1975, metres). +* - OBSLAT: Geodetic latitude of observer (IAU 1975, radians). +* - OBSLON: Longitude of observer (radians, positive eastwards). +* - EPOCH: Epoch of observation (UT1 expressed as a Modified Julian Date). +* - RA: Right Ascension of source (radians, FK5 J2000). +* - DEC: Declination of source (radians, FK5 J2000). +* - VOFF: Velocity of the user-defined reference frame, towards the +* position given by RA and DEC, measured in the heliocentric +* reference frame. +* +* If the SpecMap is 3-dimensional, source positions are provided by the +* values supplied to inputs 2 and 3 of the SpecMap (which are simply +* copied to outputs 2 and 3). Note, usable values are still required +* for the RA and DEC arguments in order to define the "user-defined" +* reference frame used by USF2HL and HLF2US. However, AST__BAD can be +* supplied for RA and DEC if the user-defined reference frame is not +* required. + +* Notes: +* - The specified conversion is appended only if the SpecMap's +* Invert attribute is zero. If it is non-zero, this function +* effectively prefixes the inverse of the conversion specified +* instead. +*/ + +/* Local Variables: */ + const char *argdesc[ MAX_ARGS ]; /* Pointers to argument descriptions */ + const char *comment; /* Pointer to comment string */ + const char *cvt_string; /* Pointer to conversion type string */ + int argdec; /* Index of DEC argument */ + int argra; /* Index of RA argument */ + int i; /* Argument index */ + int nargs; /* Number of user-supplied arguments */ + int ncvt; /* Number of coordinate conversions */ + int szargs; /* Size of arguments array */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Validate the coordinate conversion type and obtain the number of + required user-supplied arguments, and the size of the array in which + to put the user-supplied arguments (the array meay leave room after + the user-supplied arguments for various useful pre-calculated values). */ + cvt_string = CvtString( cvttype, &comment, &argra, &argdec, &nargs, + &szargs, argdesc, status ); + +/* If the coordinate conversion type was not valid, then report an + error. */ + if ( astOK && !cvt_string ) { + astError( AST__SPCIN, "AddSpecCvt(%s): Invalid spectral coordinate " + "conversion type (%d).", status, astGetClass( this ), + (int) cvttype ); + } + +/* Note the number of coordinate conversions already stored in the SpecMap. */ + if ( astOK ) { + ncvt = this->ncvt; + +/* Extend the array of conversion types and the array of pointers to + their argument lists to accommodate the new one. */ + this->cvttype = (int *) astGrow( this->cvttype, ncvt + 1, + sizeof( int ) ); + this->cvtargs = (double **) astGrow( this->cvtargs, ncvt + 1, + sizeof( double * ) ); + +/* If OK, allocate memory and store a copy of the argument list, + putting a pointer to the copy into the SpecMap. */ + if ( astOK ) { + this->cvtargs[ ncvt ] = astStore( NULL, args, + sizeof( double ) * (size_t) szargs ); + } + +/* Store the conversion type and increment the conversion count. Also put + AST__BAD in any elements of the argument array which are beyond the + end of the user-supplied arguments. These will be used to hold + intermediate values calculated on the basis of the user-supplied + arguments. */ + if ( astOK ) { + this->cvttype[ ncvt ] = cvttype; + this->ncvt++; + for( i = nargs; i < szargs; i++ ) this->cvtargs[ ncvt ][ i ] = AST__BAD; + } + } +} + +static double BaryVel( double ra, double dec, FrameDef *def, int *status ) { +/* +* Name: +* BaryVel + +* Purpose: +* Find the velocity of the earth-sun barycentre away from the source. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* double BaryVel( double ra, double dec, FrameDef *def, int *status ) + +* Class Membership: +* SpecMap method. + +* Description: +* This function finds the component of the velocity of the earth-sun +* barycentre away from a specified source position, at a given epoch, in +* the frame of rest of the centre of the Sun. + +* Parameters: +* ra +* The RA (rads, FK5 J2000) of the source. +* dec +* The Dec (rads, FK5 J2000) of the source. +* def +* Pointer to a FrameDef structure which holds the parameters which +* define the frame, together with cached intermediate results. +* status +* Pointer to the inherited status variable. + +* Returns: +* The component of the frame's velocity away from the position given by +* "ra" and "dec", in m/s, measured within the Heliographic frame of +* rest. Zero is returned if an error has already occurred. + +*/ + +/* Local Variables: */ + double dpb[ 3 ]; /* Barycentric earth position vector */ + double dph[ 3 ]; /* Heliocentric earth position vector */ + double dvh[ 3 ]; /* Heliocentric earth velocity vector */ + double v[ 3 ]; /* Source direction vector */ + +/* Check the global error status. */ + if ( !astOK ) return 0.0; + +/* Get the Cartesian vector towards the source, in the Cartesian FK5 + J2000 system. */ + palDcs2c( ra, dec, v ); + +/* If not already done so, get the Earth/Sun velocity and position vectors in + the same system. Speed is returned in units of AU/s. Store in the supplied + frame definition structure. */ + if( def->dvb[ 0 ] == AST__BAD ) { + palEvp( def->epoch, 2000.0, def->dvb, dpb, dvh, dph ); + +/* Change the barycentric velocity of the earth into the heliocentric + velocity of the barycentre. */ + def->dvb[ 0 ] = dvh[ 0 ] - def->dvb[ 0 ]; + def->dvb[ 1 ] = dvh[ 1 ] - def->dvb[ 1 ]; + def->dvb[ 2 ] = dvh[ 2 ] - def->dvb[ 2 ]; + } + +/* Return the component away from the source, of the velocity of the + barycentre relative to the sun (in m/s). */ + return -palDvdv( v, def->dvb )*149.597870E9; + +} + +static int CvtCode( const char *cvt_string, int *status ) { +/* +* Name: +* CvtCode + +* Purpose: +* Convert a conversion type from a string representation to a code value. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* int CvtCode( const char *cvt_string, int *status ) + +* Class Membership: +* SpecMap member function. + +* Description: +* This function accepts a string used to repersent one of the +* SpecMap coordinate conversions and converts it into a code +* value for internal use. + +* Parameters: +* cvt_string +* Pointer to a constant null-terminated string representing a +* spectral coordinate conversion. This is case sensitive and should +* contain no unnecessary white space. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* The equivalent conversion code. If the string was not +* recognised, the code AST__SPEC_NULL is returned, without error. + +* Notes: +* - A value of AST__SPEC_NULL will be returned if this function is +* invoked with the global error status set, or if it should fail +* for any reason. +*/ + +/* Local Variables: */ + int result; /* Result value to return */ + +/* Initialise. */ + result = AST__SPEC_NULL; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Test the string against each recognised value in turn and assign + the result. */ + if ( astChrMatch( cvt_string, "FRTOVL" ) ) { + result = AST__FRTOVL; + + } else if ( astChrMatch( cvt_string, "VLTOFR" ) ) { + result = AST__VLTOFR; + + } else if ( astChrMatch( cvt_string, "VLTOFR" ) ) { + result = AST__VLTOFR; + + } else if ( astChrMatch( cvt_string, "ENTOFR" ) ) { + result = AST__ENTOFR; + + } else if ( astChrMatch( cvt_string, "FRTOEN" ) ) { + result = AST__FRTOEN; + + } else if ( astChrMatch( cvt_string, "WNTOFR" ) ) { + result = AST__WNTOFR; + + } else if ( astChrMatch( cvt_string, "FRTOWN" ) ) { + result = AST__FRTOWN; + + } else if ( astChrMatch( cvt_string, "WVTOFR" ) ) { + result = AST__WVTOFR; + + } else if ( astChrMatch( cvt_string, "FRTOWV" ) ) { + result = AST__FRTOWV; + + } else if ( astChrMatch( cvt_string, "AWTOFR" ) ) { + result = AST__AWTOFR; + + } else if ( astChrMatch( cvt_string, "FRTOAW" ) ) { + result = AST__FRTOAW; + + } else if ( astChrMatch( cvt_string, "VRTOVL" ) ) { + result = AST__VRTOVL; + + } else if ( astChrMatch( cvt_string, "VLTOVR" ) ) { + result = AST__VLTOVR; + + } else if ( astChrMatch( cvt_string, "VOTOVL" ) ) { + result = AST__VOTOVL; + + } else if ( astChrMatch( cvt_string, "VLTOVO" ) ) { + result = AST__VLTOVO; + + } else if ( astChrMatch( cvt_string, "ZOTOVL" ) ) { + result = AST__ZOTOVL; + + } else if ( astChrMatch( cvt_string, "VLTOZO" ) ) { + result = AST__VLTOZO; + + } else if ( astChrMatch( cvt_string, "BTTOVL" ) ) { + result = AST__BTTOVL; + + } else if ( astChrMatch( cvt_string, "VLTOBT" ) ) { + result = AST__VLTOBT; + + } else if ( astChrMatch( cvt_string, "USF2HL" ) ) { + result = AST__USF2HL; + + } else if ( astChrMatch( cvt_string, "HLF2US" ) ) { + result = AST__HLF2US; + + } else if ( astChrMatch( cvt_string, "TPF2HL" ) ) { + result = AST__TPF2HL; + + } else if ( astChrMatch( cvt_string, "HLF2TP" ) ) { + result = AST__HLF2TP; + + } else if ( astChrMatch( cvt_string, "GEF2HL" ) ) { + result = AST__GEF2HL; + + } else if ( astChrMatch( cvt_string, "HLF2GE" ) ) { + result = AST__HLF2GE; + + } else if ( astChrMatch( cvt_string, "BYF2HL" ) ) { + result = AST__BYF2HL; + + } else if ( astChrMatch( cvt_string, "HLF2BY" ) ) { + result = AST__HLF2BY; + + } else if ( astChrMatch( cvt_string, "LKF2HL" ) ) { + result = AST__LKF2HL; + + } else if ( astChrMatch( cvt_string, "HLF2LK" ) ) { + result = AST__HLF2LK; + + } else if ( astChrMatch( cvt_string, "LDF2HL" ) ) { + result = AST__LDF2HL; + + } else if ( astChrMatch( cvt_string, "HLF2LD" ) ) { + result = AST__HLF2LD; + + } else if ( astChrMatch( cvt_string, "LGF2HL" ) ) { + result = AST__LGF2HL; + + } else if ( astChrMatch( cvt_string, "HLF2LG" ) ) { + result = AST__HLF2LG; + + } else if ( astChrMatch( cvt_string, "GLF2HL" ) ) { + result = AST__GLF2HL; + + } else if ( astChrMatch( cvt_string, "HLF2GL" ) ) { + result = AST__HLF2GL; + + } + +/* Return the result. */ + return result; +} + +static const char *CvtString( int cvt_code, const char **comment, + int *argra, int *argdec, int *nargs, int *szargs, + const char *arg[ MAX_ARGS ], int *status ) { +/* +* Name: +* CvtString + +* Purpose: +* Convert a conversion type from a code value to a string representation. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* const char *CvtString( int cvt_code, const char **comment, +* int *argra, int *argdec, int *nargs, +* int *szargs, const char *arg[ MAX_ARGS ], int *status ) + +* Class Membership: +* SpecMap member function. + +* Description: +* This function accepts a code value used to represent one of the +* SpecMap coordinate conversions and converts it into an +* equivalent string representation. It also returns a descriptive +* comment and information about the arguments required in order to +* perform the conversion. + +* Parameters: +* cvt_code +* The conversion code. +* comment +* Address of a location to return a pointer to a constant +* null-terminated string containing a description of the +* conversion. +* argra +* Address of an int in which to return the index of the argument +* corresponding to the source RA. Returned equal to -1 if the +* conversion does not have a source RA argument. +* argdec +* Address of an int in which to return the index of the argument +* corresponding to the source DEC. Returned equal to -1 if the +* conversion does not have a source DEC argument. +* nargs +* Address of an int in which to return the number of arguments +* required from the user in order to perform the conversion (may +* be zero). +* szargs +* Address of an int in which to return the number of arguments +* associated with the conversion. This may be bigger than "nargs" +* if the conversion can pre-calculate useful values on the basis +* of the user-supplied values. Such precalculated values are +* stored after the last user-supplied argument. +* arg +* An array in which to return a pointer to a constant +* null-terminated string for each argument (above) containing a +* description of what each argument represents. This includes both +* user-supplied arguments and pre-calculated values. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Pointer to a constant null-terminated string representation of +* the conversion code value supplied. If the code supplied is not +* valid, a NULL pointer will be returned, without error. + +* Notes: +* - A NULL pointer value will be returned if this function is +* invoked with the global error status set, or if it should fail +* for any reason. +*/ + +/* Local Variables: */ + const char *result; /* Result pointer to return */ + +/* Initialise the returned values. */ + *comment = NULL; + *nargs = 0; + *argra = -1; + *argdec = -1; + result = NULL; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Test for each valid code value in turn and assign the appropriate + return values. */ + switch ( cvt_code ) { + + case AST__FRTOVL: + *comment = "Convert frequency to rel. velocity"; + result = "FRTOVL"; + *nargs = 1; + *szargs = 1; + arg[ 0 ] = "Rest frequency (Hz)"; + break; + + case AST__VLTOFR: + *comment = "Convert rel. velocity to frequency"; + result = "VLTOFR"; + *nargs = 1; + *szargs = 1; + arg[ 0 ] = "Rest frequency (Hz)"; + break; + + case AST__ENTOFR: + *comment = "Convert energy to frequency"; + result = "ENTOFR"; + *nargs = 0; + *szargs = 0; + break; + + case AST__FRTOEN: + *comment = "Convert frequency to energy"; + result = "FRTOEN"; + *nargs = 0; + *szargs = 0; + break; + + case AST__WNTOFR: + *comment = "Convert wave number to frequency"; + result = "WNTOFR"; + *nargs = 0; + *szargs = 0; + break; + + case AST__FRTOWN: + *comment = "Convert frequency to wave number"; + result = "FRTOWN"; + *nargs = 0; + *szargs = 0; + break; + + case AST__WVTOFR: + *comment = "Convert wavelength (vacuum) to frequency"; + result = "WVTOFR"; + *nargs = 0; + *szargs = 0; + break; + + case AST__FRTOWV: + *comment = "Convert frequency to wavelength (vacuum)"; + result = "FRTOWV"; + *nargs = 0; + *szargs = 0; + break; + + case AST__AWTOFR: + *comment = "Convert wavelength (air) to frequency"; + result = "AWTOFR"; + *nargs = 0; + *szargs = 0; + break; + + case AST__FRTOAW: + *comment = "Convert frequency to wavelength (air)"; + result = "FRTOAW"; + *nargs = 0; + *szargs = 0; + break; + + case AST__VRTOVL: + *comment = "Convert radio to rel. velocity"; + result = "VRTOVL"; + *nargs = 0; + *szargs = 0; + break; + + case AST__VLTOVR: + *comment = "Convert relativistic to radio velocity"; + result = "VLTOVR"; + *nargs = 0; + *szargs = 0; + break; + + case AST__VOTOVL: + *comment = "Convert optical to rel. velocity"; + result = "VOTOVL"; + *nargs = 0; + *szargs = 0; + break; + + case AST__VLTOVO: + *comment = "Convert relativistic to optical velocity"; + result = "VLTOVO"; + *nargs = 0; + *szargs = 0; + break; + + case AST__ZOTOVL: + *comment = "Convert redshift to rel. velocity"; + result = "ZOTOVL"; + *nargs = 0; + *szargs = 0; + break; + + case AST__VLTOZO: + *comment = "Convert rel. velocity to redshift"; + result = "VLTOZO"; + *nargs = 0; + *szargs = 0; + break; + + case AST__BTTOVL: + *comment = "Convert beta factor to rel. velocity"; + result = "BTTOVL"; + *nargs = 0; + *szargs = 0; + break; + + case AST__VLTOBT: + *comment = "Convert rel. velocity to beta factor"; + result = "VLTOBT"; + *nargs = 0; + *szargs = 0; + break; + + case AST__USF2HL: + *comment = "Convert from user-defined to heliocentric frequency"; + result = "USF2HL"; + *argra = 1; + *argdec = 2; + *nargs = 3; + *szargs = 4; + arg[ 0 ] = "Velocity offset (m/s)"; + arg[ 1 ] = "RA of source (FK5 J2000, radians)"; + arg[ 2 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 3 ] = "Frequency correction factor"; + break; + + case AST__HLF2US: + *comment = "Convert from heliocentric to user-defined frequency"; + result = "HLF2US"; + *argra = 1; + *argdec = 2; + *nargs = 3; + *szargs = 4; + arg[ 0 ] = "Velocity offset (m/s)"; + arg[ 1 ] = "RA of source (FK5 J2000, radians)"; + arg[ 2 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 3 ] = "Frequency correction factor"; + break; + + case AST__TPF2HL: + *comment = "Convert from Topocentric to heliocentric frequency"; + result = "TPF2HL"; + *argra = 4; + *argdec = 5; + *nargs = 6; + *szargs = 7; + arg[ 0 ] = "Longitude (positive eastwards, radians)"; + arg[ 1 ] = "Latitude (geodetic, radians)"; + arg[ 2 ] = "Altitude (geodetic, metres)"; + arg[ 3 ] = "UT1 epoch of observaton (Modified Julian Date)"; + arg[ 4 ] = "RA of source (FK5 J2000, radians)"; + arg[ 5 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 6 ] = "Frequency correction factor"; + break; + + case AST__HLF2TP: + *comment = "Convert from Heliocentric to topocentric frequency"; + result = "HLF2TP"; + *argra = 4; + *argdec = 5; + *nargs = 6; + *szargs = 7; + arg[ 0 ] = "Longitude (positive eastwards, radians)"; + arg[ 1 ] = "Latitude (geodetic, radians)"; + arg[ 2 ] = "Altitude (geodetic, metres)"; + arg[ 3 ] = "UT1 epoch of observaton (Modified Julian Date)"; + arg[ 4 ] = "RA of source (FK5 J2000, radians)"; + arg[ 5 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 6 ] = "Frequency correction factor"; + break; + + case AST__GEF2HL: + *comment = "Convert from Geocentric to heliocentric frequency"; + result = "GEF2HL"; + *argra = 1; + *argdec = 2; + *nargs = 3; + *szargs = 4; + arg[ 0 ] = "UT1 epoch of observaton (Modified Julian Date)"; + arg[ 1 ] = "RA of source (FK5 J2000, radians)"; + arg[ 2 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 3 ] = "Frequency correction factor"; + break; + + case AST__HLF2GE: + *comment = "Convert from Heliocentric to geocentric frequency"; + result = "HLF2GE"; + *argra = 1; + *argdec = 2; + *nargs = 3; + *szargs = 4; + arg[ 0 ] = "UT1 epoch of observaton (Modified Julian Date)"; + arg[ 1 ] = "RA of source (FK5 J2000, radians)"; + arg[ 2 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 3 ] = "Frequency correction factor"; + break; + + case AST__BYF2HL: + *comment = "Convert from Barycentric to heliocentric frequency"; + result = "BYF2HL"; + *argra = 1; + *argdec = 2; + *nargs = 3; + *szargs = 4; + arg[ 0 ] = "UT1 epoch of observaton (Modified Julian Date)"; + arg[ 1 ] = "RA of source (FK5 J2000, radians)"; + arg[ 2 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 3 ] = "Frequency correction factor"; + break; + + case AST__HLF2BY: + *comment = "Convert from Heliocentric to barycentric frequency"; + result = "HLF2BY"; + *argra = 1; + *argdec = 2; + *nargs = 3; + *szargs = 4; + arg[ 0 ] = "UT1 epoch of observaton (Modified Julian Date)"; + arg[ 1 ] = "RA of source (FK5 J2000, radians)"; + arg[ 2 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 3 ] = "Frequency correction factor"; + break; + + case AST__LKF2HL: + *comment = "Convert from LSRK to heliocentric frequency"; + result = "LKF2HL"; + *argra = 0; + *argdec = 1; + *nargs = 2; + *szargs = 3; + arg[ 0 ] = "RA of source (FK5 J2000, radians)"; + arg[ 1 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 2 ] = "Frequency correction factor"; + break; + + case AST__HLF2LK: + *comment = "Convert from Heliocentric to LSRK frequency"; + result = "HLF2LK"; + *argra = 0; + *argdec = 1; + *nargs = 2; + *szargs = 3; + arg[ 0 ] = "RA of source (FK5 J2000, radians)"; + arg[ 1 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 2 ] = "Frequency correction factor"; + break; + + case AST__LDF2HL: + *comment = "Convert from LSRD to heliocentric frequency"; + result = "LDF2HL"; + *argra = 0; + *argdec = 1; + *nargs = 2; + *szargs = 3; + arg[ 0 ] = "RA of source (FK5 J2000, radians)"; + arg[ 1 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 2 ] = "Frequency correction factor"; + break; + + case AST__HLF2LD: + *comment = "Convert from Heliocentric to LSRD frequency"; + result = "HLF2LD"; + *argra = 0; + *argdec = 1; + *nargs = 2; + *szargs = 3; + arg[ 0 ] = "RA of source (FK5 J2000, radians)"; + arg[ 1 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 2 ] = "Frequency correction factor"; + break; + + case AST__LGF2HL: + *comment = "Convert from Local group to heliocentric frequency"; + result = "LGF2HL"; + *argra = 0; + *argdec = 1; + *nargs = 2; + *szargs = 3; + arg[ 0 ] = "RA of source (FK5 J2000, radians)"; + arg[ 1 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 2 ] = "Frequency correction factor"; + break; + + case AST__HLF2LG: + *comment = "Convert from Heliocentric to local group frequency"; + result = "HLF2LG"; + *argra = 0; + *argdec = 1; + *nargs = 2; + *szargs = 3; + arg[ 0 ] = "RA of source (FK5 J2000, radians)"; + arg[ 1 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 2 ] = "Frequency correction factor"; + break; + + case AST__GLF2HL: + *comment = "Convert from Galactic to heliocentric frequency"; + result = "GLF2HL"; + *argra = 0; + *argdec = 1; + *nargs = 2; + *szargs = 3; + arg[ 0 ] = "RA of source (FK5 J2000, radians)"; + arg[ 1 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 2 ] = "Frequency correction factor"; + break; + + case AST__HLF2GL: + *comment = "Convert from Heliocentric to galactic frequency"; + *argra = 0; + *argdec = 1; + result = "HLF2GL"; + *nargs = 2; + *szargs = 3; + arg[ 0 ] = "RA of source (FK5 J2000, radians)"; + arg[ 1 ] = "DEC of source (FK5 J2000, radians)"; + arg[ 2 ] = "Frequency correction factor"; + break; + + } + +/* Return the result. */ + return result; +} + +static int FrameChange( int cvt_code, int np, double *ra, double *dec, double *freq, + double *args, int forward, int *status ){ +/* +* Name: +* FrameChange + +* Purpose: +* Apply a doppler shift caused by a change of reference frame. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* int FrameChange( int cvt_code, int np, double *ra, double *dec, +* double *freq, double *args, int forward, int *status ) + +* Class Membership: +* SpecMap method. + +* Description: +* This function modifies the supplied frequency values in order to +* apply a doppler shift caused by a change of the observers rest-frame. + +* Parameters: +* cvt_code +* A code indicating the conversion to be applied. If the code does +* not correspond to a change of rest-frame, then the supplied +* frequencies are left unchanged and zero is returned as the +* function value. +* np +* The number of frequency values to transform. +* ra +* Pointer to an array of "np" RA (J2000 FK5) values at which the +* "np" frequencies are observed. These are unchanged on exit. If a +* NULL pointer is supplied, then all frequencies are assumed to be +* observed at the single RA value given by "refra" +* dec +* Pointer to an array of "np" Dec (J2000 FK5) values at which the +* "np" frequencies are observed. These are unchanged on exit. If a +* NULL pointer is supplied, then all frequencies are assumed to be +* observed at the single Dec value given by "refdec" +* freq +* Pointer to an array of "np" frequency values, measured in the +* input rest-frame. These are modified on return to hold the +* corresponding values measured in the output rest-frame. +* args +* Pointer to an array holding the conversion arguments. The number +* of arguments expected depends on the particular conversion being +* used. +* forward +* Should the conversion be applied in the forward or inverse +* direction? Non-zero for forward, zero for inverse. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Non-zero if the supplied conversion code corresponds to a change of +* reference frame. Zoer otherwise (in which case the upplied values +* will not have been changed). + +* Notes: +* - The "args" array contains RA and DEC values which give the "source" +* position (FK5 J2000). If a NULL value is supplied for the "ra" +* parameter, then these args define the position of all the frequency +* values. In addition they also define the direction of motion of +* the "user-defined" rest-frame (see "veluser"). Thus they should still +* be supplied even if "ra" is NULL. + +*/ + +/* Local Variables: */ + FrameDef def; /* Structure holding frame parameters */ + double (* cvtFunc)( double, double, FrameDef *, int * ); /* Pointer to conversion function */ + double *fcorr; /* Pointer to frequency correction factor */ + double *pdec; /* Pointer to next Dec value */ + double *pf; /* Pointer to next frequency value */ + double *pra; /* Pointer to next RA value */ + double factor; /* Frequency correction factor */ + double s; /* Velocity correction (m/s) */ + int i; /* Loop index */ + int result; /* Returned value */ + int sign; /* Sign for velocity correction */ + +/* Check inherited status. */ + if( !astOK ) return 0; + +/* Initialise */ + cvtFunc = NULL; + fcorr = NULL; + sign = 0; + +/* Set the return value to indicate that the supplied conversion code + represents a change of rest-frame. */ + result = 1; + +/* Initialise a structure which stores parameters which define the + transformation. */ + def.obsalt = AST__BAD; + def.obslat = AST__BAD; + def.obslon = AST__BAD; + def.epoch = AST__BAD; + def.refdec = AST__BAD; + def.refra = AST__BAD; + def.veluser = AST__BAD; + def.last = AST__BAD; + def.amprms[ 0 ] = AST__BAD; + def.vuser[ 0 ] = AST__BAD; + def.dvh[ 0 ] = AST__BAD; + def.dvb[ 0 ] = AST__BAD; + +/* Test for each rest-frame code value in turn and assign the appropriate + values. */ + switch ( cvt_code ) { + + case AST__USF2HL: + cvtFunc = UserVel; + def.veluser = args[ 0 ]; + def.refra = args[ 1 ]; + def.refdec = args[ 2 ]; + fcorr = args + 3; + sign = -1; + break; + + case AST__HLF2US: + cvtFunc = UserVel; + def.veluser = args[ 0 ]; + def.refra = args[ 1 ]; + def.refdec = args[ 2 ]; + fcorr = args + 3; + sign = +1; + break; + + case AST__TPF2HL: + cvtFunc = TopoVel; + def.obslon = args[ 0 ]; + def.obslat = args[ 1 ]; + def.obsalt = args[ 2 ]; + def.epoch = args[ 3 ]; + def.refra = args[ 4 ]; + def.refdec = args[ 5 ]; + fcorr = args + 6; + sign = -1; + break; + + case AST__HLF2TP: + cvtFunc = TopoVel; + def.obslon = args[ 0 ]; + def.obslat = args[ 1 ]; + def.obsalt = args[ 2 ]; + def.epoch = args[ 3 ]; + def.refra = args[ 4 ]; + def.refdec = args[ 5 ]; + fcorr = args + 6; + sign = +1; + break; + + case AST__GEF2HL: + cvtFunc = GeoVel; + def.epoch = args[ 0 ]; + def.refra = args[ 1 ]; + def.refdec = args[ 2 ]; + fcorr = args + 3; + sign = -1; + break; + + case AST__HLF2GE: + cvtFunc = GeoVel; + def.epoch = args[ 0 ]; + def.refra = args[ 1 ]; + def.refdec = args[ 2 ]; + fcorr = args + 3; + sign = +1; + break; + + case AST__BYF2HL: + cvtFunc = BaryVel; + def.epoch = args[ 0 ]; + def.refra = args[ 1 ]; + def.refdec = args[ 2 ]; + fcorr = args + 3; + sign = -1; + break; + + case AST__HLF2BY: + cvtFunc = BaryVel; + def.epoch = args[ 0 ]; + def.refra = args[ 1 ]; + def.refdec = args[ 2 ]; + fcorr = args + 3; + sign = +1; + break; + + case AST__LKF2HL: + cvtFunc = LsrkVel; + def.refra = args[ 0 ]; + def.refdec = args[ 1 ]; + fcorr = args + 2; + sign = -1; + break; + + case AST__HLF2LK: + cvtFunc = LsrkVel; + def.refra = args[ 0 ]; + def.refdec = args[ 1 ]; + fcorr = args + 2; + sign = +1; + break; + + case AST__LDF2HL: + cvtFunc = LsrdVel; + def.refra = args[ 0 ]; + def.refdec = args[ 1 ]; + fcorr = args + 2; + sign = -1; + break; + + case AST__HLF2LD: + cvtFunc = LsrdVel; + def.refra = args[ 0 ]; + def.refdec = args[ 1 ]; + fcorr = args + 2; + sign = +1; + break; + + case AST__LGF2HL: + cvtFunc = LgVel; + def.refra = args[ 0 ]; + def.refdec = args[ 1 ]; + fcorr = args + 2; + sign = -1; + break; + + case AST__HLF2LG: + cvtFunc = LgVel; + def.refra = args[ 0 ]; + def.refdec = args[ 1 ]; + fcorr = args + 2; + sign = +1; + break; + + case AST__GLF2HL: + cvtFunc = GalVel; + def.refra = args[ 0 ]; + def.refdec = args[ 1 ]; + fcorr = args + 2; + sign = -1; + break; + + case AST__HLF2GL: + cvtFunc = GalVel; + def.refra = args[ 0 ]; + def.refdec = args[ 1 ]; + fcorr = args + 2; + sign = +1; + break; + +/* If the supplied code does not represent a change of rest-frame, clear + the returned flag. */ + default: + result = 0; + } + +/* Check we have a rest-frame code. */ + if( result ) { + +/* First deal with cases where we have a single source position (given by + refra and refdec). */ + if( !ra ) { + +/* If the frequency correction factor has not been found, find it now. */ + if( *fcorr == AST__BAD ) { + +/* Get the velocity correction. This is the component of the velocity of the + output system, away from the source, as measured in the input system. */ + s = sign*cvtFunc( def.refra, def.refdec, &def, status ); + +/* Find the factor by which to correct supplied frequencies. If the + velocity correction is positive, the output frequency wil be lower than + the input frequency. */ + if( s < AST__C && s > -AST__C ) { + *fcorr = sqrt( ( AST__C - s )/( AST__C + s ) ); + } + } + +/* Correct each supplied frequency. */ + if( *fcorr != AST__BAD && *fcorr != 0.0 ) { + factor = forward ? *fcorr : 1.0 / ( *fcorr ); + pf = freq; + for( i = 0; i < np; i++, pf++ ) { + if( *pf != AST__BAD ) *pf *= factor; + } + +/* Set returned values bad if the velocity correction is un-physical. */ + } else { + pf = freq; + for( i = 0; i < np; i++ ) *(pf++) = AST__BAD; + } + +/* Now deal with cases where each frequency value has its own source + position. */ + } else { + +/* Invert the sign if we are doing a inverse transformation. */ + if( !forward ) sign = -sign; + +/* Loop round each value. */ + pf = freq; + pra = ra; + pdec = dec; + for( i = 0; i < np; i++ ) { + +/* If the ra or dec is bad, store a bad frequency. */ + if( *pra == AST__BAD || *pdec == AST__BAD || *pf == AST__BAD ) { + *pf = AST__BAD; + +/* Otherwise, produce a corrected frequency. */ + } else { + +/* Get the velocity correction. */ + s = sign*cvtFunc( *pra, *pdec, &def, status ); + +/* Correct this frequency, if possible. Otherwise set bad. */ + if( s < AST__C && s > -AST__C ) { + *pf *= sqrt( ( AST__C - s )/( AST__C + s ) ); + } else { + *pf = AST__BAD; + } + } + +/* Move on to the next position. */ + pf++; + pra++; + pdec++; + } + } + } + +/* Return the result. */ + return result; +} + +static double GalVel( double ra, double dec, FrameDef *def, int *status ) { +/* +* Name: +* GalVel + +* Purpose: +* Find the velocity of the galactic centre away from the source. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* double GalVel( double ra, double dec, FrameDef *def, int *status ) + +* Class Membership: +* SpecMap method. + +* Description: +* This function finds the component of the velocity of the galactic +* centre away from a specified source position, in the frame of rest +* of the Sun. + +* Parameters: +* ra +* The RA (rads, FK5 J2000) of the source. +* dec +* The Dec (rads, FK5 J2000) of the source. +* def +* Pointer to a FrameDef structure which holds the parameters which +* define the frame, together with cached intermediate results. +* status +* Pointer to the inherited status variable. + +* Returns: +* The component of the frame's velocity away from the position given by +* "ra" and "dec", in m/s, measured within the Heliographic frame of +* rest. Zero is returned if an error has already occurred. + +*/ + +/* Local Variables: */ + double s1, s2; + +/* Check the global error status. */ + if ( !astOK ) return 0.0; + +/* Get the component away from the source, of the velocity of the sun + relative to the dynamic LSR (in km/s). */ + s1 = (double) palRvlsrd( (float) ra, (float) dec ); + +/* Get the component away from the source, of the velocity of the + dynamic LSR relative to the galactic centre (in km/s). */ + s2 = (double) palRvgalc( (float) ra, (float) dec ); + +/* Return the total velocity of the galactic centre away from the source, + relative to the sun, in m/s. */ + return -1000.0*( s1 + s2 ); +} + +static double GeoVel( double ra, double dec, FrameDef *def, int *status ) { +/* +* Name: +* GeoVel + +* Purpose: +* Find the velocity of the earth away from the source. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* double GeoVel( double ra, double dec, FrameDef *def, int *status ) + +* Class Membership: +* SpecMap method. + +* Description: +* This function finds the component of the velocity of the earth away +* from a specified source position, at a given epoch, in the frame of +* rest of the Sun. + +* Parameters: +* ra +* The RA (rads, FK5 J2000) of the source. +* dec +* The Dec (rads, FK5 J2000) of the source. +* def +* Pointer to a FrameDef structure which holds the parameters which +* define the frame, together with cached intermediate results. +* status +* Pointer to the inherited status variable. + +* Returns: +* The component of the frame's velocity away from the position given by +* "ra" and "dec", in m/s, measured within the Heliographic frame of +* rest. Zero is returned if an error has already occurred. + +*/ + +/* Local Variables: */ + double dpb[ 3 ]; /* Barycentric earth position vector */ + double dph[ 3 ]; /* Heliocentric earth position vector */ + double dvb[ 3 ]; /* Barycentric earth velocity vector */ + double v[ 3 ]; /* Source direction vector */ + +/* Check the global error status. */ + if ( !astOK ) return 0.0; + +/* Get the Cartesian vector towards the source, in the Cartesian FK5 + J2000 system. */ + palDcs2c( ra, dec, v ); + +/* If not already done so, get the Earth/Sun velocity and position vectors in + the same system. Speed is returned in units of AU/s. Store in the supplied + frame definition structure. */ + if( def->dvh[ 0 ] == AST__BAD ) palEvp( def->epoch, 2000.0, dvb, dpb, + def->dvh, dph ); + +/* Return the component away from the source, of the velocity of the earths + centre relative to the sun (in m/s). */ + return -palDvdv( v, def->dvh )*149.597870E9; +} + +void astInitSpecMapVtab_( AstSpecMapVtab *vtab, const char *name, int *status ) { +/* +*+ +* Name: +* astInitSpecMapVtab + +* Purpose: +* Initialise a virtual function table for a SpecMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "specmap.h" +* void astInitSpecMapVtab( AstSpecMapVtab *vtab, const char *name ) + +* Class Membership: +* SpecMap vtab initialiser. + +* Description: +* This function initialises the component of a virtual function +* table which is used by the SpecMap class. + +* Parameters: +* vtab +* Pointer to the virtual function table. The components used by +* all ancestral classes will be initialised if they have not already +* been initialised. +* name +* Pointer to a constant null-terminated character string which contains +* the name of the class to which the virtual function table belongs (it +* is this pointer value that will subsequently be returned by the Object +* astClass function). +*- +*/ + +/* Local Variables: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstObjectVtab *object; /* Pointer to Object component of Vtab */ + AstMappingVtab *mapping; /* Pointer to Mapping component of Vtab */ + +/* Check the local error status. */ + if ( !astOK ) return; + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(NULL); + +/* Initialize the component of the virtual function table used by the + parent class. */ + astInitMappingVtab( (AstMappingVtab *) vtab, name ); + +/* Store a unique "magic" value in the virtual function table. This + will be used (by astIsASpecMap) to determine if an object belongs to + this class. We can conveniently use the address of the (static) + class_check variable to generate this unique value. */ + vtab->id.check = &class_check; + vtab->id.parent = &(((AstMappingVtab *) vtab)->id); + +/* Initialise member function pointers. */ +/* ------------------------------------ */ +/* Store pointers to the member functions (implemented here) that + provide virtual methods for this class. */ + vtab->SpecAdd = SpecAdd; + +/* Save the inherited pointers to methods that will be extended, and + replace them with pointers to the new member functions. */ + object = (AstObjectVtab *) vtab; + mapping = (AstMappingVtab *) vtab; + parent_getobjsize = object->GetObjSize; + object->GetObjSize = GetObjSize; + + parent_transform = mapping->Transform; + mapping->Transform = Transform; + + parent_rate = mapping->Rate; + mapping->Rate = Rate; + +/* Store replacement pointers for methods which will be over-ridden by + new member functions implemented here. */ + object->Equal = Equal; + mapping->MapMerge = MapMerge; + +/* Declare the copy constructor, destructor and class dump + function. */ + astSetCopy( vtab, Copy ); + astSetDelete( vtab, Delete ); + astSetDump( vtab, Dump, "SpecMap", + "Conversion between spectral coordinate systems" ); + +/* If we have just initialised the vtab for the current class, indicate + that the vtab is now initialised, and store a pointer to the class + identifier in the base "object" level of the vtab. */ + if( vtab == &class_vtab ) { + class_init = 1; + astSetVtabClassIdentifier( vtab, &(vtab->id) ); + } +} + +static double LgVel( double ra, double dec, FrameDef *def, int *status ) { +/* +* Name: +* LgVel + +* Purpose: +* Find the velocity of the Local Group away from the source. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* double LgVel( double ra, double dec, FrameDef *def, int *status ) + +* Class Membership: +* SpecMap method. + +* Description: +* This function finds the component of the Local Group velocity away +* from a specified source position, in the frame of rest of the Sun. + +* Parameters: +* ra +* The RA (rads, FK5 J2000) of the source. +* dec +* The Dec (rads, FK5 J2000) of the source. +* def +* Pointer to a FrameDef structure which holds the parameters which +* define the frame, together with cached intermediate results. +* status +* Pointer to the inherited status variable. + +* Returns: +* The component of the frame's velocity away from the position given by +* "ra" and "dec", in m/s, measured within the Heliographic frame of +* rest. Zero is returned if an error has already occurred. + +*/ + +/* Return the component away from the source, of the velocity of the + local group relative to the sun (in m/s). */ + return -1000.0*palRvlg( (float) ra, (float) dec ); +} + +static double LsrdVel( double ra, double dec, FrameDef *def, int *status ) { +/* +* Name: +* LsrdVel + +* Purpose: +* Find the velocity of the Dynamical LSR away from the source. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* double LsrdVel( double ra, double dec, FrameDef *def, int *status ) + +* Class Membership: +* SpecMap method. + +* Description: +* This function finds the component of the velocity of the Dynamical +* LSR away from a specified source position, in the frame of rest of +* the Sun. + +* Parameters: +* ra +* The RA (rads, FK5 J2000) of the source. +* dec +* The Dec (rads, FK5 J2000) of the source. +* def +* Pointer to a FrameDef structure which holds the parameters which +* define the frame, together with cached intermediate results. +* status +* Pointer to the inherited status variable. + +* Returns: +* The component of the frame's velocity away from the position given by +* "ra" and "dec", in m/s, measured within the Heliographic frame of +* rest. Zero is returned if an error has already occurred. + +*/ +/* Check the global error status. */ + if ( !astOK ) return 0.0; + +/* Get the component away from the source, of the velocity of the sun + relative to the dynamical LSR (in m/s). This can also be thought of as the + velocity of the LSR towards the source relative to the sun. Return the + negated value (i.e. velocity of lsrd *away from* the source. */ + return -1000.0*palRvlsrd( (float) ra, (float) dec ); +} + +static double LsrkVel( double ra, double dec, FrameDef *def, int *status ) { +/* +* Name: +* LsrkVel + +* Purpose: +* Find the velocity of the Kinematic LSR away from the source. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* double LsrkVel( double ra, double dec, FrameDef *def, int *status ) + +* Class Membership: +* SpecMap method. + +* Description: +* This function finds the component of the velocity of the Kinematic +* LSR away from a specified source position, in the frame of rest of +* the Sun. + +* Parameters: +* ra +* The RA (rads, FK5 J2000) of the source. +* dec +* The Dec (rads, FK5 J2000) of the source. +* def +* Pointer to a FrameDef structure which holds the parameters which +* define the frame, together with cached intermediate results. +* status +* Pointer to the inherited status variable. + +* Returns: +* The component of the frame's velocity away from the position given by +* "ra" and "dec", in m/s, measured within the Heliographic frame of +* rest. Zero is returned if an error has already occurred. + +*/ +/* Check the global error status. */ + if ( !astOK ) return 0.0; + +/* Get the component away from the source, of the velocity of the sun + relative to the kinematic LSR (in m/s). This can also be thought of as the + velocity of the LSR towards the source relative to the sun. Return the + negated value (i.e. velocity of lsrk *away from* the source. */ + return -1000.0*palRvlsrk( (float) ra, (float) dec ); +} + +static int MapMerge( AstMapping *this, int where, int series, int *nmap, + AstMapping ***map_list, int **invert_list, int *status ) { +/* +* Name: +* MapMerge + +* Purpose: +* Simplify a sequence of Mappings containing a SpecMap. + +* Type: +* Private function. + +* Synopsis: +* #include "mapping.h +* int MapMerge( AstMapping *this, int where, int series, int *nmap, +* AstMapping ***map_list, int **invert_list, int *status ) + +* Class Membership: +* SpecMap method (over-rides the protected astMapMerge method +* inherited from the Mapping class). + +* Description: +* This function attempts to simplify a sequence of Mappings by +* merging a nominated SpecMap in the sequence with its neighbours, +* so as to shorten the sequence if possible. +* +* In many cases, simplification will not be possible and the +* function will return -1 to indicate this, without further +* action. +* +* In most cases of interest, however, this function will either +* attempt to replace the nominated SpecMap with one which it +* considers simpler, or to merge it with the Mappings which +* immediately precede it or follow it in the sequence (both will +* normally be considered). This is sufficient to ensure the +* eventual simplification of most Mapping sequences by repeated +* application of this function. +* +* In some cases, the function may attempt more elaborate +* simplification, involving any number of other Mappings in the +* sequence. It is not restricted in the type or scope of +* simplification it may perform, but will normally only attempt +* elaborate simplification in cases where a more straightforward +* approach is not adequate. + +* Parameters: +* this +* Pointer to the nominated SpecMap which is to be merged with +* its neighbours. This should be a cloned copy of the SpecMap +* pointer contained in the array element "(*map_list)[where]" +* (see below). This pointer will not be annulled, and the +* SpecMap it identifies will not be modified by this function. +* where +* Index in the "*map_list" array (below) at which the pointer +* to the nominated SpecMap resides. +* series +* A non-zero value indicates that the sequence of Mappings to +* be simplified will be applied in series (i.e. one after the +* other), whereas a zero value indicates that they will be +* applied in parallel (i.e. on successive sub-sets of the +* input/output coordinates). +* nmap +* Address of an int which counts the number of Mappings in the +* sequence. On entry this should be set to the initial number +* of Mappings. On exit it will be updated to record the number +* of Mappings remaining after simplification. +* map_list +* Address of a pointer to a dynamically allocated array of +* Mapping pointers (produced, for example, by the astMapList +* method) which identifies the sequence of Mappings. On entry, +* the initial sequence of Mappings to be simplified should be +* supplied. +* +* On exit, the contents of this array will be modified to +* reflect any simplification carried out. Any form of +* simplification may be performed. This may involve any of: (a) +* removing Mappings by annulling any of the pointers supplied, +* (b) replacing them with pointers to new Mappings, (c) +* inserting additional Mappings and (d) changing their order. +* +* The intention is to reduce the number of Mappings in the +* sequence, if possible, and any reduction will be reflected in +* the value of "*nmap" returned. However, simplifications which +* do not reduce the length of the sequence (but improve its +* execution time, for example) may also be performed, and the +* sequence might conceivably increase in length (but normally +* only in order to split up a Mapping into pieces that can be +* more easily merged with their neighbours on subsequent +* invocations of this function). +* +* If Mappings are removed from the sequence, any gaps that +* remain will be closed up, by moving subsequent Mapping +* pointers along in the array, so that vacated elements occur +* at the end. If the sequence increases in length, the array +* will be extended (and its pointer updated) if necessary to +* accommodate any new elements. +* +* Note that any (or all) of the Mapping pointers supplied in +* this array may be annulled by this function, but the Mappings +* to which they refer are not modified in any way (although +* they may, of course, be deleted if the annulled pointer is +* the final one). +* invert_list +* Address of a pointer to a dynamically allocated array which, +* on entry, should contain values to be assigned to the Invert +* attributes of the Mappings identified in the "*map_list" +* array before they are applied (this array might have been +* produced, for example, by the astMapList method). These +* values will be used by this function instead of the actual +* Invert attributes of the Mappings supplied, which are +* ignored. +* +* On exit, the contents of this array will be updated to +* correspond with the possibly modified contents of the +* "*map_list" array. If the Mapping sequence increases in +* length, the "*invert_list" array will be extended (and its +* pointer updated) if necessary to accommodate any new +* elements. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* If simplification was possible, the function returns the index +* in the "map_list" array of the first element which was +* modified. Otherwise, it returns -1 (and makes no changes to the +* arrays supplied). + +* Notes: +* - A value of -1 will be returned if this function is invoked +* with the global error status set, or if it should fail for any +* reason. +*/ + +/* Local Variables: */ + AstMapping *new; /* Pointer to replacement Mapping */ + AstSpecMap *specmap; /* Pointer to SpecMap */ + const char *argdesc[ MAX_ARGS ]; /* Argument descriptions (junk) */ + const char *class; /* Pointer to Mapping class string */ + const char *comment; /* Pointer to comment string (junk) */ + double (*cvtargs)[ MAX_ARGS ]; /* Pointer to argument arrays */ + double tmp; /* Temporary storage */ + int *cvttype; /* Pointer to transformation type codes */ + int *szarg; /* Pointer to argument count array */ + int argdec; /* Index of DEC argument */ + int argra; /* Index of RA argument */ + int done; /* Finished (no further simplification)? */ + int iarg; /* Loop counter for arguments */ + int icvt1; /* Loop initial value */ + int icvt2; /* Loop final value */ + int icvt; /* Loop counter for transformation steps */ + int ikeep; /* Index to store step being kept */ + int imap1; /* Index of first SpecMap to merge */ + int imap2; /* Index of last SpecMap to merge */ + int imap; /* Loop counter for Mappings */ + int inc; /* Increment for transformation step loop */ + int invert; /* SpecMap applied in inverse direction? */ + int istep; /* Loop counter for transformation steps */ + int keep; /* Keep transformation step? */ + int narg; /* Number of user-supplied arguments */ + int ngone; /* Number of Mappings eliminated */ + int nin; /* Numbr of axes for SpecMaps being merged */ + int nstep0; /* Original number of transformation steps */ + int nstep; /* Total number of transformation steps */ + int result; /* Result value to return */ + int simpler; /* Simplification possible? */ + int unit; /* Replacement Mapping is a UnitMap? */ + +/* Initialise. */ + result = -1; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* SpecMaps can only be merged if they are in series (or if there is + only one Mapping present, in which case it makes no difference), so + do nothing if they are not. */ + if ( series || ( *nmap == 1 ) ) { + +/* Save the number of inputs for the SpecMap. */ + nin = astGetNin( this ); + +/* Initialise the number of transformation steps to be merged to equal + the number in the nominated SpecMap. */ + nstep = ( (AstSpecMap *) ( *map_list )[ where ] )->ncvt; + +/* Search adjacent lower-numbered Mappings until one is found which is + not a SpecMap, or is a SpecMap with a different number of axes. Accumulate + the number of transformation steps involved in any SpecMaps found. */ + imap1 = where; + while ( ( imap1 - 1 >= 0 ) && astOK ) { + class = astGetClass( ( *map_list )[ imap1 - 1 ] ); + if ( !astOK || strcmp( class, "SpecMap" ) || + astGetNin( ( *map_list )[ imap1 - 1 ] ) != nin ) break; + nstep += ( (AstSpecMap *) ( *map_list )[ imap1 - 1 ] )->ncvt; + imap1--; + } + +/* Similarly search adjacent higher-numbered Mappings. */ + imap2 = where; + while ( ( imap2 + 1 < *nmap ) && astOK ) { + class = astGetClass( ( *map_list )[ imap2 + 1 ] ); + if ( !astOK || strcmp( class, "SpecMap" ) || + astGetNin( ( *map_list )[ imap2 + 1 ] ) != nin ) break; + nstep += ( (AstSpecMap *) ( *map_list )[ imap2 + 1 ] )->ncvt; + imap2++; + } + +/* Remember the initial number of transformation steps. */ + nstep0 = nstep; + +/* Allocate memory for accumulating a list of all the transformation + steps involved in all the SpecMaps found. */ + cvttype = astMalloc( sizeof( int ) * (size_t) nstep ); + cvtargs = astMalloc( sizeof( double[ MAX_ARGS ] ) * (size_t) nstep ); + szarg = astMalloc( sizeof( int ) * (size_t) nstep ); + +/* Loop to obtain the transformation data for each SpecMap being merged. */ + nstep = 0; + for ( imap = imap1; astOK && ( imap <= imap2 ); imap++ ) { + +/* Obtain a pointer to the SpecMap and note if it is being applied in + its inverse direction. */ + specmap = (AstSpecMap *) ( *map_list )[ imap ]; + invert = ( *invert_list )[ imap ]; + +/* Set up loop limits and an increment to scan the transformation + steps in each SpecMap in either the forward or reverse direction, as + dictated by the associated "invert" value. */ + icvt1 = invert ? specmap->ncvt - 1 : 0; + icvt2 = invert ? -1 : specmap->ncvt; + inc = invert ? -1 : 1; + +/* Loop through each transformation step in the SpecMap. */ + for ( icvt = icvt1; icvt != icvt2; icvt += inc ) { + +/* Store the transformation type code and use "CvtString" to determine + the associated number of arguments. Then store these arguments. */ + cvttype[ nstep ] = specmap->cvttype[ icvt ]; + (void) CvtString( cvttype[ nstep ], &comment, &argra, &argdec, + &narg, szarg + nstep, argdesc, status ); + if ( !astOK ) break; + for ( iarg = 0; iarg < szarg[ nstep ]; iarg++ ) { + cvtargs[ nstep ][ iarg ] = specmap->cvtargs[ icvt ][ iarg ]; + } + +/* If the SpecMap is inverted, we must not only accumulate its + transformation steps in reverse, but also apply them in + reverse. For some steps this means changing arguments, for some it + means changing the transformation type code to a complementary + value, and for others it means both. Define macros to perform each + of the required changes. */ + +/* Macro to exchange a transformation type code for its inverse (and + vice versa). */ +#define SWAP_CODES( code1, code2 ) \ + if ( cvttype[ nstep ] == code1 ) { \ + cvttype[ nstep ] = code2; \ + } else if ( cvttype[ nstep ] == code2 ) { \ + cvttype[ nstep ] = code1; \ + } + +/* Macro to exchange a transformation type code for its inverse (and + vice versa), and reciprocate a specified argument. */ +#define SWAP_CODES2( code1, code2, jarg ) \ + if ( cvttype[ nstep ] == code1 ) { \ + cvttype[ nstep ] = code2; \ + tmp = cvtargs[ nstep ][ jarg ]; \ + if( tmp != AST__BAD && tmp != 0.0 ) { \ + cvtargs[ nstep ][ jarg ] = 1.0/tmp; \ + } else { \ + cvtargs[ nstep ][ jarg ] = AST__BAD; \ + } \ + } else if ( cvttype[ nstep ] == code2 ) { \ + cvttype[ nstep ] = code1; \ + tmp = cvtargs[ nstep ][ jarg ]; \ + if( tmp != AST__BAD && tmp != 0.0 ) { \ + cvtargs[ nstep ][ jarg ] = 1.0/tmp; \ + } else { \ + cvtargs[ nstep ][ jarg ] = AST__BAD; \ + } \ + } + +/* Macro to exchange a transformation type code for its inverse (and + vice versa), and negate a specified argument. */ +#define SWAP_CODES3( code1, code2, jarg ) \ + if ( cvttype[ nstep ] == code1 ) { \ + cvttype[ nstep ] = code2; \ + tmp = cvtargs[ nstep ][ jarg ]; \ + if( tmp != AST__BAD ) { \ + cvtargs[ nstep ][ jarg ] = -tmp; \ + } \ + } else if ( cvttype[ nstep ] == code2 ) { \ + cvttype[ nstep ] = code1; \ + tmp = cvtargs[ nstep ][ jarg ]; \ + if( tmp != AST__BAD ) { \ + cvtargs[ nstep ][ jarg ] = -tmp; \ + } \ + } + +/* Use these macros to apply the changes where needed. */ + if ( invert ) { + +/* Exchange transformation codes for their inverses. */ + SWAP_CODES( AST__FRTOVL, AST__VLTOFR ) + SWAP_CODES( AST__ENTOFR, AST__FRTOEN ) + SWAP_CODES( AST__WNTOFR, AST__FRTOWN ) + SWAP_CODES( AST__WVTOFR, AST__FRTOWV ) + SWAP_CODES( AST__AWTOFR, AST__FRTOAW ) + SWAP_CODES( AST__VRTOVL, AST__VLTOVR ) + SWAP_CODES( AST__VOTOVL, AST__VLTOVO ) + SWAP_CODES( AST__ZOTOVL, AST__VLTOZO ) + SWAP_CODES( AST__BTTOVL, AST__VLTOBT ) + +/* Exchange transformation codes for their inverses, and reciprocate the + frequency correction factor. */ + SWAP_CODES2( AST__TPF2HL, AST__HLF2TP, 6 ) + SWAP_CODES2( AST__USF2HL, AST__HLF2US, 3 ) + SWAP_CODES2( AST__GEF2HL, AST__HLF2GE, 3 ) + SWAP_CODES2( AST__BYF2HL, AST__HLF2BY, 3 ) + SWAP_CODES2( AST__LKF2HL, AST__HLF2LK, 2 ) + SWAP_CODES2( AST__LDF2HL, AST__HLF2LD, 2 ) + SWAP_CODES2( AST__LGF2HL, AST__HLF2LG, 2 ) + SWAP_CODES2( AST__GLF2HL, AST__HLF2GL, 2 ) + + } + +/* Undefine the local macros. */ +#undef SWAP_CODES +#undef SWAP_CODES2 +#undef SWAP_CODES3 + +/* Count the transformation steps. */ + nstep++; + } + } + +/* Loop to simplify the sequence of transformation steps until no + further improvement is possible. */ + done = 0; + while ( astOK && !done ) { + +/* Examine each remaining transformation step in turn. */ + ikeep = -1; + for ( istep = 0; istep < nstep; istep++ ) { + +/* Initially assume we will retain the current step. */ + keep = 1; + +/* The only simplifications for the conversions currently in this class act + to combine adjacent transformation steps, so only apply them while there + are at least 2 steps left. */ + if ( istep < ( nstep - 1 ) ) { + +/* Define a macro to test if two adjacent transformation type codes + have specified values. */ +#define PAIR_CVT( code1, code2 ) \ + ( ( cvttype[ istep ] == code1 ) && \ + ( cvttype[ istep + 1 ] == code2 ) ) + +/* Define a macro to test if two adjacent transformation type codes + have specified values, either way round. */ +#define PAIR_CVT2( code1, code2 ) \ + ( ( PAIR_CVT( code1, code2 ) ) || \ + ( PAIR_CVT( code2, code1 ) ) ) + +/* If a correction is followed by its inverse, and the user-supplied argument + values are unchanged (we do not need to test values stored in the + argument array which were not supplied by the user), we can eliminate them. + First check for conversions which have no user-supplied arguments. */ + if ( PAIR_CVT2( AST__ENTOFR, AST__FRTOEN ) || + PAIR_CVT2( AST__WNTOFR, AST__FRTOWN ) || + PAIR_CVT2( AST__WVTOFR, AST__FRTOWV ) || + PAIR_CVT2( AST__AWTOFR, AST__FRTOAW ) || + PAIR_CVT2( AST__VRTOVL, AST__VLTOVR ) || + PAIR_CVT2( AST__VOTOVL, AST__VLTOVO ) || + PAIR_CVT2( AST__ZOTOVL, AST__VLTOZO ) || + PAIR_CVT2( AST__BTTOVL, AST__VLTOBT ) ) { + istep++; + keep = 0; + +/* Now check for conversions which have a single user-supplied argument. */ + } else if( PAIR_CVT2( AST__FRTOVL, AST__VLTOFR ) && + astEQUAL( cvtargs[ istep ][ 0 ], + cvtargs[ istep + 1 ][ 0 ] ) ) { + istep++; + keep = 0; + +/* Now check for conversions which have two user-supplied arguments. */ + } else if( ( PAIR_CVT2( AST__LKF2HL, AST__HLF2LK ) || + PAIR_CVT2( AST__LDF2HL, AST__HLF2LD ) || + PAIR_CVT2( AST__LGF2HL, AST__HLF2LG ) || + PAIR_CVT2( AST__GLF2HL, AST__HLF2GL ) ) && + astEQUAL( cvtargs[ istep ][ 0 ], + cvtargs[ istep + 1 ][ 0 ] ) && + astEQUAL( cvtargs[ istep ][ 1 ], + cvtargs[ istep + 1 ][ 1 ] ) ) { + istep++; + keep = 0; + +/* Now check for conversions which have three user-supplied arguments. */ + } else if( ( PAIR_CVT2( AST__GEF2HL, AST__HLF2GE ) || + PAIR_CVT2( AST__BYF2HL, AST__HLF2BY ) || + PAIR_CVT2( AST__USF2HL, AST__HLF2US ) ) && + astEQUAL( cvtargs[ istep ][ 0 ], + cvtargs[ istep + 1 ][ 0 ] ) && + astEQUAL( cvtargs[ istep ][ 1 ], + cvtargs[ istep + 1 ][ 1 ] ) && + astEQUAL( cvtargs[ istep ][ 2 ], + cvtargs[ istep + 1 ][ 2 ] ) ) { + istep++; + keep = 0; + +/* Now check for conversions which have six user-supplied arguments (currently + no conversions have four or five user-supplied arguments). */ + } else if( ( PAIR_CVT2( AST__TPF2HL, AST__HLF2TP ) ) && + astEQUAL( cvtargs[ istep ][ 0 ], + cvtargs[ istep + 1 ][ 0 ] ) && + astEQUAL( cvtargs[ istep ][ 1 ], + cvtargs[ istep + 1 ][ 1 ] ) && + astEQUAL( cvtargs[ istep ][ 2 ], + cvtargs[ istep + 1 ][ 2 ] ) && + astEQUAL( cvtargs[ istep ][ 3 ], + cvtargs[ istep + 1 ][ 3 ] ) && + astEQUAL( cvtargs[ istep ][ 4 ], + cvtargs[ istep + 1 ][ 4 ] ) && + astEQUAL( cvtargs[ istep ][ 5 ], + cvtargs[ istep + 1 ][ 5 ] ) ) { + istep++; + keep = 0; + + } + +/* Undefine the local macros. */ +#undef PAIR_CVT +#undef PAIR_CVT2 + } + +/* If the current transformation (possibly modified above) is being + kept, then increment the index that identifies its new location in + the list of transformation steps. */ + if ( keep ) { + ikeep++; + +/* If the new location is different to its current location, copy the + transformation data into the new location. */ + if ( ikeep != istep ) { + cvttype[ ikeep ] = cvttype[ istep ]; + for ( iarg = 0; iarg < szarg[ istep ]; iarg++ ) { + cvtargs[ ikeep ][ iarg ] = cvtargs[ istep ][ iarg ]; + } + szarg[ ikeep ] = szarg[ istep ]; + } + } + } + +/* Note if no simplification was achieved on this iteration (i.e. the + number of transformation steps was not reduced). This is the signal + to quit. */ + done = ( ( ikeep + 1 ) >= nstep ); + +/* Note how many transformation steps now remain. */ + nstep = ikeep + 1; + } + +/* Determine how many Mappings can be eliminated by condensing all + those considered above into a single Mapping. */ + if ( astOK ) { + ngone = imap2 - imap1; + +/* Determine if the replacement Mapping can be a UnitMap (a null + Mapping). This will only be the case if all the transformation + steps were eliminated above. */ + unit = ( nstep == 0 ); + +/* Determine if simplification is possible. This will be the case if + (a) Mappings were eliminated ("ngone" is non-zero), or (b) the + number of transformation steps was reduced, or (c) the SpecMap(s) + can be replaced by a UnitMap, or (d) if there was initially only + one SpecMap present, its invert flag was set (this flag will always + be cleared in the replacement Mapping). */ + simpler = ngone || ( nstep < nstep0 ) || unit || + ( *invert_list )[ where ]; + +/* Do nothing more unless simplification is possible. */ + if ( simpler ) { + +/* If the replacement Mapping is a UnitMap, then create it. */ + if ( unit ) { + new = (AstMapping *) + astUnitMap( astGetNin( ( *map_list )[ where ] ), "", status ); + +/* Otherwise, create a replacement SpecMap and add each of the + remaining transformation steps to it. */ + } else { + new = (AstMapping *) astSpecMap( nin, 0, "", status ); + for ( istep = 0; istep < nstep; istep++ ) { + AddSpecCvt( (AstSpecMap *) new, cvttype[ istep ], + cvtargs[ istep ], status ); + } + } + +/* Annul the pointers to the Mappings being eliminated. */ + if ( astOK ) { + for ( imap = imap1; imap <= imap2; imap++ ) { + ( *map_list )[ imap ] = astAnnul( ( *map_list )[ imap ] ); + } + +/* Insert the pointer and invert value for the new Mapping. */ + ( *map_list )[ imap1 ] = new; + ( *invert_list )[ imap1 ] = 0; + +/* Move any subsequent Mapping information down to close the gap. */ + for ( imap = imap2 + 1; imap < *nmap; imap++ ) { + ( *map_list )[ imap - ngone ] = ( *map_list )[ imap ]; + ( *invert_list )[ imap - ngone ] = ( *invert_list )[ imap ]; + } + +/* Blank out any information remaining at the end of the arrays. */ + for ( imap = ( *nmap - ngone ); imap < *nmap; imap++ ) { + ( *map_list )[ imap ] = NULL; + ( *invert_list )[ imap ] = 0; + } + +/* Decrement the Mapping count and return the index of the first + Mapping which was eliminated. */ + ( *nmap ) -= ngone; + result = imap1; + +/* If an error occurred, annul the new Mapping pointer. */ + } else { + new = astAnnul( new ); + } + } + } + +/* Free the memory used for the transformation steps. */ + cvttype = astFree( cvttype ); + cvtargs = astFree( cvtargs ); + szarg = astFree( szarg ); + } + +/* If an error occurred, clear the returned value. */ + if ( !astOK ) result = -1; + +/* Return the result. */ + return result; +} + +static double Rate( AstMapping *this, double *at, int ax1, int ax2, int *status ){ +/* +* Name: +* Rate + +* Purpose: +* Calculate the rate of change of a Mapping output. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* result = Rate( AstMapping *this, double *at, int ax1, int ax2, int *status ) + +* Class Membership: +* SpecMap member function (overrides the astRate method inherited +* from the Mapping class ). + +* Description: +* This function returns the rate of change of a specified output of +* the supplied Mapping with respect to a specified input, at a +* specified input position. + +* Parameters: +* this +* Pointer to the Mapping to be applied. +* at +* The address of an array holding the axis values at the position +* at which the rate of change is to be evaluated. The number of +* elements in this array should equal the number of inputs to the +* Mapping. +* ax1 +* The index of the Mapping output for which the rate of change is to +* be found (output numbering starts at 0 for the first output). +* ax2 +* The index of the Mapping input which is to be varied in order to +* find the rate of change (input numbering starts at 0 for the first +* input). +* status +* Pointer to the inherited status variable. + +* Returned Value: +* The rate of change of Mapping output "ax1" with respect to input +* "ax2", evaluated at "at", or AST__BAD if the value cannot be +* calculated. + +* Implementation Deficiencies: +* The initial version of this implementation only deals with +* frequency->wavelength conversions. This is because the slowness of +* the numerical differentiation implemented by the astRate method in +* the parent Mapping class is cripples conversion between SpecFluxFrames. +* Such conversions only rely on rate of change of wavelength with +* respect to frequency. This implementation should be extended when +* needed. + +*/ + +/* Local Variables: */ + AstSpecMap *map; + double result; + int cvt; + +/* Check inherited status */ + if( !astOK ) return AST__BAD; + +/* Get a pointer to the SpecMap structure. */ + map = (AstSpecMap *) this; + +/* Return 1.0 if the SpecMap has no conversions. */ + if( map->ncvt == 0 ) return 1.0; + +/* Store the type of the first conversion.*/ + cvt = map->cvttype[ 0 ]; + +/* If this is a 3D SpecMap or if it has more than one component, or if + that conversion is not between frequency and wavelength, use the + astRate method inherited form the parent Mapping class. */ + if( astGetNin( map ) != 1 || map->ncvt != 1 || + ( cvt != AST__WVTOFR && cvt != AST__FRTOWV ) ) { + result = (*parent_rate)( this, at, ax1, ax2, status ); + +/* Otherwise, evaluate the known analytical expressions for the rate of + change of frequency with respect to wavelength or wavelength with + respect to frequency. */ + } else { + result = ( *at != AST__BAD ) ? -AST__C/((*at)*(*at)) : AST__BAD; + } + +/* Return the result. */ + return result; +} + +static double Refrac( double wavelen, int *status ){ +/* +* Name: +* Refrac + +* Purpose: +* Returns the refractive index of dry air at a given wavelength. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* double Refrac( double wavelen, int *status ) + +* Class Membership: +* SpecMap method. + +* Description: +* This function returns the refractive index of dry air at standard +* temperature and pressure, at a given wavelength. The formula is +* taken from the paper "Representation of Spectral Coordinates in FITS" +* (Greisen et al). + +* Parameters: +* wavelen +* The wavelength, in metres. This should be the air wavelength, +* but supplying the vacuum wavelength will make no significant +* difference. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* The refractive index. A value of 1.0 is returned if an error +* occurs, or has already occurred. + +*/ + +/* Local Variables: */ + double w2; /* Wavenumber squared */ + +/* Check the global error status. */ + if ( !astOK || wavelen == 0.0 ) return 1.0; + +/* Find the squared wave number in units of "(per um)**2". */ + w2 = 1.0E-12/( wavelen * wavelen ); + +/* Apply the rest of the algorithm as described in the FITS WCS + paper III. */ + return 1.0 + 1.0E-6*( 287.6155 + 1.62887*w2 + 0.01360*w2*w2 ); +} + +static double Rverot( double phi, double h, double ra, double da, + double st, int *status ) { +/* +* Name: +* Rverot + +* Purpose: +* Find the velocity component in a given direction due to Earth rotation. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* double Rverot( double phi, double h, double ra, double da, +* double st, int *status ) + +* Class Membership: +* SpecMap method. + +* Description: +* This function is like slaRverot, except that it takes account of the +* observers height (h), and does all calculations in double precision. + +* Parameters: +* phi +* The geodetic latitude of the observer (radians, IAU 1976). +* h +* The geodetic height above the reference spheroid of the observer +* (metres, IAU 1976). +* ra +* The geocentric apparent RA (rads) of the source. +* da +* The geocentric apparent Dec (rads) of the source. +* st +* The local apparent sidereal time (radians). +* status +* Pointer to the inherited status variable. + +* Returns: +* The component of the Earth rotation in direction [RA,DA] (km/s). +* The result is positive when the observer is receding from the +* given point on the sky. Zero is returned if an error has already +* occurred. + +*/ + +/* Local Variables: */ + double pv[ 6 ]; /* Observer position and velocity */ + double v[ 3 ]; /* Source direction vector */ + +/* Check the global error status. */ + if ( !astOK ) return 0.0; + +/* Get the Cartesian coordinates of the unit vector pointing towards the + given sky position. */ + palDcs2c( ra, da, v ); + +/* Get velocity and position of the observer. */ + palPvobs( phi, h, st, pv ); + +/* Return the component of the observer's velocity away from the sky + position, and convert from AU/s to km/s. */ + return -palDvdv( v, pv + 3 )*149.597870E6; +} + +static void SpecAdd( AstSpecMap *this, const char *cvt, const double args[], int *status ) { +/* +*++ +* Name: +c astSpecAdd +f AST_SPECADD + +* Purpose: +* Add a spectral coordinate conversion to a SpecMap. + +* Type: +* Public virtual function. + +* Synopsis: +c #include "specmap.h" +c void astSpecAdd( AstSpecMap *this, const char *cvt, const double args[] ) +f CALL AST_SPECADD( THIS, CVT, ARGS, STATUS ) + +* Class Membership: +* SpecMap method. + +* Description: +c This function adds one of the standard spectral coordinate +f This routine adds one of the standard spectral coordinate +* system conversions listed below to an existing SpecMap. +* +c When a SpecMap is first created (using astSpecMap), it simply +f When a SpecMap is first created (using AST_SPECMAP), it simply +c performs a unit (null) Mapping. By using astSpecAdd (repeatedly +f performs a unit (null) Mapping. By using AST_SPECADD (repeatedly +* if necessary), one or more coordinate conversion steps may then +* be added, which the SpecMap will perform in sequence. This allows +* multi-step conversions between a variety of spectral coordinate +* systems to be assembled out of the building blocks provided by +* this class. +* +* Normally, if a SpecMap's Invert attribute is zero (the default), +* then its forward transformation is performed by carrying out +* each of the individual coordinate conversions specified by +c astSpecAdd in the order given (i.e. with the most recently added +f AST_SPECADD in the order given (i.e. with the most recently added +* conversion applied last). +* +* This order is reversed if the SpecMap's Invert attribute is +* non-zero (or if the inverse transformation is requested by any +* other means) and each individual coordinate conversion is also +* replaced by its own inverse. This process inverts the overall +* effect of the SpecMap. In this case, the first conversion to be +* applied would be the inverse of the one most recently added. + +* Parameters: +c this +f THIS = INTEGER (Given) +* Pointer to the SpecMap. +c cvt +f CVT = CHARACTER * ( * ) (Given) +c Pointer to a null-terminated string which identifies the +f A character string which identifies the +* spectral coordinate conversion to be added to the +* SpecMap. See the "Available Conversions" section for details of +* those available. +c args +f ARGS( * ) = DOUBLE PRECISION (Given) +* An array containing argument values for the spectral +* coordinate conversion. The number of arguments required, and +* hence the number of array elements used, depends on the +* conversion specified (see the "Available Conversions" +* section). This array is ignored +c and a NULL pointer may be supplied +* if no arguments are needed. +f STATUS = INTEGER (Given and Returned) +f The global status. + +* Notes: +* - When assembling a multi-stage conversion, it can sometimes be +* difficult to determine the most economical conversion path. For +* example, when converting between reference frames, converting first +* to the heliographic reference frame as an intermediate stage is often +* sensible in formulating the problem, but may introduce unnecessary +* extra conversion steps. A solution to this is to include all the steps +* which are (logically) necessary, but then to use +c astSimplify to simplify the resulting +f AST_SIMPLIFY to simplify the resulting +* SpecMap. The simplification process will eliminate any steps +* which turn out not to be needed. +c - This function does not check to ensure that the sequence of +f - This routine does not check to ensure that the sequence of +* coordinate conversions added to a SpecMap is physically +* meaningful. + +* Available Conversions: +* The following strings (which are case-insensitive) may be supplied +c via the "cvt" parameter to indicate which spectral coordinate +f via the CVT argument to indicate which spectral coordinate +* conversion is to be added to the SpecMap. Where arguments are needed by +* the conversion, they are listed in parentheses. Values for +c these arguments should be given, via the "args" array, in the +f these arguments should be given, via the ARGS array, in the +* order indicated. Units and argument names are described at the end of +* the list of conversions. + +* - "FRTOVL" (RF): Convert frequency to relativistic velocity. +* - "VLTOFR" (RF): Convert relativistic velocity to Frequency. +* - "ENTOFR": Convert energy to frequency. +* - "FRTOEN": Convert frequency to energy. +* - "WNTOFR": Convert wave number to frequency. +* - "FRTOWN": Convert frequency to wave number. +* - "WVTOFR": Convert wavelength (vacuum) to frequency. +* - "FRTOWV": Convert frequency to wavelength (vacuum). +* - "AWTOFR": Convert wavelength (air) to frequency. +* - "FRTOAW": Convert frequency to wavelength (air). +* - "VRTOVL": Convert radio to relativistic velocity. +* - "VLTOVR": Convert relativistic to radio velocity. +* - "VOTOVL": Convert optical to relativistic velocity. +* - "VLTOVO": Convert relativistic to optical velocity. +* - "ZOTOVL": Convert redshift to relativistic velocity. +* - "VLTOZO": Convert relativistic velocity to redshift. +* - "BTTOVL": Convert beta factor to relativistic velocity. +* - "VLTOBT": Convert relativistic velocity to beta factor. +* - "USF2HL" (VOFF,RA,DEC): Convert frequency from a user-defined +* reference frame to heliocentric. +* - "HLF2US" (VOFF,RA,DEC): Convert frequency from heliocentric +* reference frame to user-defined. +* - "TPF2HL" (OBSLON,OBSLAT,OBSALT,EPOCH,RA,DEC): Convert frequency from +* topocentric reference frame to heliocentric. +* - "HLF2TP" (OBSLON,OBSLAT,OBSALT,EPOCH,RA,DEC): Convert frequency from +* heliocentric reference frame to topocentric. +* - "GEF2HL" (EPOCH,RA,DEC): Convert frequency from geocentric +* reference frame to heliocentric. +* - "HLF2GE" (EPOCH,RA,DEC): Convert frequency from +* heliocentric reference frame to geocentric. +* - "BYF2HL" (EPOCH,RA,DEC): Convert frequency from +* barycentric reference frame to heliocentric. +* - "HLF2BY" (EPOCH,RA,DEC): Convert frequency from +* heliocentric reference frame to barycentric. +* - "LKF2HL" (RA,DEC): Convert frequency from kinematic LSR +* reference frame to heliocentric. +* - "HLF2LK" (RA,DEC): Convert frequency from heliocentric +* reference frame to kinematic LSR. +* - "LDF2HL" (RA,DEC): Convert frequency from dynamical LSR +* reference frame to heliocentric. +* - "HLF2LD" (RA,DEC): Convert frequency from heliocentric +* reference frame to dynamical LSR. +* - "LGF2HL" (RA,DEC): Convert frequency from local group +* reference frame to heliocentric. +* - "HLF2LG" (RA,DEC): Convert frequency from heliocentric +* reference frame to local group. +* - "GLF2HL" (RA,DEC): Convert frequency from galactic +* reference frame to heliocentric. +* - "HLF2GL" (RA,DEC): Convert frequency from heliocentric +* reference frame to galactic. + +* The units for the values processed by the above conversions are as +* follows: +* +* - all velocities: metres per second (positive if the source receeds from +* the observer). +* - frequency: Hertz. +* - all wavelengths: metres. +* - energy: Joules. +* - wave number: cycles per metre. +* +* The arguments used in the above conversions are as follows: +* +* - RF: Rest frequency (Hz). +* - OBSALT: Geodetic altitude of observer (IAU 1975, metres). +* - OBSLAT: Geodetic latitude of observer (IAU 1975, radians). +* - OBSLON: Longitude of observer (radians - positive eastwards). +* - EPOCH: Epoch of observation (UT1 expressed as a Modified Julian Date). +* - RA: Right Ascension of source (radians, FK5 J2000). +* - DEC: Declination of source (radians, FK5 J2000). +* - VOFF: Velocity of the user-defined reference frame, towards the +* position given by RA and DEC, measured in the heliocentric +* reference frame. +* +* If the SpecMap is 3-dimensional, source positions are provided by the +* values supplied to inputs 2 and 3 of the SpecMap (which are simply +* copied to outputs 2 and 3). Note, usable values are still required +* for the RA and DEC arguments in order to define the "user-defined" +* reference frame used by USF2HL and HLF2US. However, AST__BAD can be +* supplied for RA and DEC if the user-defined reference frame is not +* required. +* +*-- +*/ + +/* Local Variables: */ + int cvttype; /* Conversion type code */ + +/* Check the inherited status. */ + if ( !astOK ) return; + +/* Validate the type string supplied and obtain the equivalent + conversion type code. */ + cvttype = CvtCode( cvt, status ); + +/* If the string was not recognised, then report an error. */ + if ( astOK && ( cvttype == AST__SPEC_NULL ) ) { + astError( AST__SPCIN, + "%s(%s): Invalid SpecMap spectral coordinate " + "conversion type \"%s\".", status, "astAddSpec", astGetClass( this ), cvt ); + } + +/* Add the new conversion to the SpecMap. */ + AddSpecCvt( this, cvttype, args, status ); +} + +static int SystemChange( int cvt_code, int np, double *values, double *args, + int forward, int *status ){ +/* +* Name: +* SystemChange + +* Purpose: +* Change values between two spectral systems. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* int SystemChange( int cvt_code, int np, double *values, double *args, +* int forward, int *status ) + +* Class Membership: +* SpecMap method. + +* Description: +* This function modifies the supplied values in order to change the +* spectral co-ordinate system (frequency, wavelength, etc) to which +* they refer. + +* Parameters: +* cvt_code +* A code indicating the conversion to be applied. If the code does +* not correspond to a change of system, then the supplied values +* are left unchanged and zero is returned as the function value. +* np +* The number of frequency values to transform. +* values +* Pointer to an array of "np" spectral values. These are modified on +* return to hold the corresponding values measured in the output +* system. +* args +* Pointer to an array holding the conversion arguments. The number +* of arguments expected depends on the particular conversion being +* used. +* forward +* Should the conversion be applied in the forward or inverse +* direction? Non-zero for forward, zero for inverse. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Non-zero if the supplied conversion code corresponds to a change of +* system. Zero otherwise (in which case the upplied values will not +* have been changed). + +*/ + +/* Local Variables: */ + double *pv; /* Pointer to next value */ + double d; /* Intermediate value */ + double f2; /* Squared frequency */ + double temp; /* Intermediate value */ + int i; /* Loop index */ + int iter; /* Iteration count */ + int result; /* Returned value */ + +/* Check inherited status. */ + if( !astOK ) return 0; + +/* Set the return value to indicate that the supplied conversion code + represents a change of system. */ + result = 1; + +/* Test for each code value in turn and assign the appropriate values. */ + switch ( cvt_code ) { + +/* Frequency to relativistic velocity. */ + case AST__FRTOVL: + if( forward ) { + if( args[ 0 ] != AST__BAD ) { + temp = args[ 0 ] * args[ 0 ]; + pv = values - 1; + for( i = 0; i < np; i++ ){ + pv++; + if( *pv != AST__BAD ) { + f2 = ( *pv ) * ( *pv ); + d = temp + f2; + if( d > 0.0 ) { + *pv = AST__C*( ( temp - f2 )/d ); + if( !astISFINITE( *pv ) ) *pv = AST__BAD; + } else { + *pv = AST__BAD; + } + } + } + } else { + pv = values; + for( i = 0; i < np; i++ ) *( pv++ ) = AST__BAD; + } + } else { + SystemChange( AST__VLTOFR, np, values, args, 1, status ); + } + break; + +/* Relativistic velocity to frequency. */ + case AST__VLTOFR: + if( forward ) { + if( args[ 0 ] != AST__BAD ) { + temp = args[ 0 ]; + pv = values - 1; + for( i = 0; i < np; i++ ){ + pv++; + if( *pv != AST__BAD ) { + d = AST__C + ( *pv ); + if( d != 0.0 ) { + d = ( AST__C - ( *pv ) )/d; + if( d >= 0.0 ) { + *pv = temp*sqrt( d ); + if( !astISFINITE( *pv ) ) *pv = AST__BAD; + } else { + *pv = AST__BAD; + } + } else { + *pv = AST__BAD; + } + } + } + } else { + pv = values; + for( i = 0; i < np; i++ ) *( pv++ ) = AST__BAD; + } + } else { + SystemChange( AST__FRTOVL, np, values, args, 1, status ); + } + break; + +/* Energy to frequency */ + case AST__ENTOFR: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD ) { + *pv /= AST__H; + } + } + } else { + SystemChange( AST__FRTOEN, np, values, args, 1, status ); + } + break; + +/* Frequency to energy */ + case AST__FRTOEN: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD ) { + *pv *= AST__H; + } + } + } else { + SystemChange( AST__ENTOFR, np, values, args, 1, status ); + } + break; + +/* Wave number to frequency */ + case AST__WNTOFR: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD ) { + *pv *= AST__C; + } + } + } else { + SystemChange( AST__FRTOWN, np, values, args, 1, status ); + } + break; + +/* Wave number to frequency */ + case AST__FRTOWN: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD ) { + *pv /= AST__C; + } + } + } else { + SystemChange( AST__WNTOFR, np, values, args, 1, status ); + } + break; + +/* Wavelength to frequency */ + case AST__WVTOFR: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD && *pv != 0.0 ) { + *pv = AST__C/( *pv ); + if( !astISFINITE( *pv ) ) *pv = AST__BAD; + } else { + *pv = AST__BAD; + } + } + } else { + SystemChange( AST__FRTOWV, np, values, args, 1, status ); + } + break; + +/* Frequency to wavelength. */ + case AST__FRTOWV: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD && *pv != 0.0 ) { + *pv = AST__C/( *pv ); + if( !astISFINITE( *pv ) ) *pv = AST__BAD; + } else { + *pv = AST__BAD; + } + } + } else { + SystemChange( AST__WVTOFR, np, values, args, 1, status ); + } + break; + +/* Wavelength in air to frequency. */ + case AST__AWTOFR: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD && *pv != 0.0 ) { + *pv = AST__C/( ( *pv )*Refrac( *pv, status ) ); + if( !astISFINITE( *pv ) ) *pv = AST__BAD; + } else { + *pv = AST__BAD; + } + } + } else { + SystemChange( AST__FRTOAW, np, values, args, 1, status ); + } + break; + +/* Frequency to wavelength in air. */ + case AST__FRTOAW: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD && *pv != 0.0 ) { + +/* Form the vacuum wavelength. */ + temp = AST__C/( *pv ); + +/* The refractive index function "Refrac" requires the wavelength in air + as its parameter. Initially assume that the wavelength in air is equal + to the vacuum wavelength to get he first estimate of the wavelength in + air. Then use this estimate to get a better refractive index in order to + form a better estimate of the air wavelength, etc. Iterate in this way a + few times. */ + *pv = temp; + for( iter = 0; iter < 3; iter++ ) { + *pv = temp/Refrac( *pv, status ); + if( !astISFINITE( *pv ) ) { + *pv = AST__BAD; + break; + } + } + + } else { + *pv = AST__BAD; + } + } + } else { + SystemChange( AST__AWTOFR, np, values, args, 1, status ); + } + break; + +/* Radio velocity to relativistic velocity */ + case AST__VRTOVL: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD ) { + temp = 1.0 - ( *pv )/AST__C; + temp *= temp; + *pv = AST__C*( 1.0 - temp )/( 1.0 + temp ); + if( !astISFINITE( *pv ) ) *pv = AST__BAD; + } + } + } else { + SystemChange( AST__VLTOVR, np, values, args, 1, status ); + } + break; + +/* Relativistic velocity to radio velocity. */ + case AST__VLTOVR: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD ) { + temp = AST__C + ( *pv ); + if( temp != 0.0 ) { + temp = (AST__C - *pv )/temp; + if( temp >= 0.0 ) { + *pv = AST__C*( 1.0 - sqrt( temp ) ); + if( !astISFINITE( *pv ) ) *pv = AST__BAD; + } else { + *pv = AST__BAD; + } + } else { + *pv = AST__BAD; + } + } + } + } else { + SystemChange( AST__VRTOVL, np, values, args, 1, status ); + } + break; + +/* Optical velocity to relativistic velocity */ + case AST__VOTOVL: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD ) { + temp = 1.0 + ( *pv )/AST__C; + temp *= temp; + *pv = AST__C*( temp - 1.0 )/( temp + 1.0 ); + if( !astISFINITE( *pv ) ) *pv = AST__BAD; + } + } + } else { + SystemChange( AST__VLTOVO, np, values, args, 1, status ); + } + break; + +/* Relativistic velocity to optical velocity. */ + case AST__VLTOVO: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD ) { + temp = AST__C - *pv; + if( temp != 0.0 ) { + temp = (AST__C + *pv )/temp; + if( temp >= 0.0 ) { + *pv = AST__C*( sqrt( temp ) - 1.0 ); + if( !astISFINITE( *pv ) ) *pv = AST__BAD; + } else { + *pv = AST__BAD; + } + } else { + *pv = AST__BAD; + } + } + } + } else { + SystemChange( AST__VOTOVL, np, values, args, 1, status ); + } + break; + +/* Redshift to relativistic velocity */ + case AST__ZOTOVL: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD ) { + temp = 1.0 + ( *pv ); + temp *= temp; + *pv = AST__C*( temp - 1.0 )/( temp + 1.0 ); + if( !astISFINITE( *pv ) ) *pv = AST__BAD; + } + } + } else { + SystemChange( AST__VLTOZO, np, values, args, 1, status ); + } + break; + +/* Relativistic velocity to redshift. */ + case AST__VLTOZO: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD ) { + temp = AST__C - *pv; + if( temp != 0.0 ) { + temp = (AST__C + *pv )/temp; + if( temp >= 0.0 ) { + *pv = sqrt( temp ) - 1.0; + if( !astISFINITE( *pv ) ) *pv = AST__BAD; + } else { + *pv = AST__BAD; + } + } else { + *pv = AST__BAD; + } + } + } + } else { + SystemChange( AST__ZOTOVL, np, values, args, 1, status ); + } + break; + +/* Beta factor to relativistic velocity */ + case AST__BTTOVL: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD ) { + *pv *= AST__C; + } + } + } else { + SystemChange( AST__VLTOBT, np, values, args, 1, status ); + } + break; + +/* Relativistic velocity to beta factor. */ + case AST__VLTOBT: + if( forward ) { + pv = values - 1; + for( i = 0; i < np; i++ ) { + pv++; + if( *pv != AST__BAD ) { + *pv /= AST__C; + } + } + } else { + SystemChange( AST__BTTOVL, np, values, args, 1, status ); + } + break; + +/* If the supplied code does not represent a change of system, clear + the returned flag. */ + default: + result = 0; + } + +/* Return the result. */ + return result; +} + +static double TopoVel( double ra, double dec, FrameDef *def, int *status ) { +/* +* Name: +* TopoVel + +* Purpose: +* Find the velocity of the observer away from the source. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* double TopoVel( double ra, double dec, FrameDef *def, int *status ) + +* Class Membership: +* SpecMap method. + +* Description: +* This function finds the component of the velocity of the observer away +* from a specified source position, at a given epoch, in the frame of +* rest of the Sun. + +* Parameters: +* ra +* The RA (rads, FK5 J2000) of the source. +* dec +* The Dec (rads, FK5 J2000) of the source. +* def +* Pointer to a FrameDef structure which holds the parameters which +* define the frame, together with cached intermediate results. +* status +* Pointer to the inherited status variable. + +* Returns: +* The component of the frame's velocity away from the position given by +* "ra" and "dec", in m/s, measured within the Heliographic frame of +* rest. Zero is returned if an error has already occurred. + +*/ + +/* Local Variables: */ + double deca; /* Apparent DEC */ + double raa; /* Apparent RA */ + double vobs; /* Velocity of observer relative to earth */ + double vearth; /* Velocity of earth realtive to sun */ + +/* Check the global error status. */ + if ( !astOK ) return 0.0; + +/* If not already done so, get the parameters defining the transformation + of mean ra and dec to apparent ra and dec, and store in the supplied frame + definition structure. */ + if( def->amprms[ 0 ] == AST__BAD ) palMappa( 2000.0, def->epoch, + def->amprms ); + +/* Convert the source position from mean ra and dec to apparent ra and dec. */ + palMapqkz( ra, dec, def->amprms, &raa, &deca ); + +/* If not already done so, get the local apparent siderial time (in radians) + and store in the supplied frame definition structure. */ + if( def->last == AST__BAD ) def->last = palGmst( def->epoch ) + + palEqeqx( def->epoch ) + + def->obslon; + +/* Get the component away from the source, of the velocity of the observer + relative to the centre of the earth (in m/s). */ + vobs = 1000.0*Rverot( def->obslat, def->obsalt, raa, deca, def->last, + status ); + +/* Get the component away from the source, of the velocity of the earth's + centre relative to the Sun, in m/s. */ + vearth = GeoVel( ra, dec, def, status ); + +/* Return the total velocity of the observer away from the source in the + frame of the sun. */ + return vobs + vearth; +} + +static AstPointSet *Transform( AstMapping *this, AstPointSet *in, + int forward, AstPointSet *out, int *status ) { +/* +* Name: +* Transform + +* Purpose: +* Apply a SpecMap to transform a set of points. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* AstPointSet *Transform( AstMapping *this, AstPointSet *in, +* int forward, AstPointSet *out, int *status ) + +* Class Membership: +* SpecMap member function (over-rides the astTransform method inherited +* from the Mapping class). + +* Description: +* This function takes a SpecMap and a set of points encapsulated +* in a PointSet and transforms the points so as to perform the +* sequence of spectral coordinate conversions specified by +* previous invocations of astSpecAdd. + +* Parameters: +* this +* Pointer to the SpecMap. +* in +* Pointer to the PointSet holding the input coordinate data. +* forward +* A non-zero value indicates that the forward coordinate transformation +* should be applied, while a zero value requests the inverse +* transformation. +* out +* Pointer to a PointSet which will hold the transformed (output) +* coordinate values. A NULL value may also be given, in which case a +* new PointSet will be created by this function. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Pointer to the output (possibly new) PointSet. + +* Notes: +* - A null pointer will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +* - The number of coordinate values per point in the input PointSet must +* match the number of coordinates for the SpecMap being applied. +* - If an output PointSet is supplied, it must have space for sufficient +* number of points and coordinate values per point to accommodate the +* result. Any excess space will be ignored. +*/ + +/* Local Variables: */ + AstPointSet *result; /* Pointer to output PointSet */ + AstSpecMap *map; /* Pointer to SpecMap to be applied */ + double **ptr_in; /* Pointer to input coordinate data */ + double **ptr_out; /* Pointer to output coordinate data */ + double *spec; /* Pointer to output spectral axis value array */ + double *alpha; /* Pointer to output RA axis value array */ + double *beta; /* Pointer to output DEC axis value array */ + int cvt; /* Loop counter for conversions */ + int end; /* Termination index for conversion loop */ + int inc; /* Increment for conversion loop */ + int map3d; /* Is the SpecMap 3-dimensional? */ + int ncoord_in; /* Number of coordinates per input point */ + int npoint; /* Number of points */ + int start; /* Starting index for conversion loop */ + +/* Check the global error status. */ + if ( !astOK ) return NULL; + +/* Obtain a pointer to the SpecMap. */ + map = (AstSpecMap *) this; + +/* Apply the parent mapping using the stored pointer to the Transform member + function inherited from the parent Mapping class. This function validates + all arguments and generates an output PointSet if necessary, but does not + actually transform any coordinate values. */ + result = (*parent_transform)( this, in, forward, out, status ); + +/* We will now extend the parent astTransform method by performing the + coordinate conversions needed to generate the output coordinate values. */ + +/* Determine the numbers of points and coordinates per point from the input + PointSet and obtain pointers for accessing the input and output coordinate + values. */ + ncoord_in = astGetNcoord( in ); + npoint = astGetNpoint( in ); + ptr_in = astGetPoints( in ); + ptr_out = astGetPoints( result ); + +/* Determine whether to apply the forward or inverse transformation, according + to the direction specified and whether the mapping has been inverted. */ + if ( astGetInvert( this ) ) forward = !forward; + +/* Transform the coordinate values. */ +/* -------------------------------- */ +/* Use "spec" as a synonym for the array of spectral axis values stored in + the output PointSet. */ + if ( astOK ) { + spec = ptr_out[ 0 ]; + +/* If this is a 3D SpecMap use "alpha" as a synonym for the array of RA axis + values and "beta" as a synonym for the array of DEC axis values stored + in the output PointSet. */ + map3d = ( ncoord_in == 3 ); + if( map3d ) { + alpha = ptr_out[ 1 ]; + beta = ptr_out[ 2 ]; + } else { + alpha = NULL; + beta = NULL; + } + +/* Initialise the output coordinate values by copying the input ones. */ + (void) memcpy( spec, ptr_in[ 0 ], sizeof( double ) * (size_t) npoint ); + if( map3d ) { + (void) memcpy( alpha, ptr_in[ 1 ], sizeof( double ) * (size_t) npoint ); + (void) memcpy( beta, ptr_in[ 2 ], sizeof( double ) * (size_t) npoint ); + } + +/* We will loop to apply each spectral coordinate conversion in turn to the + (spec) array. However, if the inverse transformation was requested, + we must loop through these transformations in reverse order, so set up + appropriate limits and an increment to control this loop. */ + start = forward ? 0 : map->ncvt - 1; + end = forward ? map->ncvt : -1; + inc = forward ? 1 : -1; + +/* Loop through the coordinate conversions in the required order. */ + for ( cvt = start; cvt != end; cvt += inc ) { + +/* Process conversions which correspond to changes of reference frames. */ + if( !FrameChange( map->cvttype[ cvt ], npoint, alpha, beta, spec, + map->cvtargs[ cvt ], forward, status ) ) { + +/* If this conversion was not a change of reference frame, it must be a + change of system. */ + SystemChange( map->cvttype[ cvt ], npoint, spec, + map->cvtargs[ cvt ], forward, status ); + } + } + } + +/* If an error has occurred and a new PointSet may have been created, then + clean up by annulling it. In any case, ensure that a NULL result is + returned.*/ + if ( !astOK ) { + if ( !out ) result = astAnnul( result ); + result = NULL; + } + +/* Return a pointer to the output PointSet. */ + return result; + +} + +static double UserVel( double ra, double dec, FrameDef *def, int *status ) { +/* +* Name: +* UserVel + +* Purpose: +* Find the component of the velocity of the user-defined rest-frame +* away from the source. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* double UserVel( double ra, double dec, FrameDef *def, int *status ) + +* Class Membership: +* SpecMap method. + +* Description: +* This function finds the component of the velocity of the user-defined +* rest-frame away from a specified position. The magnitude and direction +* of the rest-frames velocity are defined within the supplied "def" +* structure. The user-defined rest-frame is typically used to represent +* the velocity of the source within the heliocentric rest-frame. + +* Parameters: +* ra +* The RA (rads, FK5 J2000) of the source. +* dec +* The Dec (rads, FK5 J2000) of the source. +* def +* Pointer to a FrameDef structure which holds the parameters which +* define the frame, together with cached intermediate results. +* status +* Pointer to the inherited status variable. + +* Returns: +* The component of the frame's velocity away from the position given by +* "ra" and "dec", in m/s, measured within the Heliographic frame of +* rest. Zero is returned if an error has already occurred. + +* Notes: +* - The direction of the user velocity is given by def->refra and +* def->refdec (an FK5 J2000 position). The maginitude of the velocity +* is given by def->veluser, in m/s, positive when the source is moving +* away from the observer towards def->refra, def->refdec, and given +* with respect to the heliocentric rest-frame. + +*/ + +/* Local Variables: */ + double vb[ 3 ]; /* Source position vector */ + +/* Check the global error status. */ + if ( !astOK ) return 0.0; + +/* If not already done so, express the user velocity in the form of a + J2000.0 x,y,z vector. */ + if( def->vuser[ 0 ] == AST__BAD ) { + def->vuser[ 0 ] = def->veluser*cos( def->refra )*cos( def->refdec ); + def->vuser[ 1 ] = def->veluser*sin( def->refra )*cos( def->refdec ); + def->vuser[ 2 ] = def->veluser*sin( def->refdec ); + } + +/* Convert given J2000 RA,Dec to x,y,z. */ + palDcs2c( ra, dec, vb ); + +/* Return the dot product with the user velocity. Invert it to get the + velocity towards the observer (the def->veluser value is supposed to be + positive if the source is moving away from the observer). */ + return -palDvdv( def->vuser, vb ); +} + +/* Copy constructor. */ +/* ----------------- */ +static void Copy( const AstObject *objin, AstObject *objout, int *status ) { +/* +* Name: +* Copy + +* Purpose: +* Copy constructor for SpecMap objects. + +* Type: +* Private function. + +* Synopsis: +* void Copy( const AstObject *objin, AstObject *objout, int *status ) + +* Description: +* This function implements the copy constructor for SpecMap objects. + +* Parameters: +* objin +* Pointer to the object to be copied. +* objout +* Pointer to the object being constructed. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* void + +* Notes: +* - This constructor makes a deep copy. +*/ + +/* Local Variables: */ + AstSpecMap *in; /* Pointer to input SpecMap */ + AstSpecMap *out; /* Pointer to output SpecMap */ + int cvt; /* Loop counter for coordinate conversions */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain pointers to the input and output SpecMap structures. */ + in = (AstSpecMap *) objin; + out = (AstSpecMap *) objout; + +/* For safety, first clear any references to the input memory from the output + SpecMap. */ + out->cvtargs = NULL; + out->cvttype = NULL; + +/* Allocate memory for the output array of argument list pointers. */ + out->cvtargs = astMalloc( sizeof( double * ) * (size_t) in->ncvt ); + +/* If necessary, allocate memory and make a copy of the input array of + coordinate conversion codes. */ + if ( in->cvttype ) out->cvttype = astStore( NULL, in->cvttype, + sizeof( int ) + * (size_t) in->ncvt ); + +/* If OK, loop through each conversion in the input SpecMap and make a copy of + its argument list, storing the new pointer in the output argument list + array. */ + if ( astOK ) { + for ( cvt = 0; cvt < in->ncvt; cvt++ ) { + out->cvtargs[ cvt ] = astStore( NULL, in->cvtargs[ cvt ], + astSizeOf( in->cvtargs[ cvt ] ) ); + } + +/* If an error occurred while copying the argument lists, loop through the + conversions again and clean up by ensuring that the new memory allocated for + each argument list is freed. */ + if ( !astOK ) { + for ( cvt = 0; cvt < in->ncvt; cvt++ ) { + out->cvtargs[ cvt ] = astFree( out->cvtargs[ cvt ] ); + } + } + } + +/* If an error occurred, free all other memory allocated above. */ + if ( !astOK ) { + out->cvtargs = astFree( out->cvtargs ); + out->cvttype = astFree( out->cvttype ); + } +} + +/* Destructor. */ +/* ----------- */ +static void Delete( AstObject *obj, int *status ) { +/* +* Name: +* Delete + +* Purpose: +* Destructor for SpecMap objects. + +* Type: +* Private function. + +* Synopsis: +* void Delete( AstObject *obj, int *status ) + +* Description: +* This function implements the destructor for SpecMap objects. + +* Parameters: +* obj +* Pointer to the object to be deleted. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* void + +* Notes: +* This function attempts to execute even if the global error status is +* set. +*/ + +/* Local Variables: */ + AstSpecMap *this; /* Pointer to SpecMap */ + int cvt; /* Loop counter for coordinate conversions */ + +/* Obtain a pointer to the SpecMap structure. */ + this = (AstSpecMap *) obj; + +/* Loop to free the memory containing the argument list for each coordinate + conversion. */ + for ( cvt = 0; cvt < this->ncvt; cvt++ ) { + this->cvtargs[ cvt ] = astFree( this->cvtargs[ cvt ] ); + } + +/* Free the memory holding the array of conversion types and the array of + argument list pointers. */ + this->cvtargs = astFree( this->cvtargs ); + this->cvttype = astFree( this->cvttype ); +} + +/* Dump function. */ +/* -------------- */ +static void Dump( AstObject *this_object, AstChannel *channel, int *status ) { +/* +* Name: +* Dump + +* Purpose: +* Dump function for SpecMap objects. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* void Dump( AstObject *this, AstChannel *channel, int *status ) + +* Description: +* This function implements the Dump function which writes out data +* for the SpecMap class to an output Channel. + +* Parameters: +* this +* Pointer to the SpecMap whose data are being written. +* channel +* Pointer to the Channel to which the data are being written. +* status +* Pointer to the inherited status variable. +*/ + +/* Local Constants: */ +#define KEY_LEN 50 /* Maximum length of a keyword */ + +/* Local Variables: */ + AstSpecMap *this; /* Pointer to the SpecMap structure */ + char key[ KEY_LEN + 1 ]; /* Buffer for keyword string */ + const char *argdesc[ MAX_ARGS ]; /* Pointers to argument descriptions */ + const char *comment; /* Pointer to comment string */ + const char *sval; /* Pointer to string value */ + int argdec; /* Index of DEC argument */ + int argra; /* Index of RA argument */ + int iarg; /* Loop counter for arguments */ + int icvt; /* Loop counter for conversion steps */ + int ival; /* Integer value */ + int nargs; /* Number of user-supplied arguments */ + int szargs; /* Number of stored arguments */ + int set; /* Attribute value set? */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain a pointer to the SpecMap structure. */ + this = (AstSpecMap *) this_object; + +/* Write out values representing the instance variables for the SpecMap + class. Accompany these with appropriate comment strings, possibly + depending on the values being written.*/ + +/* In the case of attributes, we first use the appropriate (private) + Test... member function to see if they are set. If so, we then use + the (private) Get... function to obtain the value to be written + out. + + For attributes which are not set, we use the astGet... method to + obtain the value instead. This will supply a default value + (possibly provided by a derived class which over-rides this method) + which is more useful to a human reader as it corresponds to the + actual default attribute value. Since "set" will be zero, these + values are for information only and will not be read back. */ + +/* Number of conversion steps. */ +/* --------------------------- */ +/* Regard this as "set" if it is non-zero. */ + ival = this->ncvt; + set = ( ival != 0 ); + astWriteInt( channel, "Nspec", set, 0, ival, "Number of conversion steps" ); + +/* Write out data for each conversion step... */ + for ( icvt = 0; icvt < this->ncvt; icvt++ ) { + +/* Conversion type. */ +/* ---------------- */ +/* Change each conversion type code into an equivalent string and + obtain associated descriptive information. If the conversion code + was not recognised, report an error and give up. */ + if ( astOK ) { + sval = CvtString( this->cvttype[ icvt ], &comment, &argra, &argdec, + &nargs, &szargs, argdesc, status ); + if ( astOK && !sval ) { + astError( AST__SPCIN, + "astWrite(%s): Corrupt %s contains invalid SpecMap " + "spectral coordinate conversion code (%d).", status, + astGetClass( channel ), astGetClass( this ), + (int) this->cvttype[ icvt ] ); + break; + } + +/* Create an appropriate keyword and write out the conversion code + information. */ + (void) sprintf( key, "Spec%d", icvt + 1 ); + astWriteString( channel, key, 1, 1, sval, comment ); + +/* Write out data for each conversion argument... */ + for ( iarg = 0; iarg < szargs; iarg++ ) { + +/* Arguments. */ +/* ---------- */ +/* Create an appropriate keyword and write out the argument value, + accompanied by the descriptive comment obtained above. */ + if( this->cvtargs[ icvt ][ iarg ] != AST__BAD ) { + (void) sprintf( key, "Spec%d%c", icvt + 1, ALPHABET[ iarg ] ); + astWriteDouble( channel, key, 1, 1, this->cvtargs[ icvt ][ iarg ], + argdesc[ iarg ] ); + } + } + +/* Quit looping if an error occurs. */ + if ( !astOK ) break; + } + } + +/* Undefine macros local to this function. */ +#undef KEY_LEN +} + +/* Standard class functions. */ +/* ========================= */ +/* Implement the astIsASpecMap and astCheckSpecMap functions using the macros + defined for this purpose in the "object.h" header file. */ +astMAKE_ISA(SpecMap,Mapping) +astMAKE_CHECK(SpecMap) + +AstSpecMap *astSpecMap_( int nin, int flags, const char *options, int *status, ...) { +/* +*++ +* Name: +c astSpecMap +f AST_SPECMAP + +* Purpose: +* Create a SpecMap. + +* Type: +* Public function. + +* Synopsis: +c #include "specmap.h" +c AstSpecMap *astSpecMap( int nin, int flags, const char *options, ... ) +f RESULT = AST_SPECMAP( NIN, FLAGS, OPTIONS, STATUS ) + +* Class Membership: +* SpecMap constructor. + +* Description: +* This function creates a new SpecMap and optionally initialises +* its attributes. +* +* An SpecMap is a specialised form of Mapping which can be used to +* represent a sequence of conversions between standard spectral +* coordinate systems. This includes conversions between frequency, +* wavelength, and various forms of velocity, as well as conversions +* between different standards of rest. +* +* When a SpecMap is first created, it simply performs a unit +c (null) Mapping. Using the astSpecAdd function, +f (null) Mapping. Using the AST_SPECADD routine, +* a series of coordinate conversion steps may then be added, selected +* from the list of supported conversions. This allows multi-step +* conversions between a variety of spectral coordinate systems to +* be assembled out of the building blocks provided by this class. +* +* For details of the individual coordinate conversions available, +c see the description of the astSpecAdd function. +f see the description of the AST_SPECADD routine. +* +* Conversions are available to transform between standards of rest. +* Such conversions need to know the source position as an RA and DEC. +* This information can be supplied in the form of parameters for +* the relevant conversions, in which case the SpecMap is 1-dimensional, +* simply transforming the spectral axis values. This means that the +* same source position will always be used by the SpecMap. However, this +* may not be appropriate for an accurate description of a 3-D spectral +* cube, where changes of spatial position can produce significant +* changes in the Doppler shift introduced when transforming between +* standards of rest. For this situation, a 3-dimensional SpecMap can +* be created in which axes 2 and 3 correspond to the source RA and DEC +* The SpecMap simply copies values for axes 2 and 3 from input to +* output). + +* Parameters: +c nin +f NIN = INTEGER (Given) +* The number of inputs to the Mapping (this will also equal the +* number of outputs). This value must be either 1 or 3. In either +* case, the first input and output correspoindis the spectral axis. +* For a 3-axis SpecMap, the second and third axes give the RA and +* DEC (J2000 FK5) of the source. This positional information is +* used by conversions which transform between standards of rest, +* and replaces the "RA" and "DEC" arguments for the individual +* conversions listed in description of the "SpecAdd" +c function. +f routine. +c flags +f FLAGS = INTEGER (Given) +c This parameter is reserved for future use and should currently +f This argument is reserved for future use and should currently +* always be set to zero. +c options +f OPTIONS = CHARACTER * ( * ) (Given) +c Pointer to a null-terminated string containing an optional +c comma-separated list of attribute assignments to be used for +c initialising the new SpecMap. The syntax used is identical to +c that for the astSet function and may include "printf" format +c specifiers identified by "%" symbols in the normal way. +c If no initialisation is required, a zero-length string may be +c supplied. +f A character string containing an optional comma-separated +f list of attribute assignments to be used for initialising the +f new SpecMap. The syntax used is identical to that for the +f AST_SET routine. If no initialisation is required, a blank +f value may be supplied. +c ... +c If the "options" string contains "%" format specifiers, then +c an optional list of additional arguments may follow it in +c order to supply values to be substituted for these +c specifiers. The rules for supplying these are identical to +c those for the astSet function (and for the C "printf" +c function). +f STATUS = INTEGER (Given and Returned) +f The global status. + +* Returned Value: +c astSpecMap() +f AST_SPECMAP = INTEGER +* A pointer to the new SpecMap. + +* Notes: +* - The nature and units of the coordinate values supplied for the +* first input (i.e. the spectral input) of a SpecMap must be appropriate +* to the first conversion step applied by the SpecMap. For instance, if +* the first conversion step is "FRTOVL" (frequency to relativistic +* velocity), then the coordinate values for the first input should +* be frequency in units of Hz. Similarly, the nature and units of the +* coordinate values returned by a SpecMap will be determined by the +* last conversion step applied by the SpecMap. For instance, if the +* last conversion step is "VLTOVO" (relativistic velocity to optical +* velocity), then the coordinate values for the first output will be optical +* velocity in units of metres per second. See the description of the +c astSpecAdd function for the units expected and returned by each +f AST_SPECADD routine for the units expected and returned by each +* conversion. +* - A null Object pointer (AST__NULL) will be returned if this +c function is invoked with the AST error status set, or if it +f function is invoked with STATUS set to an error value, or if it +* should fail for any reason. +*-- +*/ + +/* Local Variables: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstSpecMap *new; /* Pointer to the new SpecMap */ + va_list args; /* Variable argument list */ + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(NULL); + +/* Check the global status. */ + if ( !astOK ) return NULL; + +/* Initialise the SpecMap, allocating memory and initialising the virtual + function table as well if necessary. */ + new = astInitSpecMap( NULL, sizeof( AstSpecMap ), !class_init, &class_vtab, + "SpecMap", nin, flags ); + +/* If successful, note that the virtual function table has been initialised. */ + if ( astOK ) { + class_init = 1; + +/* Obtain the variable argument list and pass it along with the options string + to the astVSet method to initialise the new SpecMap's attributes. */ + va_start( args, status ); + astVSet( new, options, NULL, args ); + va_end( args ); + +/* If an error occurred, clean up by deleting the new object. */ + if ( !astOK ) new = astDelete( new ); + } + +/* Return a pointer to the new SpecMap. */ + return new; +} + +AstSpecMap *astSpecMapId_( int nin, int flags, const char *options, ... ) { +/* +* Name: +* astSpecMapId_ + +* Purpose: +* Create a SpecMap. + +* Type: +* Private function. + +* Synopsis: +* #include "specmap.h" +* AstSpecMap *astSpecMapId_( int nin, int flags, const char *options, ... ) + +* Class Membership: +* SpecMap constructor. + +* Description: +* This function implements the external (public) interface to the +* astSpecMap constructor function. It returns an ID value (instead +* of a true C pointer) to external users, and must be provided +* because astSpecMap_ has a variable argument list which cannot be +* encapsulated in a macro (where this conversion would otherwise +* occur). +* +* The variable argument list also prevents this function from +* invoking astSpecMap_ directly, so it must be a re-implementation +* of it in all respects, except for the final conversion of the +* result to an ID value. + +* Parameters: +* As for astSpecMap_. + +* Returned Value: +* The ID value associated with the new SpecMap. +*/ + +/* Local Variables: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstSpecMap *new; /* Pointer to the new SpecMap */ + va_list args; /* Variable argument list */ + + int *status; /* Pointer to inherited status value */ + +/* Get a pointer to the inherited status value. */ + status = astGetStatusPtr; + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(NULL); + +/* Check the global status. */ + if ( !astOK ) return NULL; + +/* Initialise the SpecMap, allocating memory and initialising the virtual + function table as well if necessary. */ + new = astInitSpecMap( NULL, sizeof( AstSpecMap ), !class_init, &class_vtab, + "SpecMap", nin, flags ); + +/* If successful, note that the virtual function table has been initialised. */ + if ( astOK ) { + class_init = 1; + +/* Obtain the variable argument list and pass it along with the options string + to the astVSet method to initialise the new SpecMap's attributes. */ + va_start( args, options ); + astVSet( new, options, NULL, args ); + va_end( args ); + +/* If an error occurred, clean up by deleting the new object. */ + if ( !astOK ) new = astDelete( new ); + } + +/* Return an ID value for the new SpecMap. */ + return astMakeId( new ); +} + +AstSpecMap *astInitSpecMap_( void *mem, size_t size, int init, + AstSpecMapVtab *vtab, const char *name, + int nin, int flags, int *status ) { +/* +*+ +* Name: +* astInitSpecMap + +* Purpose: +* Initialise a SpecMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "specmap.h" +* AstSpecMap *astInitSpecMap( void *mem, size_t size, int init, +* AstSpecMapVtab *vtab, const char *name, +* int nin, int flags ) + +* Class Membership: +* SpecMap initialiser. + +* Description: +* This function is provided for use by class implementations to initialise +* a new SpecMap object. It allocates memory (if necessary) to accommodate +* the SpecMap plus any additional data associated with the derived class. +* It then initialises a SpecMap structure at the start of this memory. If +* the "init" flag is set, it also initialises the contents of a virtual +* function table for a SpecMap at the start of the memory passed via the +* "vtab" parameter. + +* Parameters: +* mem +* A pointer to the memory in which the SpecMap is to be initialised. +* This must be of sufficient size to accommodate the SpecMap data +* (sizeof(SpecMap)) plus any data used by the derived class. If a value +* of NULL is given, this function will allocate the memory itself using +* the "size" parameter to determine its size. +* size +* The amount of memory used by the SpecMap (plus derived class data). +* This will be used to allocate memory if a value of NULL is given for +* the "mem" parameter. This value is also stored in the SpecMap +* structure, so a valid value must be supplied even if not required for +* allocating memory. +* init +* A logical flag indicating if the SpecMap's virtual function table is +* to be initialised. If this value is non-zero, the virtual function +* table will be initialised by this function. +* vtab +* Pointer to the start of the virtual function table to be associated +* with the new SpecMap. +* name +* Pointer to a constant null-terminated character string which contains +* the name of the class to which the new object belongs (it is this +* pointer value that will subsequently be returned by the astClass +* method). +* nin +* The number of inputs and outputs for the SpecMap (either 1 or 3). +* flags +* This parameter is reserved for future use. It is currently ignored. + +* Returned Value: +* A pointer to the new SpecMap. + +* Notes: +* - A null pointer will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +*- +*/ + +/* Local Variables: */ + AstSpecMap *new; /* Pointer to the new SpecMap */ + +/* Check the global status. */ + if ( !astOK ) return NULL; + +/* Check nin is OK (1 or 3). */ + if( nin != 1 && nin != 3 ) { + astError( AST__BADNI, "astInitSpecMap(SpecMap): Supplied number of " + "SpecMap axes (%d) is illegal; it should be 1 or 2. ", status, + nin ); + } + +/* If necessary, initialise the virtual function table. */ + if ( init ) astInitSpecMapVtab( vtab, name ); + +/* Initialise a 1D Mapping structure (the parent class) as the first component + within the SpecMap structure, allocating memory if necessary. Specify that + the Mapping should be defined in both the forward and inverse directions. */ + new = (AstSpecMap *) astInitMapping( mem, size, 0, + (AstMappingVtab *) vtab, name, + nin, nin, 1, 1 ); + + if ( astOK ) { + +/* Initialise the SpecMap data. */ +/* --------------------------- */ +/* The initial state is with no conversions set, in which condition the + SpecMap simply implements a unit mapping. */ + new->ncvt = 0; + new->cvtargs = NULL; + new->cvttype = NULL; + +/* If an error occurred, clean up by deleting the new object. */ + if ( !astOK ) new = astDelete( new ); + } + +/* Return a pointer to the new object. */ + return new; +} + +AstSpecMap *astLoadSpecMap_( void *mem, size_t size, + AstSpecMapVtab *vtab, const char *name, + AstChannel *channel, int *status ) { +/* +*+ +* Name: +* astLoadSpecMap + +* Purpose: +* Load a SpecMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "specmap.h" +* AstSpecMap *astLoadSpecMap( void *mem, size_t size, +* AstSpecMapVtab *vtab, const char *name, +* AstChannel *channel ) + +* Class Membership: +* SpecMap loader. + +* Description: +* This function is provided to load a new SpecMap using data read +* from a Channel. It first loads the data used by the parent class +* (which allocates memory if necessary) and then initialises a +* SpecMap structure in this memory, using data read from the input +* Channel. +* +* If the "init" flag is set, it also initialises the contents of a +* virtual function table for a SpecMap at the start of the memory +* passed via the "vtab" parameter. + + +* Parameters: +* mem +* A pointer to the memory into which the SpecMap is to be +* loaded. This must be of sufficient size to accommodate the +* SpecMap data (sizeof(SpecMap)) plus any data used by derived +* classes. If a value of NULL is given, this function will +* allocate the memory itself using the "size" parameter to +* determine its size. +* size +* The amount of memory used by the SpecMap (plus derived class +* data). This will be used to allocate memory if a value of +* NULL is given for the "mem" parameter. This value is also +* stored in the SpecMap structure, so a valid value must be +* supplied even if not required for allocating memory. +* +* If the "vtab" parameter is NULL, the "size" value is ignored +* and sizeof(AstSpecMap) is used instead. +* vtab +* Pointer to the start of the virtual function table to be +* associated with the new SpecMap. If this is NULL, a pointer to +* the (static) virtual function table for the SpecMap class is +* used instead. +* name +* Pointer to a constant null-terminated character string which +* contains the name of the class to which the new object +* belongs (it is this pointer value that will subsequently be +* returned by the astGetClass method). +* +* If the "vtab" parameter is NULL, the "name" value is ignored +* and a pointer to the string "SpecMap" is used instead. + +* Returned Value: +* A pointer to the new SpecMap. + +* Notes: +* - A null pointer will be returned if this function is invoked +* with the global error status set, or if it should fail for any +* reason. +*- +*/ + +/* Local Constants: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ +#define KEY_LEN 50 /* Maximum length of a keyword */ + +/* Local Variables: */ + AstSpecMap *new; /* Pointer to the new SpecMap */ + char *sval; /* Pointer to string value */ + char key[ KEY_LEN + 1 ]; /* Buffer for keyword string */ + const char *argdesc[ MAX_ARGS ]; /* Pointers to argument descriptions */ + const char *comment; /* Pointer to comment string */ + int argdec; /* Index of DEC argument */ + int argra; /* Index of RA argument */ + int iarg; /* Loop counter for arguments */ + int icvt; /* Loop counter for conversion steps */ + int nargs; /* Number of user-supplied arguments */ + int szargs; /* Number of stored arguments */ + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(channel); + +/* Initialise. */ + new = NULL; + +/* Check the global error status. */ + if ( !astOK ) return new; + +/* If a NULL virtual function table has been supplied, then this is + the first loader to be invoked for this SpecMap. In this case the + SpecMap belongs to this class, so supply appropriate values to be + passed to the parent class loader (and its parent, etc.). */ + if ( !vtab ) { + size = sizeof( AstSpecMap ); + vtab = &class_vtab; + name = "SpecMap"; + +/* If required, initialise the virtual function table for this class. */ + if ( !class_init ) { + astInitSpecMapVtab( vtab, name ); + class_init = 1; + } + } + +/* Invoke the parent class loader to load data for all the ancestral + classes of the current one, returning a pointer to the resulting + partly-built SpecMap. */ + new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name, + channel ); + + if ( astOK ) { + +/* Read input data. */ +/* ================ */ +/* Request the input Channel to read all the input data appropriate to + this class into the internal "values list". */ + astReadClassData( channel, "SpecMap" ); + +/* Now read each individual data item from this list and use it to + initialise the appropriate instance variable(s) for this class. */ + +/* In the case of attributes, we first read the "raw" input value, + supplying the "unset" value as the default. If a "set" value is + obtained, we then use the appropriate (private) Set... member + function to validate and set the value properly. */ + +/* Number of conversion steps. */ +/* --------------------------- */ +/* Read the number of conversion steps and allocate memory to hold + data for each step. */ + new->ncvt = astReadInt( channel, "nspec", 0 ); + if ( new->ncvt < 0 ) new->ncvt = 0; + new->cvttype = astMalloc( sizeof( int ) * (size_t) new->ncvt ); + new->cvtargs = astMalloc( sizeof( double * ) * (size_t) new->ncvt ); + +/* If an error occurred, ensure that all allocated memory is freed. */ + if ( !astOK ) { + new->cvttype = astFree( new->cvttype ); + new->cvtargs = astFree( new->cvtargs ); + +/* Otherwise, initialise the argument pointer array. */ + } else { + for ( icvt = 0; icvt < new->ncvt; icvt++ ) { + new->cvtargs[ icvt ] = NULL; + } + +/* Read in data for each conversion step... */ + for ( icvt = 0; icvt < new->ncvt; icvt++ ) { + +/* Conversion type. */ +/* ---------------- */ +/* Create an appropriate keyword and read the string representation of + the conversion type. */ + (void) sprintf( key, "spec%d", icvt + 1 ); + sval = astReadString( channel, key, NULL ); + +/* If no value was read, report an error. */ + if ( astOK ) { + if ( !sval ) { + astError( AST__BADIN, + "astRead(%s): A spectral coordinate conversion " + "type is missing from the input SpecMap data.", status, + astGetClass( channel ) ); + +/* Otherwise, convert the string representation into the required + conversion type code. */ + } else { + new->cvttype[ icvt ] = CvtCode( sval, status ); + +/* If the string was not recognised, report an error. */ + if ( new->cvttype[ icvt ] == AST__SPEC_NULL ) { + astError( AST__BADIN, + "astRead(%s): Invalid spectral conversion " + "type \"%s\" in SpecMap data.", status, + astGetClass( channel ), sval ); + } + } + +/* Free the memory holding the string value. */ + sval = astFree( sval ); + } + +/* Obtain the number of arguments associated with the conversion and + allocate memory to hold them. */ + (void) CvtString( new->cvttype[ icvt ], &comment, &argra, + &argdec, &nargs, &szargs, argdesc, status ); + new->cvtargs[ icvt ] = astMalloc( sizeof( double ) * + (size_t) szargs ); + +/* Read in data for each argument... */ + if ( astOK ) { + for ( iarg = 0; iarg < szargs; iarg++ ) { + +/* Arguments. */ +/* ---------- */ +/* Create an appropriate keyword and read each argument value. */ + (void) sprintf( key, "spec%d%c", icvt + 1, ALPHABET[ iarg ] ); + new->cvtargs[ icvt ][ iarg ] = astReadDouble( channel, key, + AST__BAD ); + } + } + +/* Quit looping if an error occurs. */ + if ( !astOK ) break; + } + } + +/* If an error occurred, clean up by deleting the new SpecMap. */ + if ( !astOK ) new = astDelete( new ); + } + +/* Return the new SpecMap pointer. */ + return new; + +/* Undefine macros local to this function. */ +#undef KEY_LEN +} + +/* Virtual function interfaces. */ +/* ============================ */ +/* These provide the external interface to the virtual functions defined by + this class. Each simply checks the global error status and then locates and + executes the appropriate member function, using the function pointer stored + in the object's virtual function table (this pointer is located using the + astMEMBER macro defined in "object.h"). + + Note that the member function may not be the one defined here, as it may + have been over-ridden by a derived class. However, it should still have the + same interface. */ +void astSpecAdd_( AstSpecMap *this, const char *cvt, const double args[], int *status ) { + if ( !astOK ) return; + (**astMEMBER(this,SpecMap,SpecAdd))( this, cvt, args, status ); +} + + + + |