diff options
author | William Joye <wjoye@cfa.harvard.edu> | 2017-12-08 18:57:06 (GMT) |
---|---|---|
committer | William Joye <wjoye@cfa.harvard.edu> | 2017-12-08 18:57:06 (GMT) |
commit | 90a861b642f765d5657ab827aedabe3920ff9333 (patch) | |
tree | 88b93d468ca1feed91ef2958f46f3f74f1418aac /ast/specmap.c | |
parent | fba23129f50db253ed3fbbaa23d6e342bf86068e (diff) | |
download | blt-90a861b642f765d5657ab827aedabe3920ff9333.zip blt-90a861b642f765d5657ab827aedabe3920ff9333.tar.gz blt-90a861b642f765d5657ab827aedabe3920ff9333.tar.bz2 |
upgrade AST
Diffstat (limited to 'ast/specmap.c')
-rw-r--r-- | ast/specmap.c | 4671 |
1 files changed, 0 insertions, 4671 deletions
diff --git a/ast/specmap.c b/ast/specmap.c deleted file mode 100644 index 164179a..0000000 --- a/ast/specmap.c +++ /dev/null @@ -1,4671 +0,0 @@ -/* -*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 ); -} - - - - |