1 /*
2 *class++
3 * Name:
4 * SpecMap
5
6 * Purpose:
7 * Sequence of spectral coordinate conversions.
8
9 * Constructor Function:
10 c astSpecMap (also see astSpecAdd)
11 f AST_SPECMAP (also see AST_SPECADD)
12
13 * Description:
14 * A SpecMap is a specialised form of Mapping which can be used to
15 * represent a sequence of conversions between standard spectral
16 * coordinate systems.
17 *
18 * When an SpecMap is first created, it simply performs a unit
19 c (null) Mapping. Using the astSpecAdd
20 f (null) Mapping. Using the AST_SPECADD
21 c function, a series of coordinate conversion steps may then be
22 f routine, a series of coordinate conversion steps may then be
23 * added. This allows multi-step conversions between a variety of
24 * spectral coordinate systems to be assembled out of a set of building
25 * blocks.
26 *
27 * Conversions are available to transform between standards of rest.
28 * Such conversions need to know the source position as an RA and DEC.
29 * This information can be supplied in the form of parameters for
30 * the relevant conversions, in which case the SpecMap is 1-dimensional,
31 * simply transforming the spectral axis values. This means that the
32 * same source position will always be used by the SpecMap. However, this
33 * may not be appropriate for an accurate description of a 3-D spectral
34 * cube, where changes of spatial position can produce significant
35 * changes in the Doppler shift introduced when transforming between
36 * standards of rest. For this situation, a 3-dimensional SpecMap can
37 * be created in which axes 2 and 3 correspond to the source RA and DEC
38 * The SpecMap simply copies values for axes 2 and 3 from input to
39 * output), but modifies axis 1 values (the spectral axis) appropriately.
40 *
41 * For details of the individual coordinate conversions available,
42 c see the description of the astSpecAdd function.
43 f see the description of the AST_SPECADD routine.
44
45 * Inheritance:
46 * The SpecMap class inherits from the Mapping class.
47
48 * Attributes:
49 * The SpecMap class does not define any new attributes beyond those
50 * which are applicable to all Mappings.
51
52 * Functions:
53 c In addition to those functions applicable to all Mappings, the
54 c following function may also be applied to all SpecMaps:
55 f In addition to those routines applicable to all Mappings, the
56 f following routine may also be applied to all SpecMaps:
57 *
58 c - astSpecAdd: Add a spectral coordinate conversion to an SpecMap
59 f - AST_SPECADD: Add a spectral coordinate conversion to an SpecMap
60
61 * Copyright:
62 * Copyright (C) 1997-2006 Council for the Central Laboratory of the
63 * Research Councils
64 * Copyright (C) 2009 Science & Technology Facilities Council.
65 * All Rights Reserved.
66
67 * Licence:
68 * This program is free software: you can redistribute it and/or
69 * modify it under the terms of the GNU Lesser General Public
70 * License as published by the Free Software Foundation, either
71 * version 3 of the License, or (at your option) any later
72 * version.
73 *
74 * This program is distributed in the hope that it will be useful,
75 * but WITHOUT ANY WARRANTY; without even the implied warranty of
76 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
77 * GNU Lesser General Public License for more details.
78 *
79 * You should have received a copy of the GNU Lesser General
80 * License along with this program. If not, see
81 * <http://www.gnu.org/licenses/>.
82
83 * Authors:
84 * DSB: David S. Berry (Starlink)
85
86 * History:
87 * 6-NOV-2002 (DSB):
88 * Original version.
89 * 14-JUL-2003 (DSB):
90 * Added checks for NAN values produced by transformation functions.
91 * 17-SEP-2003 (DSB):
92 * - Improve FRTOAW accuracy by iterating.
93 * - Changed Refrac to use algorithm given in FITS-WCS paper 3
94 * version dated 21/9/03.
95 * 14-FEB-2006 (DSB):
96 * Override astGetObjSize.
97 * 10-MAY-2006 (DSB):
98 * Override astEqual.
99 * 15-NOV-2006 (DSB):
100 * Guard against division by zero when converting freq to wave in
101 * SystemChange.
102 * 18-JUN-2009 (DSB):
103 * Add OBSALT argument to TPF2HL and HLF2TP conversions.
104 * Change GEOLON/LAT to OBSLON/LAT for consistency with other
105 * classes.
106 * 2-OCT-2012 (DSB):
107 * Check for Infs as well as NaNs.
108 *class--
109 */
110
111 /* Module Macros. */
112 /* ============== */
113 /* Set the name of the class we are implementing. This indicates to
114 the header files that define class interfaces that they should make
115 "protected" symbols available. */
116 #define astCLASS SpecMap
117
118 /* Codes to identify spectral coordinate conversions. */
119 #define AST__SPEC_NULL 0 /* Null value */
120 #define AST__FRTOVL 1 /* Frequency to relativistic velocity */
121 #define AST__VLTOFR 2 /* Relativistic velocity to Frequency */
122 #define AST__ENTOFR 3 /* Energy to frequency */
123 #define AST__FRTOEN 4 /* Frequency to energy */
124 #define AST__WNTOFR 5 /* Wave number to frequency */
125 #define AST__FRTOWN 6 /* Frequency to wave number */
126 #define AST__WVTOFR 7 /* Wavelength (vacuum) to frequency */
127 #define AST__FRTOWV 8 /* Frequency to wavelength (vacuum) */
128 #define AST__AWTOFR 9 /* Wavelength (air) to frequency */
129 #define AST__FRTOAW 10 /* Frequency to wavelength (air) */
130 #define AST__VRTOVL 11 /* Radio to relativistic velocity */
131 #define AST__VLTOVR 12 /* Relativistic to radio velocity */
132 #define AST__VOTOVL 13 /* Optical to relativistic velocity */
133 #define AST__VLTOVO 14 /* Relativistic to optical velocity */
134 #define AST__ZOTOVL 15 /* Redshift to relativistic velocity */
135 #define AST__VLTOZO 16 /* Relativistic velocity to redshift */
136 #define AST__BTTOVL 17 /* Beta factor to relativistic velocity */
137 #define AST__VLTOBT 18 /* Relativistic velocity to beta factor */
138 #define AST__USF2HL 19 /* User-defined to heliocentric frequency */
139 #define AST__HLF2US 20 /* Heliocentric to user-defined frequency */
140 #define AST__TPF2HL 21 /* Topocentric to heliocentric frequency */
141 #define AST__HLF2TP 22 /* Heliocentric to topocentric frequency */
142 #define AST__GEF2HL 23 /* Geocentric to heliocentric frequency */
143 #define AST__HLF2GE 24 /* Heliocentric to geocentric frequency */
144 #define AST__BYF2HL 25 /* Barycentric to heliocentric frequency */
145 #define AST__HLF2BY 26 /* Heliocentric to barycentric frequency */
146 #define AST__LKF2HL 27 /* LSRK to heliocentric frequency */
147 #define AST__HLF2LK 28 /* Heliocentric to LSRK frequency */
148 #define AST__LDF2HL 29 /* LSRD to heliocentric frequency */
149 #define AST__HLF2LD 30 /* Heliocentric to LSRD frequency */
150 #define AST__LGF2HL 31 /* Local group to heliocentric frequency */
151 #define AST__HLF2LG 32 /* Heliocentric to local group frequency */
152 #define AST__GLF2HL 33 /* Galactic to heliocentric frequency */
153 #define AST__HLF2GL 34 /* Heliocentric to galactic frequency */
154
155 /* Maximum number of arguments required by a conversion. */
156 #define MAX_ARGS 7
157
158 /* The alphabet (used for generating keywords for arguments). */
159 #define ALPHABET "abcdefghijklmnopqrstuvwxyz"
160
161 /* Angle conversion */
162 #define PI 3.141592653589793238462643
163 #define PIBY2 (PI/2.0)
164 #define D2R (PI/180.0)
165 #define R2D (180.0/PI)
166
167 /* Macros which return the maximum and minimum of two values. */
168 #define MAX(aa,bb) ((aa)>(bb)?(aa):(bb))
169 #define MIN(aa,bb) ((aa)<(bb)?(aa):(bb))
170
171 /* Macro to check for equality of floating point values. We cannot
172 compare bad values directory because of the danger of floating point
173 exceptions, so bad values are dealt with explicitly. */
174 #define EQUAL(aa,bb) (((aa)==AST__BAD)?(((bb)==AST__BAD)?1:0):(((bb)==AST__BAD)?0:(fabs((aa)-(bb))<=1.0E5*MAX((fabs(aa)+fabs(bb))*DBL_EPSILON,DBL_MIN))))
175
176 /* Include files. */
177 /* ============== */
178 /* Interface definitions. */
179 /* ---------------------- */
180 #include "pal.h" /* SLALIB interface */
181
182 #include "globals.h" /* Thread-safe global data access */
183 #include "error.h" /* Error reporting facilities */
184 #include "memory.h" /* Memory allocation facilities */
185 #include "object.h" /* Base Object class */
186 #include "pointset.h" /* Sets of points/coordinates */
187 #include "mapping.h" /* Coordinate Mappings (parent class) */
188 #include "unitmap.h" /* Unit (null) Mappings */
189 #include "specmap.h" /* Interface definition for this class */
190
191 /* Error code definitions. */
192 /* ----------------------- */
193 #include "ast_err.h" /* AST error codes */
194
195 /* C header files. */
196 /* --------------- */
197 #include <ctype.h>
198 #include <stddef.h>
199 #include <stdio.h>
200 #include <string.h>
201 #include <math.h>
202
203 /* Module Variables. */
204 /* ================= */
205
206 /* Address of this static variable is used as a unique identifier for
207 member of this class. */
208 static int class_check;
209
210 /* Pointers to parent class methods which are extended by this class. */
211 static int (* parent_getobjsize)( AstObject *, int * );
212 static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
213 static double (* parent_rate)( AstMapping *, double *, int, int, int * );
214
215
216 #ifdef THREAD_SAFE
217 /* Define how to initialise thread-specific globals. */
218 #define GLOBAL_inits \
219 globals->Class_Init = 0;
220
221 /* Create the function that initialises global data for this module. */
222 astMAKE_INITGLOBALS(SpecMap)
223
224 /* Define macros for accessing each item of thread specific global data. */
225 #define class_init astGLOBAL(SpecMap,Class_Init)
226 #define class_vtab astGLOBAL(SpecMap,Class_Vtab)
227
228
229 #include <pthread.h>
230
231
232 #else
233
234
235 /* Define the class virtual function table and its initialisation flag
236 as static variables. */
237 static AstSpecMapVtab class_vtab; /* Virtual function table */
238 static int class_init = 0; /* Virtual function table initialised? */
239
240 #endif
241
242 /* Structure to hold parameters and intermediate values describing a
243 reference frame */
244 typedef struct FrameDef {
245 double obsalt; /* Observers geodetic altitude (m) */
246 double obslat; /* Observers geodetic latitude (rads) */
247 double obslon; /* Observers geodetic longitude (rads, +ve east) */
248 double epoch; /* Julian epoch of observation */
249 double refdec; /* RA of reference point (FK5 J2000) */
250 double refra; /* DEC of reference point (FK5 J2000) */
251 double veluser; /* Heliocentric velocity of user-defined system (m/s) */
252 double last; /* Local apparent sideral time */
253 double amprms[21]; /* Mean to apparent parameters */
254 double vuser[3]; /* Used-defined velocity as a FK5 J2000 vector */
255 double dvh[3]; /* Earth-sun velocity */
256 double dvb[3]; /* Barycentre-sun velocity */
257 } FrameDef;
258
259 /* External Interface Function Prototypes. */
260 /* ======================================= */
261 /* The following functions have public prototypes only (i.e. no
262 protected prototypes), so we must provide local prototypes for use
263 within this module. */
264 AstSpecMap *astSpecMapId_( int, int, const char *, ... );
265
266 /* Prototypes for Private Member Functions. */
267 /* ======================================== */
268 static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
269 static const char *CvtString( int, const char **, int *, int *, int *, int *, const char *[ MAX_ARGS ], int * );
270 static double BaryVel( double, double, FrameDef *, int * );
271 static double GalVel( double, double, FrameDef *, int * );
272 static double GeoVel( double, double, FrameDef *, int * );
273 static double LgVel( double, double, FrameDef *, int * );
274 static double LsrdVel( double, double, FrameDef *, int * );
275 static double LsrkVel( double, double, FrameDef *, int * );
276 static double Rate( AstMapping *, double *, int, int, int * );
277 static double Refrac( double, int * );
278 static double Rverot( double, double, double, double, double, int * );
279 static double TopoVel( double, double, FrameDef *, int * );
280 static double UserVel( double, double, FrameDef *, int * );
281 static int CvtCode( const char *, int * );
282 static int Equal( AstObject *, AstObject *, int * );
283 static int FrameChange( int, int, double *, double *, double *, double *, int, int * );
284 static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * );
285 static int SystemChange( int, int, double *, double *, int, int * );
286 static void AddSpecCvt( AstSpecMap *, int, const double *, int * );
287 static void Copy( const AstObject *, AstObject *, int * );
288 static void Delete( AstObject *, int * );
289 static void Dump( AstObject *, AstChannel *, int * );
290 static void SpecAdd( AstSpecMap *, const char *, const double[], int * );
291
292 static int GetObjSize( AstObject *, int * );
293 /* Member functions. */
294 /* ================= */
Equal(AstObject * this_object,AstObject * that_object,int * status)295 static int Equal( AstObject *this_object, AstObject *that_object, int *status ) {
296 /*
297 * Name:
298 * Equal
299
300 * Purpose:
301 * Test if two SpecMaps are equivalent.
302
303 * Type:
304 * Private function.
305
306 * Synopsis:
307 * #include "specmap.h"
308 * int Equal( AstObject *this, AstObject *that, int *status )
309
310 * Class Membership:
311 * SpecMap member function (over-rides the astEqual protected
312 * method inherited from the astMapping class).
313
314 * Description:
315 * This function returns a boolean result (0 or 1) to indicate whether
316 * two SpecMaps are equivalent.
317
318 * Parameters:
319 * this
320 * Pointer to the first Object (a SpecMap).
321 * that
322 * Pointer to the second Object.
323 * status
324 * Pointer to the inherited status variable.
325
326 * Returned Value:
327 * One if the SpecMaps are equivalent, zero otherwise.
328
329 * Notes:
330 * - A value of zero will be returned if this function is invoked
331 * with the global status set, or if it should fail for any reason.
332 */
333
334 /* Local Variables: */
335 AstSpecMap *that;
336 AstSpecMap *this;
337 const char *argdesc[ MAX_ARGS ];
338 const char *comment;
339 int argdec;
340 int argra;
341 int i, j;
342 int nargs;
343 int nin;
344 int nout;
345 int result;
346 int szargs;
347
348 /* Initialise. */
349 result = 0;
350
351 /* Check the global error status. */
352 if ( !astOK ) return result;
353
354 /* Obtain pointers to the two SpecMap structures. */
355 this = (AstSpecMap *) this_object;
356 that = (AstSpecMap *) that_object;
357
358 /* Check the second object is a SpecMap. We know the first is a
359 SpecMap since we have arrived at this implementation of the virtual
360 function. */
361 if( astIsASpecMap( that ) ) {
362
363 /* Get the number of inputs and outputs and check they are the same for both. */
364 nin = astGetNin( this );
365 nout = astGetNout( this );
366 if( astGetNin( that ) == nin && astGetNout( that ) == nout ) {
367
368 /* If the Invert flags for the two SpecMaps differ, it may still be possible
369 for them to be equivalent. First compare the SpecMaps if their Invert
370 flags are the same. In this case all the attributes of the two SpecMaps
371 must be identical. */
372 if( astGetInvert( this ) == astGetInvert( that ) ) {
373 if( this->ncvt == that->ncvt ) {
374 result = 1;
375 for( i = 0; i < this->ncvt && result; i++ ) {
376 if( this->cvttype[ i ] != that->cvttype[ i ] ) {
377 result = 0;
378 } else {
379 CvtString( this->cvttype[ i ], &comment, &argra,
380 &argdec, &nargs, &szargs, argdesc, status );
381 for( j = 0; j < nargs; j++ ) {
382 if( !astEQUAL( this->cvtargs[ i ][ j ],
383 that->cvtargs[ i ][ j ] ) ){
384 result = 0;
385 break;
386 }
387 }
388 }
389 }
390 }
391
392 /* If the Invert flags for the two SpecMaps differ, the attributes of the two
393 SpecMaps must be inversely related to each other. */
394 } else {
395
396 /* In the specific case of a SpecMap, Invert flags must be equal. */
397 result = 0;
398
399 }
400 }
401 }
402
403 /* If an error occurred, clear the result value. */
404 if ( !astOK ) result = 0;
405
406 /* Return the result, */
407 return result;
408 }
409
GetObjSize(AstObject * this_object,int * status)410 static int GetObjSize( AstObject *this_object, int *status ) {
411 /*
412 * Name:
413 * GetObjSize
414
415 * Purpose:
416 * Return the in-memory size of an Object.
417
418 * Type:
419 * Private function.
420
421 * Synopsis:
422 * #include "specmap.h"
423 * int GetObjSize( AstObject *this, int *status )
424
425 * Class Membership:
426 * SpecMap member function (over-rides the astGetObjSize protected
427 * method inherited from the parent class).
428
429 * Description:
430 * This function returns the in-memory size of the supplied SpecMap,
431 * in bytes.
432
433 * Parameters:
434 * this
435 * Pointer to the SpecMap.
436 * status
437 * Pointer to the inherited status variable.
438
439 * Returned Value:
440 * The Object size, in bytes.
441
442 * Notes:
443 * - A value of zero will be returned if this function is invoked
444 * with the global status set, or if it should fail for any reason.
445 */
446
447 /* Local Variables: */
448 AstSpecMap *this; /* Pointer to SpecMap structure */
449 int result; /* Result value to return */
450 int cvt; /* Loop counter for coordinate conversions */
451
452 /* Initialise. */
453 result = 0;
454
455 /* Check the global error status. */
456 if ( !astOK ) return result;
457
458 /* Obtain a pointers to the SpecMap structure. */
459 this = (AstSpecMap *) this_object;
460
461 /* Invoke the GetObjSize method inherited from the parent class, and then
462 add on any components of the class structure defined by thsi class
463 which are stored in dynamically allocated memory. */
464 result = (*parent_getobjsize)( this_object, status );
465
466 for ( cvt = 0; cvt < this->ncvt; cvt++ ) {
467 result += astTSizeOf( this->cvtargs[ cvt ] );
468 }
469
470 result += astTSizeOf( this->cvtargs );
471 result += astTSizeOf( this->cvttype );
472
473 /* If an error occurred, clear the result value. */
474 if ( !astOK ) result = 0;
475
476 /* Return the result, */
477 return result;
478 }
479
AddSpecCvt(AstSpecMap * this,int cvttype,const double * args,int * status)480 static void AddSpecCvt( AstSpecMap *this, int cvttype, const double *args, int *status ) {
481 /*
482 * Name:
483 * AddSpecCvt
484
485 * Purpose:
486 * Add a coordinate conversion step to an SpecMap.
487
488 * Type:
489 * Private function.
490
491 * Synopsis:
492 * #include "specmap.h"
493 * void AddSpecCvt( AstSpecMap *this, int cvttype, const double *args )
494
495 * Class Membership:
496 * SpecMap member function.
497
498 * Description:
499 * This function allows one of the supported spectral coordinate
500 * conversions to be appended to a SpecMap. When a SpecMap is first
501 * created (using astSpecMap), it simply performs a unit mapping. By
502 * using AddSpecCvt repeatedly, a series of coordinate conversions may
503 * then be specified which the SpecMap will subsequently perform in
504 * sequence. This allows a complex coordinate conversion to be
505 * assembled out of the basic building blocks. The SpecMap will also
506 * perform the inverse coordinate conversion (applying the individual
507 * conversion steps in reverse) if required.
508
509 * Parameters:
510 * this
511 * Pointer to the SpecMap.
512 * cvttype
513 * A code to identify which spectral coordinate conversion is to be
514 * appended. See the "Coordinate Conversions" section for details
515 * of those available.
516 * args
517 * Pointer to an array of double containing the argument values
518 * required to fully specify the required coordinate
519 * conversion. The number of arguments depends on the conversion
520 * (see the "Coordinate Conversions" section for details). This
521 * value is ignored and may be NULL if no arguments are required.
522
523 * Returned Value:
524 * void.
525
526 * Coordinate Conversions:
527 * The following values may be supplied for the "cvttype" parameter
528 * in order to specify the coordinate conversion to be performed.
529 * The argument(s) required to fully specify each conversion are
530 * indicated in parentheses after each value, and described at the end
531 * of the list. Values for these should be given in the array pointed
532 * at by "args".
533 *
534 * AST__FRTOVL( RF )
535 * Convert frequency to relativistic velocity.
536 * AST__VLTOFR( RF )
537 * Convert relativistic velocity to Frequency.
538 * AST__ENTOFR
539 * Convert energy to frequency.
540 * AST__FRTOEN
541 * Convert frequency to energy.
542 * AST__WNTOFR
543 * Convert wave number to frequency.
544 * AST__FRTOWN
545 * Convert frequency to wave number.
546 * AST__WVTOFR
547 * Convert wavelength (vacuum) to frequency.
548 * AST__FRTOWV
549 * Convert frequency to wavelength (vacuum).
550 * AST__AWTOFR
551 * Convert wavelength (air) to frequency.
552 * AST__FRTOAW
553 * Convert frequency to wavelength (air).
554 * AST__VRTOVL
555 * Convert radio to relativistic velocity.
556 * AST__VLTOVR
557 * Convert relativistic to radio velocity.
558 * AST__VOTOVL
559 * Convert optical to relativistic velocity.
560 * AST__VLTOVO
561 * Convert relativistic to optical velocity.
562 * AST__ZOTOVL
563 * Convert redshift to relativistic velocity.
564 * AST__VLTOZO
565 * Convert relativistic velocity to redshift.
566 * AST__BTTOVL
567 * Convert beta factor to relativistic velocity.
568 * AST__VLTOBT
569 * Convert relativistic velocity to beta factor.
570 * AST_USF2HL( VOFF, RA, DEC )
571 * Convert frequency from a user-defined reference frame to
572 * heliocentric.
573 * AST__HLF2US( VOFF, RA, DEC )
574 * Convert frequency from heliocentric reference frame to
575 * user-defined.
576 * AST__TPF2HL( OBSLON, OBSLAT, OBSALT, EPOCH, RA, DEC )
577 * Convert from Topocentric to heliocentric frequency
578 * AST__HLF2TP( OBSLON, OBSLAT, OBSALT, EPOCH, RA, DEC )
579 * Convert from Heliocentric to topocentric frequency.
580 * AST__GEF2HL( EPOCH, RA, DEC )
581 * Convert from Geocentric to heliocentric frequency.
582 * AST__HLF2GE( EPOCH, RA, DEC )
583 * Convert from Heliocentric to geocentric frequency.
584 * AST__BYF2HL( EPOCH, RA, DEC )
585 * Convert from Barycentric to heliocentric frequency.
586 * AST__HLF2BY( EPOCH, RA, DEC )
587 * Convert from Heliocentric to barycentric frequency.
588 * AST__LKF2HL( RA, DEC )
589 * Convert from LSRK to heliocentric frequency.
590 * AST__HLF2LK( RA, DEC )
591 * Convert from Heliocentric to LSRK frequency.
592 * AST__LDF2HL( RA, DEC )
593 * Convert from LSRD to heliocentric frequency.
594 * AST__HLF2LD( RA, DEC )
595 * Convert from Heliocentric to LSRD frequency.
596 * AST__LGF2HL( RA, DEC )
597 * Convert from Local group to heliocentric frequency.
598 * AST__HLF2LG( RA, DEC )
599 * Convert from Heliocentric to local group frequency.
600 * AST__GLF2HL( RA, DEC )
601 * Convert from Galactic to heliocentric frequency.
602 * AST__HLF2GL( RA, DEC )
603 * Convert from Heliocentric to galactic frequency.
604 *
605 * The units for the values processed by the above conversions are as
606 * follows:
607 *
608 * - all velocities: metres per second.
609 * - frequency: Hertz.
610 * - all wavelengths: metres.
611 * - energy: Joules.
612 * - wave number: cycles per metre.
613 *
614 * The arguments used in the above conversions are as follows:
615 *
616 * - RF: Rest frequency (Hz).
617 * - OBSALT: Geodetic altitude of observer (IAU 1975, metres).
618 * - OBSLAT: Geodetic latitude of observer (IAU 1975, radians).
619 * - OBSLON: Longitude of observer (radians, positive eastwards).
620 * - EPOCH: Epoch of observation (UT1 expressed as a Modified Julian Date).
621 * - RA: Right Ascension of source (radians, FK5 J2000).
622 * - DEC: Declination of source (radians, FK5 J2000).
623 * - VOFF: Velocity of the user-defined reference frame, towards the
624 * position given by RA and DEC, measured in the heliocentric
625 * reference frame.
626 *
627 * If the SpecMap is 3-dimensional, source positions are provided by the
628 * values supplied to inputs 2 and 3 of the SpecMap (which are simply
629 * copied to outputs 2 and 3). Note, usable values are still required
630 * for the RA and DEC arguments in order to define the "user-defined"
631 * reference frame used by USF2HL and HLF2US. However, AST__BAD can be
632 * supplied for RA and DEC if the user-defined reference frame is not
633 * required.
634
635 * Notes:
636 * - The specified conversion is appended only if the SpecMap's
637 * Invert attribute is zero. If it is non-zero, this function
638 * effectively prefixes the inverse of the conversion specified
639 * instead.
640 */
641
642 /* Local Variables: */
643 const char *argdesc[ MAX_ARGS ]; /* Pointers to argument descriptions */
644 const char *comment; /* Pointer to comment string */
645 const char *cvt_string; /* Pointer to conversion type string */
646 int argdec; /* Index of DEC argument */
647 int argra; /* Index of RA argument */
648 int i; /* Argument index */
649 int nargs; /* Number of user-supplied arguments */
650 int ncvt; /* Number of coordinate conversions */
651 int szargs; /* Size of arguments array */
652
653 /* Check the global error status. */
654 if ( !astOK ) return;
655
656 /* Validate the coordinate conversion type and obtain the number of
657 required user-supplied arguments, and the size of the array in which
658 to put the user-supplied arguments (the array meay leave room after
659 the user-supplied arguments for various useful pre-calculated values). */
660 cvt_string = CvtString( cvttype, &comment, &argra, &argdec, &nargs,
661 &szargs, argdesc, status );
662
663 /* If the coordinate conversion type was not valid, then report an
664 error. */
665 if ( astOK && !cvt_string ) {
666 astError( AST__SPCIN, "AddSpecCvt(%s): Invalid spectral coordinate "
667 "conversion type (%d).", status, astGetClass( this ),
668 (int) cvttype );
669 }
670
671 /* Note the number of coordinate conversions already stored in the SpecMap. */
672 if ( astOK ) {
673 ncvt = this->ncvt;
674
675 /* Extend the array of conversion types and the array of pointers to
676 their argument lists to accommodate the new one. */
677 this->cvttype = (int *) astGrow( this->cvttype, ncvt + 1,
678 sizeof( int ) );
679 this->cvtargs = (double **) astGrow( this->cvtargs, ncvt + 1,
680 sizeof( double * ) );
681
682 /* If OK, allocate memory and store a copy of the argument list,
683 putting a pointer to the copy into the SpecMap. */
684 if ( astOK ) {
685 this->cvtargs[ ncvt ] = astStore( NULL, args,
686 sizeof( double ) * (size_t) szargs );
687 }
688
689 /* Store the conversion type and increment the conversion count. Also put
690 AST__BAD in any elements of the argument array which are beyond the
691 end of the user-supplied arguments. These will be used to hold
692 intermediate values calculated on the basis of the user-supplied
693 arguments. */
694 if ( astOK ) {
695 this->cvttype[ ncvt ] = cvttype;
696 this->ncvt++;
697 for( i = nargs; i < szargs; i++ ) this->cvtargs[ ncvt ][ i ] = AST__BAD;
698 }
699 }
700 }
701
BaryVel(double ra,double dec,FrameDef * def,int * status)702 static double BaryVel( double ra, double dec, FrameDef *def, int *status ) {
703 /*
704 * Name:
705 * BaryVel
706
707 * Purpose:
708 * Find the velocity of the earth-sun barycentre away from the source.
709
710 * Type:
711 * Private function.
712
713 * Synopsis:
714 * #include "specmap.h"
715 * double BaryVel( double ra, double dec, FrameDef *def, int *status )
716
717 * Class Membership:
718 * SpecMap method.
719
720 * Description:
721 * This function finds the component of the velocity of the earth-sun
722 * barycentre away from a specified source position, at a given epoch, in
723 * the frame of rest of the centre of the Sun.
724
725 * Parameters:
726 * ra
727 * The RA (rads, FK5 J2000) of the source.
728 * dec
729 * The Dec (rads, FK5 J2000) of the source.
730 * def
731 * Pointer to a FrameDef structure which holds the parameters which
732 * define the frame, together with cached intermediate results.
733 * status
734 * Pointer to the inherited status variable.
735
736 * Returns:
737 * The component of the frame's velocity away from the position given by
738 * "ra" and "dec", in m/s, measured within the Heliographic frame of
739 * rest. Zero is returned if an error has already occurred.
740
741 */
742
743 /* Local Variables: */
744 double dpb[ 3 ]; /* Barycentric earth position vector */
745 double dph[ 3 ]; /* Heliocentric earth position vector */
746 double dvh[ 3 ]; /* Heliocentric earth velocity vector */
747 double v[ 3 ]; /* Source direction vector */
748
749 /* Check the global error status. */
750 if ( !astOK ) return 0.0;
751
752 /* Get the Cartesian vector towards the source, in the Cartesian FK5
753 J2000 system. */
754 palDcs2c( ra, dec, v );
755
756 /* If not already done so, get the Earth/Sun velocity and position vectors in
757 the same system. Speed is returned in units of AU/s. Store in the supplied
758 frame definition structure. */
759 if( def->dvb[ 0 ] == AST__BAD ) {
760 palEvp( def->epoch, 2000.0, def->dvb, dpb, dvh, dph );
761
762 /* Change the barycentric velocity of the earth into the heliocentric
763 velocity of the barycentre. */
764 def->dvb[ 0 ] = dvh[ 0 ] - def->dvb[ 0 ];
765 def->dvb[ 1 ] = dvh[ 1 ] - def->dvb[ 1 ];
766 def->dvb[ 2 ] = dvh[ 2 ] - def->dvb[ 2 ];
767 }
768
769 /* Return the component away from the source, of the velocity of the
770 barycentre relative to the sun (in m/s). */
771 return -palDvdv( v, def->dvb )*149.597870E9;
772
773 }
774
CvtCode(const char * cvt_string,int * status)775 static int CvtCode( const char *cvt_string, int *status ) {
776 /*
777 * Name:
778 * CvtCode
779
780 * Purpose:
781 * Convert a conversion type from a string representation to a code value.
782
783 * Type:
784 * Private function.
785
786 * Synopsis:
787 * #include "specmap.h"
788 * int CvtCode( const char *cvt_string, int *status )
789
790 * Class Membership:
791 * SpecMap member function.
792
793 * Description:
794 * This function accepts a string used to repersent one of the
795 * SpecMap coordinate conversions and converts it into a code
796 * value for internal use.
797
798 * Parameters:
799 * cvt_string
800 * Pointer to a constant null-terminated string representing a
801 * spectral coordinate conversion. This is case sensitive and should
802 * contain no unnecessary white space.
803 * status
804 * Pointer to the inherited status variable.
805
806 * Returned Value:
807 * The equivalent conversion code. If the string was not
808 * recognised, the code AST__SPEC_NULL is returned, without error.
809
810 * Notes:
811 * - A value of AST__SPEC_NULL will be returned if this function is
812 * invoked with the global error status set, or if it should fail
813 * for any reason.
814 */
815
816 /* Local Variables: */
817 int result; /* Result value to return */
818
819 /* Initialise. */
820 result = AST__SPEC_NULL;
821
822 /* Check the global error status. */
823 if ( !astOK ) return result;
824
825 /* Test the string against each recognised value in turn and assign
826 the result. */
827 if ( astChrMatch( cvt_string, "FRTOVL" ) ) {
828 result = AST__FRTOVL;
829
830 } else if ( astChrMatch( cvt_string, "VLTOFR" ) ) {
831 result = AST__VLTOFR;
832
833 } else if ( astChrMatch( cvt_string, "VLTOFR" ) ) {
834 result = AST__VLTOFR;
835
836 } else if ( astChrMatch( cvt_string, "ENTOFR" ) ) {
837 result = AST__ENTOFR;
838
839 } else if ( astChrMatch( cvt_string, "FRTOEN" ) ) {
840 result = AST__FRTOEN;
841
842 } else if ( astChrMatch( cvt_string, "WNTOFR" ) ) {
843 result = AST__WNTOFR;
844
845 } else if ( astChrMatch( cvt_string, "FRTOWN" ) ) {
846 result = AST__FRTOWN;
847
848 } else if ( astChrMatch( cvt_string, "WVTOFR" ) ) {
849 result = AST__WVTOFR;
850
851 } else if ( astChrMatch( cvt_string, "FRTOWV" ) ) {
852 result = AST__FRTOWV;
853
854 } else if ( astChrMatch( cvt_string, "AWTOFR" ) ) {
855 result = AST__AWTOFR;
856
857 } else if ( astChrMatch( cvt_string, "FRTOAW" ) ) {
858 result = AST__FRTOAW;
859
860 } else if ( astChrMatch( cvt_string, "VRTOVL" ) ) {
861 result = AST__VRTOVL;
862
863 } else if ( astChrMatch( cvt_string, "VLTOVR" ) ) {
864 result = AST__VLTOVR;
865
866 } else if ( astChrMatch( cvt_string, "VOTOVL" ) ) {
867 result = AST__VOTOVL;
868
869 } else if ( astChrMatch( cvt_string, "VLTOVO" ) ) {
870 result = AST__VLTOVO;
871
872 } else if ( astChrMatch( cvt_string, "ZOTOVL" ) ) {
873 result = AST__ZOTOVL;
874
875 } else if ( astChrMatch( cvt_string, "VLTOZO" ) ) {
876 result = AST__VLTOZO;
877
878 } else if ( astChrMatch( cvt_string, "BTTOVL" ) ) {
879 result = AST__BTTOVL;
880
881 } else if ( astChrMatch( cvt_string, "VLTOBT" ) ) {
882 result = AST__VLTOBT;
883
884 } else if ( astChrMatch( cvt_string, "USF2HL" ) ) {
885 result = AST__USF2HL;
886
887 } else if ( astChrMatch( cvt_string, "HLF2US" ) ) {
888 result = AST__HLF2US;
889
890 } else if ( astChrMatch( cvt_string, "TPF2HL" ) ) {
891 result = AST__TPF2HL;
892
893 } else if ( astChrMatch( cvt_string, "HLF2TP" ) ) {
894 result = AST__HLF2TP;
895
896 } else if ( astChrMatch( cvt_string, "GEF2HL" ) ) {
897 result = AST__GEF2HL;
898
899 } else if ( astChrMatch( cvt_string, "HLF2GE" ) ) {
900 result = AST__HLF2GE;
901
902 } else if ( astChrMatch( cvt_string, "BYF2HL" ) ) {
903 result = AST__BYF2HL;
904
905 } else if ( astChrMatch( cvt_string, "HLF2BY" ) ) {
906 result = AST__HLF2BY;
907
908 } else if ( astChrMatch( cvt_string, "LKF2HL" ) ) {
909 result = AST__LKF2HL;
910
911 } else if ( astChrMatch( cvt_string, "HLF2LK" ) ) {
912 result = AST__HLF2LK;
913
914 } else if ( astChrMatch( cvt_string, "LDF2HL" ) ) {
915 result = AST__LDF2HL;
916
917 } else if ( astChrMatch( cvt_string, "HLF2LD" ) ) {
918 result = AST__HLF2LD;
919
920 } else if ( astChrMatch( cvt_string, "LGF2HL" ) ) {
921 result = AST__LGF2HL;
922
923 } else if ( astChrMatch( cvt_string, "HLF2LG" ) ) {
924 result = AST__HLF2LG;
925
926 } else if ( astChrMatch( cvt_string, "GLF2HL" ) ) {
927 result = AST__GLF2HL;
928
929 } else if ( astChrMatch( cvt_string, "HLF2GL" ) ) {
930 result = AST__HLF2GL;
931
932 }
933
934 /* Return the result. */
935 return result;
936 }
937
CvtString(int cvt_code,const char ** comment,int * argra,int * argdec,int * nargs,int * szargs,const char * arg[MAX_ARGS],int * status)938 static const char *CvtString( int cvt_code, const char **comment,
939 int *argra, int *argdec, int *nargs, int *szargs,
940 const char *arg[ MAX_ARGS ], int *status ) {
941 /*
942 * Name:
943 * CvtString
944
945 * Purpose:
946 * Convert a conversion type from a code value to a string representation.
947
948 * Type:
949 * Private function.
950
951 * Synopsis:
952 * #include "specmap.h"
953 * const char *CvtString( int cvt_code, const char **comment,
954 * int *argra, int *argdec, int *nargs,
955 * int *szargs, const char *arg[ MAX_ARGS ], int *status )
956
957 * Class Membership:
958 * SpecMap member function.
959
960 * Description:
961 * This function accepts a code value used to represent one of the
962 * SpecMap coordinate conversions and converts it into an
963 * equivalent string representation. It also returns a descriptive
964 * comment and information about the arguments required in order to
965 * perform the conversion.
966
967 * Parameters:
968 * cvt_code
969 * The conversion code.
970 * comment
971 * Address of a location to return a pointer to a constant
972 * null-terminated string containing a description of the
973 * conversion.
974 * argra
975 * Address of an int in which to return the index of the argument
976 * corresponding to the source RA. Returned equal to -1 if the
977 * conversion does not have a source RA argument.
978 * argdec
979 * Address of an int in which to return the index of the argument
980 * corresponding to the source DEC. Returned equal to -1 if the
981 * conversion does not have a source DEC argument.
982 * nargs
983 * Address of an int in which to return the number of arguments
984 * required from the user in order to perform the conversion (may
985 * be zero).
986 * szargs
987 * Address of an int in which to return the number of arguments
988 * associated with the conversion. This may be bigger than "nargs"
989 * if the conversion can pre-calculate useful values on the basis
990 * of the user-supplied values. Such precalculated values are
991 * stored after the last user-supplied argument.
992 * arg
993 * An array in which to return a pointer to a constant
994 * null-terminated string for each argument (above) containing a
995 * description of what each argument represents. This includes both
996 * user-supplied arguments and pre-calculated values.
997 * status
998 * Pointer to the inherited status variable.
999
1000 * Returned Value:
1001 * Pointer to a constant null-terminated string representation of
1002 * the conversion code value supplied. If the code supplied is not
1003 * valid, a NULL pointer will be returned, without error.
1004
1005 * Notes:
1006 * - A NULL pointer value will be returned if this function is
1007 * invoked with the global error status set, or if it should fail
1008 * for any reason.
1009 */
1010
1011 /* Local Variables: */
1012 const char *result; /* Result pointer to return */
1013
1014 /* Initialise the returned values. */
1015 *comment = NULL;
1016 *nargs = 0;
1017 *argra = -1;
1018 *argdec = -1;
1019 result = NULL;
1020
1021 /* Check the global error status. */
1022 if ( !astOK ) return result;
1023
1024 /* Test for each valid code value in turn and assign the appropriate
1025 return values. */
1026 switch ( cvt_code ) {
1027
1028 case AST__FRTOVL:
1029 *comment = "Convert frequency to rel. velocity";
1030 result = "FRTOVL";
1031 *nargs = 1;
1032 *szargs = 1;
1033 arg[ 0 ] = "Rest frequency (Hz)";
1034 break;
1035
1036 case AST__VLTOFR:
1037 *comment = "Convert rel. velocity to frequency";
1038 result = "VLTOFR";
1039 *nargs = 1;
1040 *szargs = 1;
1041 arg[ 0 ] = "Rest frequency (Hz)";
1042 break;
1043
1044 case AST__ENTOFR:
1045 *comment = "Convert energy to frequency";
1046 result = "ENTOFR";
1047 *nargs = 0;
1048 *szargs = 0;
1049 break;
1050
1051 case AST__FRTOEN:
1052 *comment = "Convert frequency to energy";
1053 result = "FRTOEN";
1054 *nargs = 0;
1055 *szargs = 0;
1056 break;
1057
1058 case AST__WNTOFR:
1059 *comment = "Convert wave number to frequency";
1060 result = "WNTOFR";
1061 *nargs = 0;
1062 *szargs = 0;
1063 break;
1064
1065 case AST__FRTOWN:
1066 *comment = "Convert frequency to wave number";
1067 result = "FRTOWN";
1068 *nargs = 0;
1069 *szargs = 0;
1070 break;
1071
1072 case AST__WVTOFR:
1073 *comment = "Convert wavelength (vacuum) to frequency";
1074 result = "WVTOFR";
1075 *nargs = 0;
1076 *szargs = 0;
1077 break;
1078
1079 case AST__FRTOWV:
1080 *comment = "Convert frequency to wavelength (vacuum)";
1081 result = "FRTOWV";
1082 *nargs = 0;
1083 *szargs = 0;
1084 break;
1085
1086 case AST__AWTOFR:
1087 *comment = "Convert wavelength (air) to frequency";
1088 result = "AWTOFR";
1089 *nargs = 0;
1090 *szargs = 0;
1091 break;
1092
1093 case AST__FRTOAW:
1094 *comment = "Convert frequency to wavelength (air)";
1095 result = "FRTOAW";
1096 *nargs = 0;
1097 *szargs = 0;
1098 break;
1099
1100 case AST__VRTOVL:
1101 *comment = "Convert radio to rel. velocity";
1102 result = "VRTOVL";
1103 *nargs = 0;
1104 *szargs = 0;
1105 break;
1106
1107 case AST__VLTOVR:
1108 *comment = "Convert relativistic to radio velocity";
1109 result = "VLTOVR";
1110 *nargs = 0;
1111 *szargs = 0;
1112 break;
1113
1114 case AST__VOTOVL:
1115 *comment = "Convert optical to rel. velocity";
1116 result = "VOTOVL";
1117 *nargs = 0;
1118 *szargs = 0;
1119 break;
1120
1121 case AST__VLTOVO:
1122 *comment = "Convert relativistic to optical velocity";
1123 result = "VLTOVO";
1124 *nargs = 0;
1125 *szargs = 0;
1126 break;
1127
1128 case AST__ZOTOVL:
1129 *comment = "Convert redshift to rel. velocity";
1130 result = "ZOTOVL";
1131 *nargs = 0;
1132 *szargs = 0;
1133 break;
1134
1135 case AST__VLTOZO:
1136 *comment = "Convert rel. velocity to redshift";
1137 result = "VLTOZO";
1138 *nargs = 0;
1139 *szargs = 0;
1140 break;
1141
1142 case AST__BTTOVL:
1143 *comment = "Convert beta factor to rel. velocity";
1144 result = "BTTOVL";
1145 *nargs = 0;
1146 *szargs = 0;
1147 break;
1148
1149 case AST__VLTOBT:
1150 *comment = "Convert rel. velocity to beta factor";
1151 result = "VLTOBT";
1152 *nargs = 0;
1153 *szargs = 0;
1154 break;
1155
1156 case AST__USF2HL:
1157 *comment = "Convert from user-defined to heliocentric frequency";
1158 result = "USF2HL";
1159 *argra = 1;
1160 *argdec = 2;
1161 *nargs = 3;
1162 *szargs = 4;
1163 arg[ 0 ] = "Velocity offset (m/s)";
1164 arg[ 1 ] = "RA of source (FK5 J2000, radians)";
1165 arg[ 2 ] = "DEC of source (FK5 J2000, radians)";
1166 arg[ 3 ] = "Frequency correction factor";
1167 break;
1168
1169 case AST__HLF2US:
1170 *comment = "Convert from heliocentric to user-defined frequency";
1171 result = "HLF2US";
1172 *argra = 1;
1173 *argdec = 2;
1174 *nargs = 3;
1175 *szargs = 4;
1176 arg[ 0 ] = "Velocity offset (m/s)";
1177 arg[ 1 ] = "RA of source (FK5 J2000, radians)";
1178 arg[ 2 ] = "DEC of source (FK5 J2000, radians)";
1179 arg[ 3 ] = "Frequency correction factor";
1180 break;
1181
1182 case AST__TPF2HL:
1183 *comment = "Convert from Topocentric to heliocentric frequency";
1184 result = "TPF2HL";
1185 *argra = 4;
1186 *argdec = 5;
1187 *nargs = 6;
1188 *szargs = 7;
1189 arg[ 0 ] = "Longitude (positive eastwards, radians)";
1190 arg[ 1 ] = "Latitude (geodetic, radians)";
1191 arg[ 2 ] = "Altitude (geodetic, metres)";
1192 arg[ 3 ] = "UT1 epoch of observaton (Modified Julian Date)";
1193 arg[ 4 ] = "RA of source (FK5 J2000, radians)";
1194 arg[ 5 ] = "DEC of source (FK5 J2000, radians)";
1195 arg[ 6 ] = "Frequency correction factor";
1196 break;
1197
1198 case AST__HLF2TP:
1199 *comment = "Convert from Heliocentric to topocentric frequency";
1200 result = "HLF2TP";
1201 *argra = 4;
1202 *argdec = 5;
1203 *nargs = 6;
1204 *szargs = 7;
1205 arg[ 0 ] = "Longitude (positive eastwards, radians)";
1206 arg[ 1 ] = "Latitude (geodetic, radians)";
1207 arg[ 2 ] = "Altitude (geodetic, metres)";
1208 arg[ 3 ] = "UT1 epoch of observaton (Modified Julian Date)";
1209 arg[ 4 ] = "RA of source (FK5 J2000, radians)";
1210 arg[ 5 ] = "DEC of source (FK5 J2000, radians)";
1211 arg[ 6 ] = "Frequency correction factor";
1212 break;
1213
1214 case AST__GEF2HL:
1215 *comment = "Convert from Geocentric to heliocentric frequency";
1216 result = "GEF2HL";
1217 *argra = 1;
1218 *argdec = 2;
1219 *nargs = 3;
1220 *szargs = 4;
1221 arg[ 0 ] = "UT1 epoch of observaton (Modified Julian Date)";
1222 arg[ 1 ] = "RA of source (FK5 J2000, radians)";
1223 arg[ 2 ] = "DEC of source (FK5 J2000, radians)";
1224 arg[ 3 ] = "Frequency correction factor";
1225 break;
1226
1227 case AST__HLF2GE:
1228 *comment = "Convert from Heliocentric to geocentric frequency";
1229 result = "HLF2GE";
1230 *argra = 1;
1231 *argdec = 2;
1232 *nargs = 3;
1233 *szargs = 4;
1234 arg[ 0 ] = "UT1 epoch of observaton (Modified Julian Date)";
1235 arg[ 1 ] = "RA of source (FK5 J2000, radians)";
1236 arg[ 2 ] = "DEC of source (FK5 J2000, radians)";
1237 arg[ 3 ] = "Frequency correction factor";
1238 break;
1239
1240 case AST__BYF2HL:
1241 *comment = "Convert from Barycentric to heliocentric frequency";
1242 result = "BYF2HL";
1243 *argra = 1;
1244 *argdec = 2;
1245 *nargs = 3;
1246 *szargs = 4;
1247 arg[ 0 ] = "UT1 epoch of observaton (Modified Julian Date)";
1248 arg[ 1 ] = "RA of source (FK5 J2000, radians)";
1249 arg[ 2 ] = "DEC of source (FK5 J2000, radians)";
1250 arg[ 3 ] = "Frequency correction factor";
1251 break;
1252
1253 case AST__HLF2BY:
1254 *comment = "Convert from Heliocentric to barycentric frequency";
1255 result = "HLF2BY";
1256 *argra = 1;
1257 *argdec = 2;
1258 *nargs = 3;
1259 *szargs = 4;
1260 arg[ 0 ] = "UT1 epoch of observaton (Modified Julian Date)";
1261 arg[ 1 ] = "RA of source (FK5 J2000, radians)";
1262 arg[ 2 ] = "DEC of source (FK5 J2000, radians)";
1263 arg[ 3 ] = "Frequency correction factor";
1264 break;
1265
1266 case AST__LKF2HL:
1267 *comment = "Convert from LSRK to heliocentric frequency";
1268 result = "LKF2HL";
1269 *argra = 0;
1270 *argdec = 1;
1271 *nargs = 2;
1272 *szargs = 3;
1273 arg[ 0 ] = "RA of source (FK5 J2000, radians)";
1274 arg[ 1 ] = "DEC of source (FK5 J2000, radians)";
1275 arg[ 2 ] = "Frequency correction factor";
1276 break;
1277
1278 case AST__HLF2LK:
1279 *comment = "Convert from Heliocentric to LSRK frequency";
1280 result = "HLF2LK";
1281 *argra = 0;
1282 *argdec = 1;
1283 *nargs = 2;
1284 *szargs = 3;
1285 arg[ 0 ] = "RA of source (FK5 J2000, radians)";
1286 arg[ 1 ] = "DEC of source (FK5 J2000, radians)";
1287 arg[ 2 ] = "Frequency correction factor";
1288 break;
1289
1290 case AST__LDF2HL:
1291 *comment = "Convert from LSRD to heliocentric frequency";
1292 result = "LDF2HL";
1293 *argra = 0;
1294 *argdec = 1;
1295 *nargs = 2;
1296 *szargs = 3;
1297 arg[ 0 ] = "RA of source (FK5 J2000, radians)";
1298 arg[ 1 ] = "DEC of source (FK5 J2000, radians)";
1299 arg[ 2 ] = "Frequency correction factor";
1300 break;
1301
1302 case AST__HLF2LD:
1303 *comment = "Convert from Heliocentric to LSRD frequency";
1304 result = "HLF2LD";
1305 *argra = 0;
1306 *argdec = 1;
1307 *nargs = 2;
1308 *szargs = 3;
1309 arg[ 0 ] = "RA of source (FK5 J2000, radians)";
1310 arg[ 1 ] = "DEC of source (FK5 J2000, radians)";
1311 arg[ 2 ] = "Frequency correction factor";
1312 break;
1313
1314 case AST__LGF2HL:
1315 *comment = "Convert from Local group to heliocentric frequency";
1316 result = "LGF2HL";
1317 *argra = 0;
1318 *argdec = 1;
1319 *nargs = 2;
1320 *szargs = 3;
1321 arg[ 0 ] = "RA of source (FK5 J2000, radians)";
1322 arg[ 1 ] = "DEC of source (FK5 J2000, radians)";
1323 arg[ 2 ] = "Frequency correction factor";
1324 break;
1325
1326 case AST__HLF2LG:
1327 *comment = "Convert from Heliocentric to local group frequency";
1328 result = "HLF2LG";
1329 *argra = 0;
1330 *argdec = 1;
1331 *nargs = 2;
1332 *szargs = 3;
1333 arg[ 0 ] = "RA of source (FK5 J2000, radians)";
1334 arg[ 1 ] = "DEC of source (FK5 J2000, radians)";
1335 arg[ 2 ] = "Frequency correction factor";
1336 break;
1337
1338 case AST__GLF2HL:
1339 *comment = "Convert from Galactic to heliocentric frequency";
1340 result = "GLF2HL";
1341 *argra = 0;
1342 *argdec = 1;
1343 *nargs = 2;
1344 *szargs = 3;
1345 arg[ 0 ] = "RA of source (FK5 J2000, radians)";
1346 arg[ 1 ] = "DEC of source (FK5 J2000, radians)";
1347 arg[ 2 ] = "Frequency correction factor";
1348 break;
1349
1350 case AST__HLF2GL:
1351 *comment = "Convert from Heliocentric to galactic frequency";
1352 *argra = 0;
1353 *argdec = 1;
1354 result = "HLF2GL";
1355 *nargs = 2;
1356 *szargs = 3;
1357 arg[ 0 ] = "RA of source (FK5 J2000, radians)";
1358 arg[ 1 ] = "DEC of source (FK5 J2000, radians)";
1359 arg[ 2 ] = "Frequency correction factor";
1360 break;
1361
1362 }
1363
1364 /* Return the result. */
1365 return result;
1366 }
1367
FrameChange(int cvt_code,int np,double * ra,double * dec,double * freq,double * args,int forward,int * status)1368 static int FrameChange( int cvt_code, int np, double *ra, double *dec, double *freq,
1369 double *args, int forward, int *status ){
1370 /*
1371 * Name:
1372 * FrameChange
1373
1374 * Purpose:
1375 * Apply a doppler shift caused by a change of reference frame.
1376
1377 * Type:
1378 * Private function.
1379
1380 * Synopsis:
1381 * #include "specmap.h"
1382 * int FrameChange( int cvt_code, int np, double *ra, double *dec,
1383 * double *freq, double *args, int forward, int *status )
1384
1385 * Class Membership:
1386 * SpecMap method.
1387
1388 * Description:
1389 * This function modifies the supplied frequency values in order to
1390 * apply a doppler shift caused by a change of the observers rest-frame.
1391
1392 * Parameters:
1393 * cvt_code
1394 * A code indicating the conversion to be applied. If the code does
1395 * not correspond to a change of rest-frame, then the supplied
1396 * frequencies are left unchanged and zero is returned as the
1397 * function value.
1398 * np
1399 * The number of frequency values to transform.
1400 * ra
1401 * Pointer to an array of "np" RA (J2000 FK5) values at which the
1402 * "np" frequencies are observed. These are unchanged on exit. If a
1403 * NULL pointer is supplied, then all frequencies are assumed to be
1404 * observed at the single RA value given by "refra"
1405 * dec
1406 * Pointer to an array of "np" Dec (J2000 FK5) values at which the
1407 * "np" frequencies are observed. These are unchanged on exit. If a
1408 * NULL pointer is supplied, then all frequencies are assumed to be
1409 * observed at the single Dec value given by "refdec"
1410 * freq
1411 * Pointer to an array of "np" frequency values, measured in the
1412 * input rest-frame. These are modified on return to hold the
1413 * corresponding values measured in the output rest-frame.
1414 * args
1415 * Pointer to an array holding the conversion arguments. The number
1416 * of arguments expected depends on the particular conversion being
1417 * used.
1418 * forward
1419 * Should the conversion be applied in the forward or inverse
1420 * direction? Non-zero for forward, zero for inverse.
1421 * status
1422 * Pointer to the inherited status variable.
1423
1424 * Returned Value:
1425 * Non-zero if the supplied conversion code corresponds to a change of
1426 * reference frame. Zoer otherwise (in which case the upplied values
1427 * will not have been changed).
1428
1429 * Notes:
1430 * - The "args" array contains RA and DEC values which give the "source"
1431 * position (FK5 J2000). If a NULL value is supplied for the "ra"
1432 * parameter, then these args define the position of all the frequency
1433 * values. In addition they also define the direction of motion of
1434 * the "user-defined" rest-frame (see "veluser"). Thus they should still
1435 * be supplied even if "ra" is NULL.
1436
1437 */
1438
1439 /* Local Variables: */
1440 FrameDef def; /* Structure holding frame parameters */
1441 double (* cvtFunc)( double, double, FrameDef *, int * ); /* Pointer to conversion function */
1442 double *fcorr; /* Pointer to frequency correction factor */
1443 double *pdec; /* Pointer to next Dec value */
1444 double *pf; /* Pointer to next frequency value */
1445 double *pra; /* Pointer to next RA value */
1446 double factor; /* Frequency correction factor */
1447 double s; /* Velocity correction (m/s) */
1448 int i; /* Loop index */
1449 int result; /* Returned value */
1450 int sign; /* Sign for velocity correction */
1451
1452 /* Check inherited status. */
1453 if( !astOK ) return 0;
1454
1455 /* Initialise */
1456 cvtFunc = NULL;
1457 fcorr = NULL;
1458 sign = 0;
1459
1460 /* Set the return value to indicate that the supplied conversion code
1461 represents a change of rest-frame. */
1462 result = 1;
1463
1464 /* Initialise a structure which stores parameters which define the
1465 transformation. */
1466 def.obsalt = AST__BAD;
1467 def.obslat = AST__BAD;
1468 def.obslon = AST__BAD;
1469 def.epoch = AST__BAD;
1470 def.refdec = AST__BAD;
1471 def.refra = AST__BAD;
1472 def.veluser = AST__BAD;
1473 def.last = AST__BAD;
1474 def.amprms[ 0 ] = AST__BAD;
1475 def.vuser[ 0 ] = AST__BAD;
1476 def.dvh[ 0 ] = AST__BAD;
1477 def.dvb[ 0 ] = AST__BAD;
1478
1479 /* Test for each rest-frame code value in turn and assign the appropriate
1480 values. */
1481 switch ( cvt_code ) {
1482
1483 case AST__USF2HL:
1484 cvtFunc = UserVel;
1485 def.veluser = args[ 0 ];
1486 def.refra = args[ 1 ];
1487 def.refdec = args[ 2 ];
1488 fcorr = args + 3;
1489 sign = -1;
1490 break;
1491
1492 case AST__HLF2US:
1493 cvtFunc = UserVel;
1494 def.veluser = args[ 0 ];
1495 def.refra = args[ 1 ];
1496 def.refdec = args[ 2 ];
1497 fcorr = args + 3;
1498 sign = +1;
1499 break;
1500
1501 case AST__TPF2HL:
1502 cvtFunc = TopoVel;
1503 def.obslon = args[ 0 ];
1504 def.obslat = args[ 1 ];
1505 def.obsalt = args[ 2 ];
1506 def.epoch = args[ 3 ];
1507 def.refra = args[ 4 ];
1508 def.refdec = args[ 5 ];
1509 fcorr = args + 6;
1510 sign = -1;
1511 break;
1512
1513 case AST__HLF2TP:
1514 cvtFunc = TopoVel;
1515 def.obslon = args[ 0 ];
1516 def.obslat = args[ 1 ];
1517 def.obsalt = args[ 2 ];
1518 def.epoch = args[ 3 ];
1519 def.refra = args[ 4 ];
1520 def.refdec = args[ 5 ];
1521 fcorr = args + 6;
1522 sign = +1;
1523 break;
1524
1525 case AST__GEF2HL:
1526 cvtFunc = GeoVel;
1527 def.epoch = args[ 0 ];
1528 def.refra = args[ 1 ];
1529 def.refdec = args[ 2 ];
1530 fcorr = args + 3;
1531 sign = -1;
1532 break;
1533
1534 case AST__HLF2GE:
1535 cvtFunc = GeoVel;
1536 def.epoch = args[ 0 ];
1537 def.refra = args[ 1 ];
1538 def.refdec = args[ 2 ];
1539 fcorr = args + 3;
1540 sign = +1;
1541 break;
1542
1543 case AST__BYF2HL:
1544 cvtFunc = BaryVel;
1545 def.epoch = args[ 0 ];
1546 def.refra = args[ 1 ];
1547 def.refdec = args[ 2 ];
1548 fcorr = args + 3;
1549 sign = -1;
1550 break;
1551
1552 case AST__HLF2BY:
1553 cvtFunc = BaryVel;
1554 def.epoch = args[ 0 ];
1555 def.refra = args[ 1 ];
1556 def.refdec = args[ 2 ];
1557 fcorr = args + 3;
1558 sign = +1;
1559 break;
1560
1561 case AST__LKF2HL:
1562 cvtFunc = LsrkVel;
1563 def.refra = args[ 0 ];
1564 def.refdec = args[ 1 ];
1565 fcorr = args + 2;
1566 sign = -1;
1567 break;
1568
1569 case AST__HLF2LK:
1570 cvtFunc = LsrkVel;
1571 def.refra = args[ 0 ];
1572 def.refdec = args[ 1 ];
1573 fcorr = args + 2;
1574 sign = +1;
1575 break;
1576
1577 case AST__LDF2HL:
1578 cvtFunc = LsrdVel;
1579 def.refra = args[ 0 ];
1580 def.refdec = args[ 1 ];
1581 fcorr = args + 2;
1582 sign = -1;
1583 break;
1584
1585 case AST__HLF2LD:
1586 cvtFunc = LsrdVel;
1587 def.refra = args[ 0 ];
1588 def.refdec = args[ 1 ];
1589 fcorr = args + 2;
1590 sign = +1;
1591 break;
1592
1593 case AST__LGF2HL:
1594 cvtFunc = LgVel;
1595 def.refra = args[ 0 ];
1596 def.refdec = args[ 1 ];
1597 fcorr = args + 2;
1598 sign = -1;
1599 break;
1600
1601 case AST__HLF2LG:
1602 cvtFunc = LgVel;
1603 def.refra = args[ 0 ];
1604 def.refdec = args[ 1 ];
1605 fcorr = args + 2;
1606 sign = +1;
1607 break;
1608
1609 case AST__GLF2HL:
1610 cvtFunc = GalVel;
1611 def.refra = args[ 0 ];
1612 def.refdec = args[ 1 ];
1613 fcorr = args + 2;
1614 sign = -1;
1615 break;
1616
1617 case AST__HLF2GL:
1618 cvtFunc = GalVel;
1619 def.refra = args[ 0 ];
1620 def.refdec = args[ 1 ];
1621 fcorr = args + 2;
1622 sign = +1;
1623 break;
1624
1625 /* If the supplied code does not represent a change of rest-frame, clear
1626 the returned flag. */
1627 default:
1628 result = 0;
1629 }
1630
1631 /* Check we have a rest-frame code. */
1632 if( result ) {
1633
1634 /* First deal with cases where we have a single source position (given by
1635 refra and refdec). */
1636 if( !ra ) {
1637
1638 /* If the frequency correction factor has not been found, find it now. */
1639 if( *fcorr == AST__BAD ) {
1640
1641 /* Get the velocity correction. This is the component of the velocity of the
1642 output system, away from the source, as measured in the input system. */
1643 s = sign*cvtFunc( def.refra, def.refdec, &def, status );
1644
1645 /* Find the factor by which to correct supplied frequencies. If the
1646 velocity correction is positive, the output frequency wil be lower than
1647 the input frequency. */
1648 if( s < AST__C && s > -AST__C ) {
1649 *fcorr = sqrt( ( AST__C - s )/( AST__C + s ) );
1650 }
1651 }
1652
1653 /* Correct each supplied frequency. */
1654 if( *fcorr != AST__BAD && *fcorr != 0.0 ) {
1655 factor = forward ? *fcorr : 1.0 / ( *fcorr );
1656 pf = freq;
1657 for( i = 0; i < np; i++, pf++ ) {
1658 if( *pf != AST__BAD ) *pf *= factor;
1659 }
1660
1661 /* Set returned values bad if the velocity correction is un-physical. */
1662 } else {
1663 pf = freq;
1664 for( i = 0; i < np; i++ ) *(pf++) = AST__BAD;
1665 }
1666
1667 /* Now deal with cases where each frequency value has its own source
1668 position. */
1669 } else {
1670
1671 /* Invert the sign if we are doing a inverse transformation. */
1672 if( !forward ) sign = -sign;
1673
1674 /* Loop round each value. */
1675 pf = freq;
1676 pra = ra;
1677 pdec = dec;
1678 for( i = 0; i < np; i++ ) {
1679
1680 /* If the ra or dec is bad, store a bad frequency. */
1681 if( *pra == AST__BAD || *pdec == AST__BAD || *pf == AST__BAD ) {
1682 *pf = AST__BAD;
1683
1684 /* Otherwise, produce a corrected frequency. */
1685 } else {
1686
1687 /* Get the velocity correction. */
1688 s = sign*cvtFunc( *pra, *pdec, &def, status );
1689
1690 /* Correct this frequency, if possible. Otherwise set bad. */
1691 if( s < AST__C && s > -AST__C ) {
1692 *pf *= sqrt( ( AST__C - s )/( AST__C + s ) );
1693 } else {
1694 *pf = AST__BAD;
1695 }
1696 }
1697
1698 /* Move on to the next position. */
1699 pf++;
1700 pra++;
1701 pdec++;
1702 }
1703 }
1704 }
1705
1706 /* Return the result. */
1707 return result;
1708 }
1709
GalVel(double ra,double dec,FrameDef * def,int * status)1710 static double GalVel( double ra, double dec, FrameDef *def, int *status ) {
1711 /*
1712 * Name:
1713 * GalVel
1714
1715 * Purpose:
1716 * Find the velocity of the galactic centre away from the source.
1717
1718 * Type:
1719 * Private function.
1720
1721 * Synopsis:
1722 * #include "specmap.h"
1723 * double GalVel( double ra, double dec, FrameDef *def, int *status )
1724
1725 * Class Membership:
1726 * SpecMap method.
1727
1728 * Description:
1729 * This function finds the component of the velocity of the galactic
1730 * centre away from a specified source position, in the frame of rest
1731 * of the Sun.
1732
1733 * Parameters:
1734 * ra
1735 * The RA (rads, FK5 J2000) of the source.
1736 * dec
1737 * The Dec (rads, FK5 J2000) of the source.
1738 * def
1739 * Pointer to a FrameDef structure which holds the parameters which
1740 * define the frame, together with cached intermediate results.
1741 * status
1742 * Pointer to the inherited status variable.
1743
1744 * Returns:
1745 * The component of the frame's velocity away from the position given by
1746 * "ra" and "dec", in m/s, measured within the Heliographic frame of
1747 * rest. Zero is returned if an error has already occurred.
1748
1749 */
1750
1751 /* Local Variables: */
1752 double s1, s2;
1753
1754 /* Check the global error status. */
1755 if ( !astOK ) return 0.0;
1756
1757 /* Get the component away from the source, of the velocity of the sun
1758 relative to the dynamic LSR (in km/s). */
1759 s1 = (double) palRvlsrd( (float) ra, (float) dec );
1760
1761 /* Get the component away from the source, of the velocity of the
1762 dynamic LSR relative to the galactic centre (in km/s). */
1763 s2 = (double) palRvgalc( (float) ra, (float) dec );
1764
1765 /* Return the total velocity of the galactic centre away from the source,
1766 relative to the sun, in m/s. */
1767 return -1000.0*( s1 + s2 );
1768 }
1769
GeoVel(double ra,double dec,FrameDef * def,int * status)1770 static double GeoVel( double ra, double dec, FrameDef *def, int *status ) {
1771 /*
1772 * Name:
1773 * GeoVel
1774
1775 * Purpose:
1776 * Find the velocity of the earth away from the source.
1777
1778 * Type:
1779 * Private function.
1780
1781 * Synopsis:
1782 * #include "specmap.h"
1783 * double GeoVel( double ra, double dec, FrameDef *def, int *status )
1784
1785 * Class Membership:
1786 * SpecMap method.
1787
1788 * Description:
1789 * This function finds the component of the velocity of the earth away
1790 * from a specified source position, at a given epoch, in the frame of
1791 * rest of the Sun.
1792
1793 * Parameters:
1794 * ra
1795 * The RA (rads, FK5 J2000) of the source.
1796 * dec
1797 * The Dec (rads, FK5 J2000) of the source.
1798 * def
1799 * Pointer to a FrameDef structure which holds the parameters which
1800 * define the frame, together with cached intermediate results.
1801 * status
1802 * Pointer to the inherited status variable.
1803
1804 * Returns:
1805 * The component of the frame's velocity away from the position given by
1806 * "ra" and "dec", in m/s, measured within the Heliographic frame of
1807 * rest. Zero is returned if an error has already occurred.
1808
1809 */
1810
1811 /* Local Variables: */
1812 double dpb[ 3 ]; /* Barycentric earth position vector */
1813 double dph[ 3 ]; /* Heliocentric earth position vector */
1814 double dvb[ 3 ]; /* Barycentric earth velocity vector */
1815 double v[ 3 ]; /* Source direction vector */
1816
1817 /* Check the global error status. */
1818 if ( !astOK ) return 0.0;
1819
1820 /* Get the Cartesian vector towards the source, in the Cartesian FK5
1821 J2000 system. */
1822 palDcs2c( ra, dec, v );
1823
1824 /* If not already done so, get the Earth/Sun velocity and position vectors in
1825 the same system. Speed is returned in units of AU/s. Store in the supplied
1826 frame definition structure. */
1827 if( def->dvh[ 0 ] == AST__BAD ) palEvp( def->epoch, 2000.0, dvb, dpb,
1828 def->dvh, dph );
1829
1830 /* Return the component away from the source, of the velocity of the earths
1831 centre relative to the sun (in m/s). */
1832 return -palDvdv( v, def->dvh )*149.597870E9;
1833 }
1834
astInitSpecMapVtab_(AstSpecMapVtab * vtab,const char * name,int * status)1835 void astInitSpecMapVtab_( AstSpecMapVtab *vtab, const char *name, int *status ) {
1836 /*
1837 *+
1838 * Name:
1839 * astInitSpecMapVtab
1840
1841 * Purpose:
1842 * Initialise a virtual function table for a SpecMap.
1843
1844 * Type:
1845 * Protected function.
1846
1847 * Synopsis:
1848 * #include "specmap.h"
1849 * void astInitSpecMapVtab( AstSpecMapVtab *vtab, const char *name )
1850
1851 * Class Membership:
1852 * SpecMap vtab initialiser.
1853
1854 * Description:
1855 * This function initialises the component of a virtual function
1856 * table which is used by the SpecMap class.
1857
1858 * Parameters:
1859 * vtab
1860 * Pointer to the virtual function table. The components used by
1861 * all ancestral classes will be initialised if they have not already
1862 * been initialised.
1863 * name
1864 * Pointer to a constant null-terminated character string which contains
1865 * the name of the class to which the virtual function table belongs (it
1866 * is this pointer value that will subsequently be returned by the Object
1867 * astClass function).
1868 *-
1869 */
1870
1871 /* Local Variables: */
1872 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
1873 AstObjectVtab *object; /* Pointer to Object component of Vtab */
1874 AstMappingVtab *mapping; /* Pointer to Mapping component of Vtab */
1875
1876 /* Check the local error status. */
1877 if ( !astOK ) return;
1878
1879 /* Get a pointer to the thread specific global data structure. */
1880 astGET_GLOBALS(NULL);
1881
1882 /* Initialize the component of the virtual function table used by the
1883 parent class. */
1884 astInitMappingVtab( (AstMappingVtab *) vtab, name );
1885
1886 /* Store a unique "magic" value in the virtual function table. This
1887 will be used (by astIsASpecMap) to determine if an object belongs to
1888 this class. We can conveniently use the address of the (static)
1889 class_check variable to generate this unique value. */
1890 vtab->id.check = &class_check;
1891 vtab->id.parent = &(((AstMappingVtab *) vtab)->id);
1892
1893 /* Initialise member function pointers. */
1894 /* ------------------------------------ */
1895 /* Store pointers to the member functions (implemented here) that
1896 provide virtual methods for this class. */
1897 vtab->SpecAdd = SpecAdd;
1898
1899 /* Save the inherited pointers to methods that will be extended, and
1900 replace them with pointers to the new member functions. */
1901 object = (AstObjectVtab *) vtab;
1902 mapping = (AstMappingVtab *) vtab;
1903 parent_getobjsize = object->GetObjSize;
1904 object->GetObjSize = GetObjSize;
1905
1906 parent_transform = mapping->Transform;
1907 mapping->Transform = Transform;
1908
1909 parent_rate = mapping->Rate;
1910 mapping->Rate = Rate;
1911
1912 /* Store replacement pointers for methods which will be over-ridden by
1913 new member functions implemented here. */
1914 object->Equal = Equal;
1915 mapping->MapMerge = MapMerge;
1916
1917 /* Declare the copy constructor, destructor and class dump
1918 function. */
1919 astSetCopy( vtab, Copy );
1920 astSetDelete( vtab, Delete );
1921 astSetDump( vtab, Dump, "SpecMap",
1922 "Conversion between spectral coordinate systems" );
1923
1924 /* If we have just initialised the vtab for the current class, indicate
1925 that the vtab is now initialised, and store a pointer to the class
1926 identifier in the base "object" level of the vtab. */
1927 if( vtab == &class_vtab ) {
1928 class_init = 1;
1929 astSetVtabClassIdentifier( vtab, &(vtab->id) );
1930 }
1931 }
1932
LgVel(double ra,double dec,FrameDef * def,int * status)1933 static double LgVel( double ra, double dec, FrameDef *def, int *status ) {
1934 /*
1935 * Name:
1936 * LgVel
1937
1938 * Purpose:
1939 * Find the velocity of the Local Group away from the source.
1940
1941 * Type:
1942 * Private function.
1943
1944 * Synopsis:
1945 * #include "specmap.h"
1946 * double LgVel( double ra, double dec, FrameDef *def, int *status )
1947
1948 * Class Membership:
1949 * SpecMap method.
1950
1951 * Description:
1952 * This function finds the component of the Local Group velocity away
1953 * from a specified source position, in the frame of rest of the Sun.
1954
1955 * Parameters:
1956 * ra
1957 * The RA (rads, FK5 J2000) of the source.
1958 * dec
1959 * The Dec (rads, FK5 J2000) of the source.
1960 * def
1961 * Pointer to a FrameDef structure which holds the parameters which
1962 * define the frame, together with cached intermediate results.
1963 * status
1964 * Pointer to the inherited status variable.
1965
1966 * Returns:
1967 * The component of the frame's velocity away from the position given by
1968 * "ra" and "dec", in m/s, measured within the Heliographic frame of
1969 * rest. Zero is returned if an error has already occurred.
1970
1971 */
1972
1973 /* Return the component away from the source, of the velocity of the
1974 local group relative to the sun (in m/s). */
1975 return -1000.0*palRvlg( (float) ra, (float) dec );
1976 }
1977
LsrdVel(double ra,double dec,FrameDef * def,int * status)1978 static double LsrdVel( double ra, double dec, FrameDef *def, int *status ) {
1979 /*
1980 * Name:
1981 * LsrdVel
1982
1983 * Purpose:
1984 * Find the velocity of the Dynamical LSR away from the source.
1985
1986 * Type:
1987 * Private function.
1988
1989 * Synopsis:
1990 * #include "specmap.h"
1991 * double LsrdVel( double ra, double dec, FrameDef *def, int *status )
1992
1993 * Class Membership:
1994 * SpecMap method.
1995
1996 * Description:
1997 * This function finds the component of the velocity of the Dynamical
1998 * LSR away from a specified source position, in the frame of rest of
1999 * the Sun.
2000
2001 * Parameters:
2002 * ra
2003 * The RA (rads, FK5 J2000) of the source.
2004 * dec
2005 * The Dec (rads, FK5 J2000) of the source.
2006 * def
2007 * Pointer to a FrameDef structure which holds the parameters which
2008 * define the frame, together with cached intermediate results.
2009 * status
2010 * Pointer to the inherited status variable.
2011
2012 * Returns:
2013 * The component of the frame's velocity away from the position given by
2014 * "ra" and "dec", in m/s, measured within the Heliographic frame of
2015 * rest. Zero is returned if an error has already occurred.
2016
2017 */
2018 /* Check the global error status. */
2019 if ( !astOK ) return 0.0;
2020
2021 /* Get the component away from the source, of the velocity of the sun
2022 relative to the dynamical LSR (in m/s). This can also be thought of as the
2023 velocity of the LSR towards the source relative to the sun. Return the
2024 negated value (i.e. velocity of lsrd *away from* the source. */
2025 return -1000.0*palRvlsrd( (float) ra, (float) dec );
2026 }
2027
LsrkVel(double ra,double dec,FrameDef * def,int * status)2028 static double LsrkVel( double ra, double dec, FrameDef *def, int *status ) {
2029 /*
2030 * Name:
2031 * LsrkVel
2032
2033 * Purpose:
2034 * Find the velocity of the Kinematic LSR away from the source.
2035
2036 * Type:
2037 * Private function.
2038
2039 * Synopsis:
2040 * #include "specmap.h"
2041 * double LsrkVel( double ra, double dec, FrameDef *def, int *status )
2042
2043 * Class Membership:
2044 * SpecMap method.
2045
2046 * Description:
2047 * This function finds the component of the velocity of the Kinematic
2048 * LSR away from a specified source position, in the frame of rest of
2049 * the Sun.
2050
2051 * Parameters:
2052 * ra
2053 * The RA (rads, FK5 J2000) of the source.
2054 * dec
2055 * The Dec (rads, FK5 J2000) of the source.
2056 * def
2057 * Pointer to a FrameDef structure which holds the parameters which
2058 * define the frame, together with cached intermediate results.
2059 * status
2060 * Pointer to the inherited status variable.
2061
2062 * Returns:
2063 * The component of the frame's velocity away from the position given by
2064 * "ra" and "dec", in m/s, measured within the Heliographic frame of
2065 * rest. Zero is returned if an error has already occurred.
2066
2067 */
2068 /* Check the global error status. */
2069 if ( !astOK ) return 0.0;
2070
2071 /* Get the component away from the source, of the velocity of the sun
2072 relative to the kinematic LSR (in m/s). This can also be thought of as the
2073 velocity of the LSR towards the source relative to the sun. Return the
2074 negated value (i.e. velocity of lsrk *away from* the source. */
2075 return -1000.0*palRvlsrk( (float) ra, (float) dec );
2076 }
2077
MapMerge(AstMapping * this,int where,int series,int * nmap,AstMapping *** map_list,int ** invert_list,int * status)2078 static int MapMerge( AstMapping *this, int where, int series, int *nmap,
2079 AstMapping ***map_list, int **invert_list, int *status ) {
2080 /*
2081 * Name:
2082 * MapMerge
2083
2084 * Purpose:
2085 * Simplify a sequence of Mappings containing a SpecMap.
2086
2087 * Type:
2088 * Private function.
2089
2090 * Synopsis:
2091 * #include "mapping.h
2092 * int MapMerge( AstMapping *this, int where, int series, int *nmap,
2093 * AstMapping ***map_list, int **invert_list, int *status )
2094
2095 * Class Membership:
2096 * SpecMap method (over-rides the protected astMapMerge method
2097 * inherited from the Mapping class).
2098
2099 * Description:
2100 * This function attempts to simplify a sequence of Mappings by
2101 * merging a nominated SpecMap in the sequence with its neighbours,
2102 * so as to shorten the sequence if possible.
2103 *
2104 * In many cases, simplification will not be possible and the
2105 * function will return -1 to indicate this, without further
2106 * action.
2107 *
2108 * In most cases of interest, however, this function will either
2109 * attempt to replace the nominated SpecMap with one which it
2110 * considers simpler, or to merge it with the Mappings which
2111 * immediately precede it or follow it in the sequence (both will
2112 * normally be considered). This is sufficient to ensure the
2113 * eventual simplification of most Mapping sequences by repeated
2114 * application of this function.
2115 *
2116 * In some cases, the function may attempt more elaborate
2117 * simplification, involving any number of other Mappings in the
2118 * sequence. It is not restricted in the type or scope of
2119 * simplification it may perform, but will normally only attempt
2120 * elaborate simplification in cases where a more straightforward
2121 * approach is not adequate.
2122
2123 * Parameters:
2124 * this
2125 * Pointer to the nominated SpecMap which is to be merged with
2126 * its neighbours. This should be a cloned copy of the SpecMap
2127 * pointer contained in the array element "(*map_list)[where]"
2128 * (see below). This pointer will not be annulled, and the
2129 * SpecMap it identifies will not be modified by this function.
2130 * where
2131 * Index in the "*map_list" array (below) at which the pointer
2132 * to the nominated SpecMap resides.
2133 * series
2134 * A non-zero value indicates that the sequence of Mappings to
2135 * be simplified will be applied in series (i.e. one after the
2136 * other), whereas a zero value indicates that they will be
2137 * applied in parallel (i.e. on successive sub-sets of the
2138 * input/output coordinates).
2139 * nmap
2140 * Address of an int which counts the number of Mappings in the
2141 * sequence. On entry this should be set to the initial number
2142 * of Mappings. On exit it will be updated to record the number
2143 * of Mappings remaining after simplification.
2144 * map_list
2145 * Address of a pointer to a dynamically allocated array of
2146 * Mapping pointers (produced, for example, by the astMapList
2147 * method) which identifies the sequence of Mappings. On entry,
2148 * the initial sequence of Mappings to be simplified should be
2149 * supplied.
2150 *
2151 * On exit, the contents of this array will be modified to
2152 * reflect any simplification carried out. Any form of
2153 * simplification may be performed. This may involve any of: (a)
2154 * removing Mappings by annulling any of the pointers supplied,
2155 * (b) replacing them with pointers to new Mappings, (c)
2156 * inserting additional Mappings and (d) changing their order.
2157 *
2158 * The intention is to reduce the number of Mappings in the
2159 * sequence, if possible, and any reduction will be reflected in
2160 * the value of "*nmap" returned. However, simplifications which
2161 * do not reduce the length of the sequence (but improve its
2162 * execution time, for example) may also be performed, and the
2163 * sequence might conceivably increase in length (but normally
2164 * only in order to split up a Mapping into pieces that can be
2165 * more easily merged with their neighbours on subsequent
2166 * invocations of this function).
2167 *
2168 * If Mappings are removed from the sequence, any gaps that
2169 * remain will be closed up, by moving subsequent Mapping
2170 * pointers along in the array, so that vacated elements occur
2171 * at the end. If the sequence increases in length, the array
2172 * will be extended (and its pointer updated) if necessary to
2173 * accommodate any new elements.
2174 *
2175 * Note that any (or all) of the Mapping pointers supplied in
2176 * this array may be annulled by this function, but the Mappings
2177 * to which they refer are not modified in any way (although
2178 * they may, of course, be deleted if the annulled pointer is
2179 * the final one).
2180 * invert_list
2181 * Address of a pointer to a dynamically allocated array which,
2182 * on entry, should contain values to be assigned to the Invert
2183 * attributes of the Mappings identified in the "*map_list"
2184 * array before they are applied (this array might have been
2185 * produced, for example, by the astMapList method). These
2186 * values will be used by this function instead of the actual
2187 * Invert attributes of the Mappings supplied, which are
2188 * ignored.
2189 *
2190 * On exit, the contents of this array will be updated to
2191 * correspond with the possibly modified contents of the
2192 * "*map_list" array. If the Mapping sequence increases in
2193 * length, the "*invert_list" array will be extended (and its
2194 * pointer updated) if necessary to accommodate any new
2195 * elements.
2196 * status
2197 * Pointer to the inherited status variable.
2198
2199 * Returned Value:
2200 * If simplification was possible, the function returns the index
2201 * in the "map_list" array of the first element which was
2202 * modified. Otherwise, it returns -1 (and makes no changes to the
2203 * arrays supplied).
2204
2205 * Notes:
2206 * - A value of -1 will be returned if this function is invoked
2207 * with the global error status set, or if it should fail for any
2208 * reason.
2209 */
2210
2211 /* Local Variables: */
2212 AstMapping *new; /* Pointer to replacement Mapping */
2213 AstSpecMap *specmap; /* Pointer to SpecMap */
2214 const char *argdesc[ MAX_ARGS ]; /* Argument descriptions (junk) */
2215 const char *class; /* Pointer to Mapping class string */
2216 const char *comment; /* Pointer to comment string (junk) */
2217 double (*cvtargs)[ MAX_ARGS ]; /* Pointer to argument arrays */
2218 double tmp; /* Temporary storage */
2219 int *cvttype; /* Pointer to transformation type codes */
2220 int *szarg; /* Pointer to argument count array */
2221 int argdec; /* Index of DEC argument */
2222 int argra; /* Index of RA argument */
2223 int done; /* Finished (no further simplification)? */
2224 int iarg; /* Loop counter for arguments */
2225 int icvt1; /* Loop initial value */
2226 int icvt2; /* Loop final value */
2227 int icvt; /* Loop counter for transformation steps */
2228 int ikeep; /* Index to store step being kept */
2229 int imap1; /* Index of first SpecMap to merge */
2230 int imap2; /* Index of last SpecMap to merge */
2231 int imap; /* Loop counter for Mappings */
2232 int inc; /* Increment for transformation step loop */
2233 int invert; /* SpecMap applied in inverse direction? */
2234 int istep; /* Loop counter for transformation steps */
2235 int keep; /* Keep transformation step? */
2236 int narg; /* Number of user-supplied arguments */
2237 int ngone; /* Number of Mappings eliminated */
2238 int nin; /* Numbr of axes for SpecMaps being merged */
2239 int nstep0; /* Original number of transformation steps */
2240 int nstep; /* Total number of transformation steps */
2241 int result; /* Result value to return */
2242 int simpler; /* Simplification possible? */
2243 int unit; /* Replacement Mapping is a UnitMap? */
2244
2245 /* Initialise. */
2246 result = -1;
2247
2248 /* Check the global error status. */
2249 if ( !astOK ) return result;
2250
2251 /* SpecMaps can only be merged if they are in series (or if there is
2252 only one Mapping present, in which case it makes no difference), so
2253 do nothing if they are not. */
2254 if ( series || ( *nmap == 1 ) ) {
2255
2256 /* Save the number of inputs for the SpecMap. */
2257 nin = astGetNin( this );
2258
2259 /* Initialise the number of transformation steps to be merged to equal
2260 the number in the nominated SpecMap. */
2261 nstep = ( (AstSpecMap *) ( *map_list )[ where ] )->ncvt;
2262
2263 /* Search adjacent lower-numbered Mappings until one is found which is
2264 not a SpecMap, or is a SpecMap with a different number of axes. Accumulate
2265 the number of transformation steps involved in any SpecMaps found. */
2266 imap1 = where;
2267 while ( ( imap1 - 1 >= 0 ) && astOK ) {
2268 class = astGetClass( ( *map_list )[ imap1 - 1 ] );
2269 if ( !astOK || strcmp( class, "SpecMap" ) ||
2270 astGetNin( ( *map_list )[ imap1 - 1 ] ) != nin ) break;
2271 nstep += ( (AstSpecMap *) ( *map_list )[ imap1 - 1 ] )->ncvt;
2272 imap1--;
2273 }
2274
2275 /* Similarly search adjacent higher-numbered Mappings. */
2276 imap2 = where;
2277 while ( ( imap2 + 1 < *nmap ) && astOK ) {
2278 class = astGetClass( ( *map_list )[ imap2 + 1 ] );
2279 if ( !astOK || strcmp( class, "SpecMap" ) ||
2280 astGetNin( ( *map_list )[ imap2 + 1 ] ) != nin ) break;
2281 nstep += ( (AstSpecMap *) ( *map_list )[ imap2 + 1 ] )->ncvt;
2282 imap2++;
2283 }
2284
2285 /* Remember the initial number of transformation steps. */
2286 nstep0 = nstep;
2287
2288 /* Allocate memory for accumulating a list of all the transformation
2289 steps involved in all the SpecMaps found. */
2290 cvttype = astMalloc( sizeof( int ) * (size_t) nstep );
2291 cvtargs = astMalloc( sizeof( double[ MAX_ARGS ] ) * (size_t) nstep );
2292 szarg = astMalloc( sizeof( int ) * (size_t) nstep );
2293
2294 /* Loop to obtain the transformation data for each SpecMap being merged. */
2295 nstep = 0;
2296 for ( imap = imap1; astOK && ( imap <= imap2 ); imap++ ) {
2297
2298 /* Obtain a pointer to the SpecMap and note if it is being applied in
2299 its inverse direction. */
2300 specmap = (AstSpecMap *) ( *map_list )[ imap ];
2301 invert = ( *invert_list )[ imap ];
2302
2303 /* Set up loop limits and an increment to scan the transformation
2304 steps in each SpecMap in either the forward or reverse direction, as
2305 dictated by the associated "invert" value. */
2306 icvt1 = invert ? specmap->ncvt - 1 : 0;
2307 icvt2 = invert ? -1 : specmap->ncvt;
2308 inc = invert ? -1 : 1;
2309
2310 /* Loop through each transformation step in the SpecMap. */
2311 for ( icvt = icvt1; icvt != icvt2; icvt += inc ) {
2312
2313 /* Store the transformation type code and use "CvtString" to determine
2314 the associated number of arguments. Then store these arguments. */
2315 cvttype[ nstep ] = specmap->cvttype[ icvt ];
2316 (void) CvtString( cvttype[ nstep ], &comment, &argra, &argdec,
2317 &narg, szarg + nstep, argdesc, status );
2318 if ( !astOK ) break;
2319 for ( iarg = 0; iarg < szarg[ nstep ]; iarg++ ) {
2320 cvtargs[ nstep ][ iarg ] = specmap->cvtargs[ icvt ][ iarg ];
2321 }
2322
2323 /* If the SpecMap is inverted, we must not only accumulate its
2324 transformation steps in reverse, but also apply them in
2325 reverse. For some steps this means changing arguments, for some it
2326 means changing the transformation type code to a complementary
2327 value, and for others it means both. Define macros to perform each
2328 of the required changes. */
2329
2330 /* Macro to exchange a transformation type code for its inverse (and
2331 vice versa). */
2332 #define SWAP_CODES( code1, code2 ) \
2333 if ( cvttype[ nstep ] == code1 ) { \
2334 cvttype[ nstep ] = code2; \
2335 } else if ( cvttype[ nstep ] == code2 ) { \
2336 cvttype[ nstep ] = code1; \
2337 }
2338
2339 /* Macro to exchange a transformation type code for its inverse (and
2340 vice versa), and reciprocate a specified argument. */
2341 #define SWAP_CODES2( code1, code2, jarg ) \
2342 if ( cvttype[ nstep ] == code1 ) { \
2343 cvttype[ nstep ] = code2; \
2344 tmp = cvtargs[ nstep ][ jarg ]; \
2345 if( tmp != AST__BAD && tmp != 0.0 ) { \
2346 cvtargs[ nstep ][ jarg ] = 1.0/tmp; \
2347 } else { \
2348 cvtargs[ nstep ][ jarg ] = AST__BAD; \
2349 } \
2350 } else if ( cvttype[ nstep ] == code2 ) { \
2351 cvttype[ nstep ] = code1; \
2352 tmp = cvtargs[ nstep ][ jarg ]; \
2353 if( tmp != AST__BAD && tmp != 0.0 ) { \
2354 cvtargs[ nstep ][ jarg ] = 1.0/tmp; \
2355 } else { \
2356 cvtargs[ nstep ][ jarg ] = AST__BAD; \
2357 } \
2358 }
2359
2360 /* Macro to exchange a transformation type code for its inverse (and
2361 vice versa), and negate a specified argument. */
2362 #define SWAP_CODES3( code1, code2, jarg ) \
2363 if ( cvttype[ nstep ] == code1 ) { \
2364 cvttype[ nstep ] = code2; \
2365 tmp = cvtargs[ nstep ][ jarg ]; \
2366 if( tmp != AST__BAD ) { \
2367 cvtargs[ nstep ][ jarg ] = -tmp; \
2368 } \
2369 } else if ( cvttype[ nstep ] == code2 ) { \
2370 cvttype[ nstep ] = code1; \
2371 tmp = cvtargs[ nstep ][ jarg ]; \
2372 if( tmp != AST__BAD ) { \
2373 cvtargs[ nstep ][ jarg ] = -tmp; \
2374 } \
2375 }
2376
2377 /* Use these macros to apply the changes where needed. */
2378 if ( invert ) {
2379
2380 /* Exchange transformation codes for their inverses. */
2381 SWAP_CODES( AST__FRTOVL, AST__VLTOFR )
2382 SWAP_CODES( AST__ENTOFR, AST__FRTOEN )
2383 SWAP_CODES( AST__WNTOFR, AST__FRTOWN )
2384 SWAP_CODES( AST__WVTOFR, AST__FRTOWV )
2385 SWAP_CODES( AST__AWTOFR, AST__FRTOAW )
2386 SWAP_CODES( AST__VRTOVL, AST__VLTOVR )
2387 SWAP_CODES( AST__VOTOVL, AST__VLTOVO )
2388 SWAP_CODES( AST__ZOTOVL, AST__VLTOZO )
2389 SWAP_CODES( AST__BTTOVL, AST__VLTOBT )
2390
2391 /* Exchange transformation codes for their inverses, and reciprocate the
2392 frequency correction factor. */
2393 SWAP_CODES2( AST__TPF2HL, AST__HLF2TP, 6 )
2394 SWAP_CODES2( AST__USF2HL, AST__HLF2US, 3 )
2395 SWAP_CODES2( AST__GEF2HL, AST__HLF2GE, 3 )
2396 SWAP_CODES2( AST__BYF2HL, AST__HLF2BY, 3 )
2397 SWAP_CODES2( AST__LKF2HL, AST__HLF2LK, 2 )
2398 SWAP_CODES2( AST__LDF2HL, AST__HLF2LD, 2 )
2399 SWAP_CODES2( AST__LGF2HL, AST__HLF2LG, 2 )
2400 SWAP_CODES2( AST__GLF2HL, AST__HLF2GL, 2 )
2401
2402 }
2403
2404 /* Undefine the local macros. */
2405 #undef SWAP_CODES
2406 #undef SWAP_CODES2
2407 #undef SWAP_CODES3
2408
2409 /* Count the transformation steps. */
2410 nstep++;
2411 }
2412 }
2413
2414 /* Loop to simplify the sequence of transformation steps until no
2415 further improvement is possible. */
2416 done = 0;
2417 while ( astOK && !done ) {
2418
2419 /* Examine each remaining transformation step in turn. */
2420 ikeep = -1;
2421 for ( istep = 0; istep < nstep; istep++ ) {
2422
2423 /* Initially assume we will retain the current step. */
2424 keep = 1;
2425
2426 /* The only simplifications for the conversions currently in this class act
2427 to combine adjacent transformation steps, so only apply them while there
2428 are at least 2 steps left. */
2429 if ( istep < ( nstep - 1 ) ) {
2430
2431 /* Define a macro to test if two adjacent transformation type codes
2432 have specified values. */
2433 #define PAIR_CVT( code1, code2 ) \
2434 ( ( cvttype[ istep ] == code1 ) && \
2435 ( cvttype[ istep + 1 ] == code2 ) )
2436
2437 /* Define a macro to test if two adjacent transformation type codes
2438 have specified values, either way round. */
2439 #define PAIR_CVT2( code1, code2 ) \
2440 ( ( PAIR_CVT( code1, code2 ) ) || \
2441 ( PAIR_CVT( code2, code1 ) ) )
2442
2443 /* If a correction is followed by its inverse, and the user-supplied argument
2444 values are unchanged (we do not need to test values stored in the
2445 argument array which were not supplied by the user), we can eliminate them.
2446 First check for conversions which have no user-supplied arguments. */
2447 if ( PAIR_CVT2( AST__ENTOFR, AST__FRTOEN ) ||
2448 PAIR_CVT2( AST__WNTOFR, AST__FRTOWN ) ||
2449 PAIR_CVT2( AST__WVTOFR, AST__FRTOWV ) ||
2450 PAIR_CVT2( AST__AWTOFR, AST__FRTOAW ) ||
2451 PAIR_CVT2( AST__VRTOVL, AST__VLTOVR ) ||
2452 PAIR_CVT2( AST__VOTOVL, AST__VLTOVO ) ||
2453 PAIR_CVT2( AST__ZOTOVL, AST__VLTOZO ) ||
2454 PAIR_CVT2( AST__BTTOVL, AST__VLTOBT ) ) {
2455 istep++;
2456 keep = 0;
2457
2458 /* Now check for conversions which have a single user-supplied argument. */
2459 } else if( PAIR_CVT2( AST__FRTOVL, AST__VLTOFR ) &&
2460 EQUAL( cvtargs[ istep ][ 0 ],
2461 cvtargs[ istep + 1 ][ 0 ] ) ) {
2462 istep++;
2463 keep = 0;
2464
2465 /* Now check for conversions which have two user-supplied arguments. */
2466 } else if( ( PAIR_CVT2( AST__LKF2HL, AST__HLF2LK ) ||
2467 PAIR_CVT2( AST__LDF2HL, AST__HLF2LD ) ||
2468 PAIR_CVT2( AST__LGF2HL, AST__HLF2LG ) ||
2469 PAIR_CVT2( AST__GLF2HL, AST__HLF2GL ) ) &&
2470 EQUAL( cvtargs[ istep ][ 0 ],
2471 cvtargs[ istep + 1 ][ 0 ] ) &&
2472 EQUAL( cvtargs[ istep ][ 1 ],
2473 cvtargs[ istep + 1 ][ 1 ] ) ) {
2474 istep++;
2475 keep = 0;
2476
2477 /* Now check for conversions which have three user-supplied arguments. */
2478 } else if( ( PAIR_CVT2( AST__GEF2HL, AST__HLF2GE ) ||
2479 PAIR_CVT2( AST__BYF2HL, AST__HLF2BY ) ||
2480 PAIR_CVT2( AST__USF2HL, AST__HLF2US ) ) &&
2481 EQUAL( cvtargs[ istep ][ 0 ],
2482 cvtargs[ istep + 1 ][ 0 ] ) &&
2483 EQUAL( cvtargs[ istep ][ 1 ],
2484 cvtargs[ istep + 1 ][ 1 ] ) &&
2485 EQUAL( cvtargs[ istep ][ 2 ],
2486 cvtargs[ istep + 1 ][ 2 ] ) ) {
2487 istep++;
2488 keep = 0;
2489
2490 /* Now check for conversions which have six user-supplied arguments (currently
2491 no conversions have four or five user-supplied arguments). */
2492 } else if( ( PAIR_CVT2( AST__TPF2HL, AST__HLF2TP ) ) &&
2493 EQUAL( cvtargs[ istep ][ 0 ],
2494 cvtargs[ istep + 1 ][ 0 ] ) &&
2495 EQUAL( cvtargs[ istep ][ 1 ],
2496 cvtargs[ istep + 1 ][ 1 ] ) &&
2497 EQUAL( cvtargs[ istep ][ 2 ],
2498 cvtargs[ istep + 1 ][ 2 ] ) &&
2499 EQUAL( cvtargs[ istep ][ 3 ],
2500 cvtargs[ istep + 1 ][ 3 ] ) &&
2501 EQUAL( cvtargs[ istep ][ 4 ],
2502 cvtargs[ istep + 1 ][ 4 ] ) &&
2503 EQUAL( cvtargs[ istep ][ 5 ],
2504 cvtargs[ istep + 1 ][ 5 ] ) ) {
2505 istep++;
2506 keep = 0;
2507
2508 }
2509
2510 /* Undefine the local macros. */
2511 #undef PAIR_CVT
2512 #undef PAIR_CVT2
2513 }
2514
2515 /* If the current transformation (possibly modified above) is being
2516 kept, then increment the index that identifies its new location in
2517 the list of transformation steps. */
2518 if ( keep ) {
2519 ikeep++;
2520
2521 /* If the new location is different to its current location, copy the
2522 transformation data into the new location. */
2523 if ( ikeep != istep ) {
2524 cvttype[ ikeep ] = cvttype[ istep ];
2525 for ( iarg = 0; iarg < szarg[ istep ]; iarg++ ) {
2526 cvtargs[ ikeep ][ iarg ] = cvtargs[ istep ][ iarg ];
2527 }
2528 szarg[ ikeep ] = szarg[ istep ];
2529 }
2530 }
2531 }
2532
2533 /* Note if no simplification was achieved on this iteration (i.e. the
2534 number of transformation steps was not reduced). This is the signal
2535 to quit. */
2536 done = ( ( ikeep + 1 ) >= nstep );
2537
2538 /* Note how many transformation steps now remain. */
2539 nstep = ikeep + 1;
2540 }
2541
2542 /* Determine how many Mappings can be eliminated by condensing all
2543 those considered above into a single Mapping. */
2544 if ( astOK ) {
2545 ngone = imap2 - imap1;
2546
2547 /* Determine if the replacement Mapping can be a UnitMap (a null
2548 Mapping). This will only be the case if all the transformation
2549 steps were eliminated above. */
2550 unit = ( nstep == 0 );
2551
2552 /* Determine if simplification is possible. This will be the case if
2553 (a) Mappings were eliminated ("ngone" is non-zero), or (b) the
2554 number of transformation steps was reduced, or (c) the SpecMap(s)
2555 can be replaced by a UnitMap, or (d) if there was initially only
2556 one SpecMap present, its invert flag was set (this flag will always
2557 be cleared in the replacement Mapping). */
2558 simpler = ngone || ( nstep < nstep0 ) || unit ||
2559 ( *invert_list )[ where ];
2560
2561 /* Do nothing more unless simplification is possible. */
2562 if ( simpler ) {
2563
2564 /* If the replacement Mapping is a UnitMap, then create it. */
2565 if ( unit ) {
2566 new = (AstMapping *)
2567 astUnitMap( astGetNin( ( *map_list )[ where ] ), "", status );
2568
2569 /* Otherwise, create a replacement SpecMap and add each of the
2570 remaining transformation steps to it. */
2571 } else {
2572 new = (AstMapping *) astSpecMap( nin, 0, "", status );
2573 for ( istep = 0; istep < nstep; istep++ ) {
2574 AddSpecCvt( (AstSpecMap *) new, cvttype[ istep ],
2575 cvtargs[ istep ], status );
2576 }
2577 }
2578
2579 /* Annul the pointers to the Mappings being eliminated. */
2580 if ( astOK ) {
2581 for ( imap = imap1; imap <= imap2; imap++ ) {
2582 ( *map_list )[ imap ] = astAnnul( ( *map_list )[ imap ] );
2583 }
2584
2585 /* Insert the pointer and invert value for the new Mapping. */
2586 ( *map_list )[ imap1 ] = new;
2587 ( *invert_list )[ imap1 ] = 0;
2588
2589 /* Move any subsequent Mapping information down to close the gap. */
2590 for ( imap = imap2 + 1; imap < *nmap; imap++ ) {
2591 ( *map_list )[ imap - ngone ] = ( *map_list )[ imap ];
2592 ( *invert_list )[ imap - ngone ] = ( *invert_list )[ imap ];
2593 }
2594
2595 /* Blank out any information remaining at the end of the arrays. */
2596 for ( imap = ( *nmap - ngone ); imap < *nmap; imap++ ) {
2597 ( *map_list )[ imap ] = NULL;
2598 ( *invert_list )[ imap ] = 0;
2599 }
2600
2601 /* Decrement the Mapping count and return the index of the first
2602 Mapping which was eliminated. */
2603 ( *nmap ) -= ngone;
2604 result = imap1;
2605
2606 /* If an error occurred, annul the new Mapping pointer. */
2607 } else {
2608 new = astAnnul( new );
2609 }
2610 }
2611 }
2612
2613 /* Free the memory used for the transformation steps. */
2614 cvttype = astFree( cvttype );
2615 cvtargs = astFree( cvtargs );
2616 szarg = astFree( szarg );
2617 }
2618
2619 /* If an error occurred, clear the returned value. */
2620 if ( !astOK ) result = -1;
2621
2622 /* Return the result. */
2623 return result;
2624 }
2625
Rate(AstMapping * this,double * at,int ax1,int ax2,int * status)2626 static double Rate( AstMapping *this, double *at, int ax1, int ax2, int *status ){
2627 /*
2628 * Name:
2629 * Rate
2630
2631 * Purpose:
2632 * Calculate the rate of change of a Mapping output.
2633
2634 * Type:
2635 * Private function.
2636
2637 * Synopsis:
2638 * #include "specmap.h"
2639 * result = Rate( AstMapping *this, double *at, int ax1, int ax2, int *status )
2640
2641 * Class Membership:
2642 * SpecMap member function (overrides the astRate method inherited
2643 * from the Mapping class ).
2644
2645 * Description:
2646 * This function returns the rate of change of a specified output of
2647 * the supplied Mapping with respect to a specified input, at a
2648 * specified input position.
2649
2650 * Parameters:
2651 * this
2652 * Pointer to the Mapping to be applied.
2653 * at
2654 * The address of an array holding the axis values at the position
2655 * at which the rate of change is to be evaluated. The number of
2656 * elements in this array should equal the number of inputs to the
2657 * Mapping.
2658 * ax1
2659 * The index of the Mapping output for which the rate of change is to
2660 * be found (output numbering starts at 0 for the first output).
2661 * ax2
2662 * The index of the Mapping input which is to be varied in order to
2663 * find the rate of change (input numbering starts at 0 for the first
2664 * input).
2665 * status
2666 * Pointer to the inherited status variable.
2667
2668 * Returned Value:
2669 * The rate of change of Mapping output "ax1" with respect to input
2670 * "ax2", evaluated at "at", or AST__BAD if the value cannot be
2671 * calculated.
2672
2673 * Implementation Deficiencies:
2674 * The initial version of this implementation only deals with
2675 * frequency->wavelength conversions. This is because the slowness of
2676 * the numerical differentiation implemented by the astRate method in
2677 * the parent Mapping class is cripples conversion between SpecFluxFrames.
2678 * Such conversions only rely on rate of change of wavelength with
2679 * respect to frequency. This implementation should be extended when
2680 * needed.
2681
2682 */
2683
2684 /* Local Variables: */
2685 AstSpecMap *map;
2686 double result;
2687 int cvt;
2688
2689 /* Check inherited status */
2690 if( !astOK ) return AST__BAD;
2691
2692 /* Get a pointer to the SpecMap structure. */
2693 map = (AstSpecMap *) this;
2694
2695 /* Return 1.0 if the SpecMap has no conversions. */
2696 if( map->ncvt == 0 ) return 1.0;
2697
2698 /* Store the type of the first conversion.*/
2699 cvt = map->cvttype[ 0 ];
2700
2701 /* If this is a 3D SpecMap or if it has more than one component, or if
2702 that conversion is not between frequency and wavelength, use the
2703 astRate method inherited form the parent Mapping class. */
2704 if( astGetNin( map ) != 1 || map->ncvt != 1 ||
2705 ( cvt != AST__WVTOFR && cvt != AST__FRTOWV ) ) {
2706 result = (*parent_rate)( this, at, ax1, ax2, status );
2707
2708 /* Otherwise, evaluate the known analytical expressions for the rate of
2709 change of frequency with respect to wavelength or wavelength with
2710 respect to frequency. */
2711 } else {
2712 result = ( *at != AST__BAD ) ? -AST__C/((*at)*(*at)) : AST__BAD;
2713 }
2714
2715 /* Return the result. */
2716 return result;
2717 }
2718
Refrac(double wavelen,int * status)2719 static double Refrac( double wavelen, int *status ){
2720 /*
2721 * Name:
2722 * Refrac
2723
2724 * Purpose:
2725 * Returns the refractive index of dry air at a given wavelength.
2726
2727 * Type:
2728 * Private function.
2729
2730 * Synopsis:
2731 * #include "specmap.h"
2732 * double Refrac( double wavelen, int *status )
2733
2734 * Class Membership:
2735 * SpecMap method.
2736
2737 * Description:
2738 * This function returns the refractive index of dry air at standard
2739 * temperature and pressure, at a given wavelength. The formula is
2740 * taken from the paper "Representation of Spectral Coordinates in FITS"
2741 * (Greisen et al).
2742
2743 * Parameters:
2744 * wavelen
2745 * The wavelength, in metres. This should be the air wavelength,
2746 * but supplying the vacuum wavelength will make no significant
2747 * difference.
2748 * status
2749 * Pointer to the inherited status variable.
2750
2751 * Returned Value:
2752 * The refractive index. A value of 1.0 is returned if an error
2753 * occurs, or has already occurred.
2754
2755 */
2756
2757 /* Local Variables: */
2758 double w2; /* Wavenumber squared */
2759
2760 /* Check the global error status. */
2761 if ( !astOK || wavelen == 0.0 ) return 1.0;
2762
2763 /* Find the squared wave number in units of "(per um)**2". */
2764 w2 = 1.0E-12/( wavelen * wavelen );
2765
2766 /* Apply the rest of the algorithm as described in the FITS WCS
2767 paper III. */
2768 return 1.0 + 1.0E-6*( 287.6155 + 1.62887*w2 + 0.01360*w2*w2 );
2769 }
2770
Rverot(double phi,double h,double ra,double da,double st,int * status)2771 static double Rverot( double phi, double h, double ra, double da,
2772 double st, int *status ) {
2773 /*
2774 * Name:
2775 * Rverot
2776
2777 * Purpose:
2778 * Find the velocity component in a given direction due to Earth rotation.
2779
2780 * Type:
2781 * Private function.
2782
2783 * Synopsis:
2784 * #include "specmap.h"
2785 * double Rverot( double phi, double h, double ra, double da,
2786 * double st, int *status )
2787
2788 * Class Membership:
2789 * SpecMap method.
2790
2791 * Description:
2792 * This function is like slaRverot, except that it takes account of the
2793 * observers height (h), and does all calculations in double precision.
2794
2795 * Parameters:
2796 * phi
2797 * The geodetic latitude of the observer (radians, IAU 1976).
2798 * h
2799 * The geodetic height above the reference spheroid of the observer
2800 * (metres, IAU 1976).
2801 * ra
2802 * The geocentric apparent RA (rads) of the source.
2803 * da
2804 * The geocentric apparent Dec (rads) of the source.
2805 * st
2806 * The local apparent sidereal time (radians).
2807 * status
2808 * Pointer to the inherited status variable.
2809
2810 * Returns:
2811 * The component of the Earth rotation in direction [RA,DA] (km/s).
2812 * The result is positive when the observer is receding from the
2813 * given point on the sky. Zero is returned if an error has already
2814 * occurred.
2815
2816 */
2817
2818 /* Local Variables: */
2819 double pv[ 6 ]; /* Observer position and velocity */
2820 double v[ 3 ]; /* Source direction vector */
2821
2822 /* Check the global error status. */
2823 if ( !astOK ) return 0.0;
2824
2825 /* Get the Cartesian coordinates of the unit vector pointing towards the
2826 given sky position. */
2827 palDcs2c( ra, da, v );
2828
2829 /* Get velocity and position of the observer. */
2830 palPvobs( phi, h, st, pv );
2831
2832 /* Return the component of the observer's velocity away from the sky
2833 position, and convert from AU/s to km/s. */
2834 return -palDvdv( v, pv + 3 )*149.597870E6;
2835 }
2836
SpecAdd(AstSpecMap * this,const char * cvt,const double args[],int * status)2837 static void SpecAdd( AstSpecMap *this, const char *cvt, const double args[], int *status ) {
2838 /*
2839 *++
2840 * Name:
2841 c astSpecAdd
2842 f AST_SPECADD
2843
2844 * Purpose:
2845 * Add a spectral coordinate conversion to a SpecMap.
2846
2847 * Type:
2848 * Public virtual function.
2849
2850 * Synopsis:
2851 c #include "specmap.h"
2852 c void astSpecAdd( AstSpecMap *this, const char *cvt, const double args[] )
2853 f CALL AST_SPECADD( THIS, CVT, ARGS, STATUS )
2854
2855 * Class Membership:
2856 * SpecMap method.
2857
2858 * Description:
2859 c This function adds one of the standard spectral coordinate
2860 f This routine adds one of the standard spectral coordinate
2861 * system conversions listed below to an existing SpecMap.
2862 *
2863 c When a SpecMap is first created (using astSpecMap), it simply
2864 f When a SpecMap is first created (using AST_SPECMAP), it simply
2865 c performs a unit (null) Mapping. By using astSpecAdd (repeatedly
2866 f performs a unit (null) Mapping. By using AST_SPECADD (repeatedly
2867 * if necessary), one or more coordinate conversion steps may then
2868 * be added, which the SpecMap will perform in sequence. This allows
2869 * multi-step conversions between a variety of spectral coordinate
2870 * systems to be assembled out of the building blocks provided by
2871 * this class.
2872 *
2873 * Normally, if a SpecMap's Invert attribute is zero (the default),
2874 * then its forward transformation is performed by carrying out
2875 * each of the individual coordinate conversions specified by
2876 c astSpecAdd in the order given (i.e. with the most recently added
2877 f AST_SPECADD in the order given (i.e. with the most recently added
2878 * conversion applied last).
2879 *
2880 * This order is reversed if the SpecMap's Invert attribute is
2881 * non-zero (or if the inverse transformation is requested by any
2882 * other means) and each individual coordinate conversion is also
2883 * replaced by its own inverse. This process inverts the overall
2884 * effect of the SpecMap. In this case, the first conversion to be
2885 * applied would be the inverse of the one most recently added.
2886
2887 * Parameters:
2888 c this
2889 f THIS = INTEGER (Given)
2890 * Pointer to the SpecMap.
2891 c cvt
2892 f CVT = CHARACTER * ( * ) (Given)
2893 c Pointer to a null-terminated string which identifies the
2894 f A character string which identifies the
2895 * spectral coordinate conversion to be added to the
2896 * SpecMap. See the "Available Conversions" section for details of
2897 * those available.
2898 c args
2899 f ARGS( * ) = DOUBLE PRECISION (Given)
2900 * An array containing argument values for the spectral
2901 * coordinate conversion. The number of arguments required, and
2902 * hence the number of array elements used, depends on the
2903 * conversion specified (see the "Available Conversions"
2904 * section). This array is ignored
2905 c and a NULL pointer may be supplied
2906 * if no arguments are needed.
2907 f STATUS = INTEGER (Given and Returned)
2908 f The global status.
2909
2910 * Notes:
2911 * - When assembling a multi-stage conversion, it can sometimes be
2912 * difficult to determine the most economical conversion path. For
2913 * example, when converting between reference frames, converting first
2914 * to the heliographic reference frame as an intermediate stage is often
2915 * sensible in formulating the problem, but may introduce unnecessary
2916 * extra conversion steps. A solution to this is to include all the steps
2917 * which are (logically) necessary, but then to use
2918 c astSimplify to simplify the resulting
2919 f AST_SIMPLIFY to simplify the resulting
2920 * SpecMap. The simplification process will eliminate any steps
2921 * which turn out not to be needed.
2922 c - This function does not check to ensure that the sequence of
2923 f - This routine does not check to ensure that the sequence of
2924 * coordinate conversions added to a SpecMap is physically
2925 * meaningful.
2926
2927 * Available Conversions:
2928 * The following strings (which are case-insensitive) may be supplied
2929 c via the "cvt" parameter to indicate which spectral coordinate
2930 f via the CVT argument to indicate which spectral coordinate
2931 * conversion is to be added to the SpecMap. Where arguments are needed by
2932 * the conversion, they are listed in parentheses. Values for
2933 c these arguments should be given, via the "args" array, in the
2934 f these arguments should be given, via the ARGS array, in the
2935 * order indicated. Units and argument names are described at the end of
2936 * the list of conversions.
2937
2938 * - "FRTOVL" (RF): Convert frequency to relativistic velocity.
2939 * - "VLTOFR" (RF): Convert relativistic velocity to Frequency.
2940 * - "ENTOFR": Convert energy to frequency.
2941 * - "FRTOEN": Convert frequency to energy.
2942 * - "WNTOFR": Convert wave number to frequency.
2943 * - "FRTOWN": Convert frequency to wave number.
2944 * - "WVTOFR": Convert wavelength (vacuum) to frequency.
2945 * - "FRTOWV": Convert frequency to wavelength (vacuum).
2946 * - "AWTOFR": Convert wavelength (air) to frequency.
2947 * - "FRTOAW": Convert frequency to wavelength (air).
2948 * - "VRTOVL": Convert radio to relativistic velocity.
2949 * - "VLTOVR": Convert relativistic to radio velocity.
2950 * - "VOTOVL": Convert optical to relativistic velocity.
2951 * - "VLTOVO": Convert relativistic to optical velocity.
2952 * - "ZOTOVL": Convert redshift to relativistic velocity.
2953 * - "VLTOZO": Convert relativistic velocity to redshift.
2954 * - "BTTOVL": Convert beta factor to relativistic velocity.
2955 * - "VLTOBT": Convert relativistic velocity to beta factor.
2956 * - "USF2HL" (VOFF,RA,DEC): Convert frequency from a user-defined
2957 * reference frame to heliocentric.
2958 * - "HLF2US" (VOFF,RA,DEC): Convert frequency from heliocentric
2959 * reference frame to user-defined.
2960 * - "TPF2HL" (OBSLON,OBSLAT,OBSALT,EPOCH,RA,DEC): Convert frequency from
2961 * topocentric reference frame to heliocentric.
2962 * - "HLF2TP" (OBSLON,OBSLAT,OBSALT,EPOCH,RA,DEC): Convert frequency from
2963 * heliocentric reference frame to topocentric.
2964 * - "GEF2HL" (EPOCH,RA,DEC): Convert frequency from geocentric
2965 * reference frame to heliocentric.
2966 * - "HLF2GE" (EPOCH,RA,DEC): Convert frequency from
2967 * heliocentric reference frame to geocentric.
2968 * - "BYF2HL" (EPOCH,RA,DEC): Convert frequency from
2969 * barycentric reference frame to heliocentric.
2970 * - "HLF2BY" (EPOCH,RA,DEC): Convert frequency from
2971 * heliocentric reference frame to barycentric.
2972 * - "LKF2HL" (RA,DEC): Convert frequency from kinematic LSR
2973 * reference frame to heliocentric.
2974 * - "HLF2LK" (RA,DEC): Convert frequency from heliocentric
2975 * reference frame to kinematic LSR.
2976 * - "LDF2HL" (RA,DEC): Convert frequency from dynamical LSR
2977 * reference frame to heliocentric.
2978 * - "HLF2LD" (RA,DEC): Convert frequency from heliocentric
2979 * reference frame to dynamical LSR.
2980 * - "LGF2HL" (RA,DEC): Convert frequency from local group
2981 * reference frame to heliocentric.
2982 * - "HLF2LG" (RA,DEC): Convert frequency from heliocentric
2983 * reference frame to local group.
2984 * - "GLF2HL" (RA,DEC): Convert frequency from galactic
2985 * reference frame to heliocentric.
2986 * - "HLF2GL" (RA,DEC): Convert frequency from heliocentric
2987 * reference frame to galactic.
2988
2989 * The units for the values processed by the above conversions are as
2990 * follows:
2991 *
2992 * - all velocities: metres per second (positive if the source receeds from
2993 * the observer).
2994 * - frequency: Hertz.
2995 * - all wavelengths: metres.
2996 * - energy: Joules.
2997 * - wave number: cycles per metre.
2998 *
2999 * The arguments used in the above conversions are as follows:
3000 *
3001 * - RF: Rest frequency (Hz).
3002 * - OBSALT: Geodetic altitude of observer (IAU 1975, metres).
3003 * - OBSLAT: Geodetic latitude of observer (IAU 1975, radians).
3004 * - OBSLON: Longitude of observer (radians - positive eastwards).
3005 * - EPOCH: Epoch of observation (UT1 expressed as a Modified Julian Date).
3006 * - RA: Right Ascension of source (radians, FK5 J2000).
3007 * - DEC: Declination of source (radians, FK5 J2000).
3008 * - VOFF: Velocity of the user-defined reference frame, towards the
3009 * position given by RA and DEC, measured in the heliocentric
3010 * reference frame.
3011 *
3012 * If the SpecMap is 3-dimensional, source positions are provided by the
3013 * values supplied to inputs 2 and 3 of the SpecMap (which are simply
3014 * copied to outputs 2 and 3). Note, usable values are still required
3015 * for the RA and DEC arguments in order to define the "user-defined"
3016 * reference frame used by USF2HL and HLF2US. However, AST__BAD can be
3017 * supplied for RA and DEC if the user-defined reference frame is not
3018 * required.
3019 *
3020 *--
3021 */
3022
3023 /* Local Variables: */
3024 int cvttype; /* Conversion type code */
3025
3026 /* Check the inherited status. */
3027 if ( !astOK ) return;
3028
3029 /* Validate the type string supplied and obtain the equivalent
3030 conversion type code. */
3031 cvttype = CvtCode( cvt, status );
3032
3033 /* If the string was not recognised, then report an error. */
3034 if ( astOK && ( cvttype == AST__SPEC_NULL ) ) {
3035 astError( AST__SPCIN,
3036 "%s(%s): Invalid SpecMap spectral coordinate "
3037 "conversion type \"%s\".", status, "astAddSpec", astGetClass( this ), cvt );
3038 }
3039
3040 /* Add the new conversion to the SpecMap. */
3041 AddSpecCvt( this, cvttype, args, status );
3042 }
3043
SystemChange(int cvt_code,int np,double * values,double * args,int forward,int * status)3044 static int SystemChange( int cvt_code, int np, double *values, double *args,
3045 int forward, int *status ){
3046 /*
3047 * Name:
3048 * SystemChange
3049
3050 * Purpose:
3051 * Change values between two spectral systems.
3052
3053 * Type:
3054 * Private function.
3055
3056 * Synopsis:
3057 * #include "specmap.h"
3058 * int SystemChange( int cvt_code, int np, double *values, double *args,
3059 * int forward, int *status )
3060
3061 * Class Membership:
3062 * SpecMap method.
3063
3064 * Description:
3065 * This function modifies the supplied values in order to change the
3066 * spectral co-ordinate system (frequency, wavelength, etc) to which
3067 * they refer.
3068
3069 * Parameters:
3070 * cvt_code
3071 * A code indicating the conversion to be applied. If the code does
3072 * not correspond to a change of system, then the supplied values
3073 * are left unchanged and zero is returned as the function value.
3074 * np
3075 * The number of frequency values to transform.
3076 * values
3077 * Pointer to an array of "np" spectral values. These are modified on
3078 * return to hold the corresponding values measured in the output
3079 * system.
3080 * args
3081 * Pointer to an array holding the conversion arguments. The number
3082 * of arguments expected depends on the particular conversion being
3083 * used.
3084 * forward
3085 * Should the conversion be applied in the forward or inverse
3086 * direction? Non-zero for forward, zero for inverse.
3087 * status
3088 * Pointer to the inherited status variable.
3089
3090 * Returned Value:
3091 * Non-zero if the supplied conversion code corresponds to a change of
3092 * system. Zero otherwise (in which case the upplied values will not
3093 * have been changed).
3094
3095 */
3096
3097 /* Local Variables: */
3098 double *pv; /* Pointer to next value */
3099 double d; /* Intermediate value */
3100 double f2; /* Squared frequency */
3101 double temp; /* Intermediate value */
3102 int i; /* Loop index */
3103 int iter; /* Iteration count */
3104 int result; /* Returned value */
3105
3106 /* Check inherited status. */
3107 if( !astOK ) return 0;
3108
3109 /* Set the return value to indicate that the supplied conversion code
3110 represents a change of system. */
3111 result = 1;
3112
3113 /* Test for each code value in turn and assign the appropriate values. */
3114 switch ( cvt_code ) {
3115
3116 /* Frequency to relativistic velocity. */
3117 case AST__FRTOVL:
3118 if( forward ) {
3119 if( args[ 0 ] != AST__BAD ) {
3120 temp = args[ 0 ] * args[ 0 ];
3121 pv = values - 1;
3122 for( i = 0; i < np; i++ ){
3123 pv++;
3124 if( *pv != AST__BAD ) {
3125 f2 = ( *pv ) * ( *pv );
3126 d = temp + f2;
3127 if( d > 0.0 ) {
3128 *pv = AST__C*( ( temp - f2 )/d );
3129 if( !astISFINITE( *pv ) ) *pv = AST__BAD;
3130 } else {
3131 *pv = AST__BAD;
3132 }
3133 }
3134 }
3135 } else {
3136 pv = values;
3137 for( i = 0; i < np; i++ ) *( pv++ ) = AST__BAD;
3138 }
3139 } else {
3140 SystemChange( AST__VLTOFR, np, values, args, 1, status );
3141 }
3142 break;
3143
3144 /* Relativistic velocity to frequency. */
3145 case AST__VLTOFR:
3146 if( forward ) {
3147 if( args[ 0 ] != AST__BAD ) {
3148 temp = args[ 0 ];
3149 pv = values - 1;
3150 for( i = 0; i < np; i++ ){
3151 pv++;
3152 if( *pv != AST__BAD ) {
3153 d = AST__C + ( *pv );
3154 if( d != 0.0 ) {
3155 d = ( AST__C - ( *pv ) )/d;
3156 if( d >= 0.0 ) {
3157 *pv = temp*sqrt( d );
3158 if( !astISFINITE( *pv ) ) *pv = AST__BAD;
3159 } else {
3160 *pv = AST__BAD;
3161 }
3162 } else {
3163 *pv = AST__BAD;
3164 }
3165 }
3166 }
3167 } else {
3168 pv = values;
3169 for( i = 0; i < np; i++ ) *( pv++ ) = AST__BAD;
3170 }
3171 } else {
3172 SystemChange( AST__FRTOVL, np, values, args, 1, status );
3173 }
3174 break;
3175
3176 /* Energy to frequency */
3177 case AST__ENTOFR:
3178 if( forward ) {
3179 pv = values - 1;
3180 for( i = 0; i < np; i++ ) {
3181 pv++;
3182 if( *pv != AST__BAD ) {
3183 *pv /= AST__H;
3184 }
3185 }
3186 } else {
3187 SystemChange( AST__FRTOEN, np, values, args, 1, status );
3188 }
3189 break;
3190
3191 /* Frequency to energy */
3192 case AST__FRTOEN:
3193 if( forward ) {
3194 pv = values - 1;
3195 for( i = 0; i < np; i++ ) {
3196 pv++;
3197 if( *pv != AST__BAD ) {
3198 *pv *= AST__H;
3199 }
3200 }
3201 } else {
3202 SystemChange( AST__ENTOFR, np, values, args, 1, status );
3203 }
3204 break;
3205
3206 /* Wave number to frequency */
3207 case AST__WNTOFR:
3208 if( forward ) {
3209 pv = values - 1;
3210 for( i = 0; i < np; i++ ) {
3211 pv++;
3212 if( *pv != AST__BAD ) {
3213 *pv *= AST__C;
3214 }
3215 }
3216 } else {
3217 SystemChange( AST__FRTOWN, np, values, args, 1, status );
3218 }
3219 break;
3220
3221 /* Wave number to frequency */
3222 case AST__FRTOWN:
3223 if( forward ) {
3224 pv = values - 1;
3225 for( i = 0; i < np; i++ ) {
3226 pv++;
3227 if( *pv != AST__BAD ) {
3228 *pv /= AST__C;
3229 }
3230 }
3231 } else {
3232 SystemChange( AST__WNTOFR, np, values, args, 1, status );
3233 }
3234 break;
3235
3236 /* Wavelength to frequency */
3237 case AST__WVTOFR:
3238 if( forward ) {
3239 pv = values - 1;
3240 for( i = 0; i < np; i++ ) {
3241 pv++;
3242 if( *pv != AST__BAD && *pv != 0.0 ) {
3243 *pv = AST__C/( *pv );
3244 if( !astISFINITE( *pv ) ) *pv = AST__BAD;
3245 } else {
3246 *pv = AST__BAD;
3247 }
3248 }
3249 } else {
3250 SystemChange( AST__FRTOWV, np, values, args, 1, status );
3251 }
3252 break;
3253
3254 /* Frequency to wavelength. */
3255 case AST__FRTOWV:
3256 if( forward ) {
3257 pv = values - 1;
3258 for( i = 0; i < np; i++ ) {
3259 pv++;
3260 if( *pv != AST__BAD && *pv != 0.0 ) {
3261 *pv = AST__C/( *pv );
3262 if( !astISFINITE( *pv ) ) *pv = AST__BAD;
3263 } else {
3264 *pv = AST__BAD;
3265 }
3266 }
3267 } else {
3268 SystemChange( AST__WVTOFR, np, values, args, 1, status );
3269 }
3270 break;
3271
3272 /* Wavelength in air to frequency. */
3273 case AST__AWTOFR:
3274 if( forward ) {
3275 pv = values - 1;
3276 for( i = 0; i < np; i++ ) {
3277 pv++;
3278 if( *pv != AST__BAD && *pv != 0.0 ) {
3279 *pv = AST__C/( ( *pv )*Refrac( *pv, status ) );
3280 if( !astISFINITE( *pv ) ) *pv = AST__BAD;
3281 } else {
3282 *pv = AST__BAD;
3283 }
3284 }
3285 } else {
3286 SystemChange( AST__FRTOAW, np, values, args, 1, status );
3287 }
3288 break;
3289
3290 /* Frequency to wavelength in air. */
3291 case AST__FRTOAW:
3292 if( forward ) {
3293 pv = values - 1;
3294 for( i = 0; i < np; i++ ) {
3295 pv++;
3296 if( *pv != AST__BAD && *pv != 0.0 ) {
3297
3298 /* Form the vacuum wavelength. */
3299 temp = AST__C/( *pv );
3300
3301 /* The refractive index function "Refrac" requires the wavelength in air
3302 as its parameter. Initially assume that the wavelength in air is equal
3303 to the vacuum wavelength to get he first estimate of the wavelength in
3304 air. Then use this estimate to get a better refractive index in order to
3305 form a better estimate of the air wavelength, etc. Iterate in this way a
3306 few times. */
3307 *pv = temp;
3308 for( iter = 0; iter < 3; iter++ ) {
3309 *pv = temp/Refrac( *pv, status );
3310 if( !astISFINITE( *pv ) ) {
3311 *pv = AST__BAD;
3312 break;
3313 }
3314 }
3315
3316 } else {
3317 *pv = AST__BAD;
3318 }
3319 }
3320 } else {
3321 SystemChange( AST__AWTOFR, np, values, args, 1, status );
3322 }
3323 break;
3324
3325 /* Radio velocity to relativistic velocity */
3326 case AST__VRTOVL:
3327 if( forward ) {
3328 pv = values - 1;
3329 for( i = 0; i < np; i++ ) {
3330 pv++;
3331 if( *pv != AST__BAD ) {
3332 temp = 1.0 - ( *pv )/AST__C;
3333 temp *= temp;
3334 *pv = AST__C*( 1.0 - temp )/( 1.0 + temp );
3335 if( !astISFINITE( *pv ) ) *pv = AST__BAD;
3336 }
3337 }
3338 } else {
3339 SystemChange( AST__VLTOVR, np, values, args, 1, status );
3340 }
3341 break;
3342
3343 /* Relativistic velocity to radio velocity. */
3344 case AST__VLTOVR:
3345 if( forward ) {
3346 pv = values - 1;
3347 for( i = 0; i < np; i++ ) {
3348 pv++;
3349 if( *pv != AST__BAD ) {
3350 temp = AST__C + ( *pv );
3351 if( temp != 0.0 ) {
3352 temp = (AST__C - *pv )/temp;
3353 if( temp >= 0.0 ) {
3354 *pv = AST__C*( 1.0 - sqrt( temp ) );
3355 if( !astISFINITE( *pv ) ) *pv = AST__BAD;
3356 } else {
3357 *pv = AST__BAD;
3358 }
3359 } else {
3360 *pv = AST__BAD;
3361 }
3362 }
3363 }
3364 } else {
3365 SystemChange( AST__VRTOVL, np, values, args, 1, status );
3366 }
3367 break;
3368
3369 /* Optical velocity to relativistic velocity */
3370 case AST__VOTOVL:
3371 if( forward ) {
3372 pv = values - 1;
3373 for( i = 0; i < np; i++ ) {
3374 pv++;
3375 if( *pv != AST__BAD ) {
3376 temp = 1.0 + ( *pv )/AST__C;
3377 temp *= temp;
3378 *pv = AST__C*( temp - 1.0 )/( temp + 1.0 );
3379 if( !astISFINITE( *pv ) ) *pv = AST__BAD;
3380 }
3381 }
3382 } else {
3383 SystemChange( AST__VLTOVO, np, values, args, 1, status );
3384 }
3385 break;
3386
3387 /* Relativistic velocity to optical velocity. */
3388 case AST__VLTOVO:
3389 if( forward ) {
3390 pv = values - 1;
3391 for( i = 0; i < np; i++ ) {
3392 pv++;
3393 if( *pv != AST__BAD ) {
3394 temp = AST__C - *pv;
3395 if( temp != 0.0 ) {
3396 temp = (AST__C + *pv )/temp;
3397 if( temp >= 0.0 ) {
3398 *pv = AST__C*( sqrt( temp ) - 1.0 );
3399 if( !astISFINITE( *pv ) ) *pv = AST__BAD;
3400 } else {
3401 *pv = AST__BAD;
3402 }
3403 } else {
3404 *pv = AST__BAD;
3405 }
3406 }
3407 }
3408 } else {
3409 SystemChange( AST__VOTOVL, np, values, args, 1, status );
3410 }
3411 break;
3412
3413 /* Redshift to relativistic velocity */
3414 case AST__ZOTOVL:
3415 if( forward ) {
3416 pv = values - 1;
3417 for( i = 0; i < np; i++ ) {
3418 pv++;
3419 if( *pv != AST__BAD ) {
3420 temp = 1.0 + ( *pv );
3421 temp *= temp;
3422 *pv = AST__C*( temp - 1.0 )/( temp + 1.0 );
3423 if( !astISFINITE( *pv ) ) *pv = AST__BAD;
3424 }
3425 }
3426 } else {
3427 SystemChange( AST__VLTOZO, np, values, args, 1, status );
3428 }
3429 break;
3430
3431 /* Relativistic velocity to redshift. */
3432 case AST__VLTOZO:
3433 if( forward ) {
3434 pv = values - 1;
3435 for( i = 0; i < np; i++ ) {
3436 pv++;
3437 if( *pv != AST__BAD ) {
3438 temp = AST__C - *pv;
3439 if( temp != 0.0 ) {
3440 temp = (AST__C + *pv )/temp;
3441 if( temp >= 0.0 ) {
3442 *pv = sqrt( temp ) - 1.0;
3443 if( !astISFINITE( *pv ) ) *pv = AST__BAD;
3444 } else {
3445 *pv = AST__BAD;
3446 }
3447 } else {
3448 *pv = AST__BAD;
3449 }
3450 }
3451 }
3452 } else {
3453 SystemChange( AST__ZOTOVL, np, values, args, 1, status );
3454 }
3455 break;
3456
3457 /* Beta factor to relativistic velocity */
3458 case AST__BTTOVL:
3459 if( forward ) {
3460 pv = values - 1;
3461 for( i = 0; i < np; i++ ) {
3462 pv++;
3463 if( *pv != AST__BAD ) {
3464 *pv *= AST__C;
3465 }
3466 }
3467 } else {
3468 SystemChange( AST__VLTOBT, np, values, args, 1, status );
3469 }
3470 break;
3471
3472 /* Relativistic velocity to beta factor. */
3473 case AST__VLTOBT:
3474 if( forward ) {
3475 pv = values - 1;
3476 for( i = 0; i < np; i++ ) {
3477 pv++;
3478 if( *pv != AST__BAD ) {
3479 *pv /= AST__C;
3480 }
3481 }
3482 } else {
3483 SystemChange( AST__BTTOVL, np, values, args, 1, status );
3484 }
3485 break;
3486
3487 /* If the supplied code does not represent a change of system, clear
3488 the returned flag. */
3489 default:
3490 result = 0;
3491 }
3492
3493 /* Return the result. */
3494 return result;
3495 }
3496
TopoVel(double ra,double dec,FrameDef * def,int * status)3497 static double TopoVel( double ra, double dec, FrameDef *def, int *status ) {
3498 /*
3499 * Name:
3500 * TopoVel
3501
3502 * Purpose:
3503 * Find the velocity of the observer away from the source.
3504
3505 * Type:
3506 * Private function.
3507
3508 * Synopsis:
3509 * #include "specmap.h"
3510 * double TopoVel( double ra, double dec, FrameDef *def, int *status )
3511
3512 * Class Membership:
3513 * SpecMap method.
3514
3515 * Description:
3516 * This function finds the component of the velocity of the observer away
3517 * from a specified source position, at a given epoch, in the frame of
3518 * rest of the Sun.
3519
3520 * Parameters:
3521 * ra
3522 * The RA (rads, FK5 J2000) of the source.
3523 * dec
3524 * The Dec (rads, FK5 J2000) of the source.
3525 * def
3526 * Pointer to a FrameDef structure which holds the parameters which
3527 * define the frame, together with cached intermediate results.
3528 * status
3529 * Pointer to the inherited status variable.
3530
3531 * Returns:
3532 * The component of the frame's velocity away from the position given by
3533 * "ra" and "dec", in m/s, measured within the Heliographic frame of
3534 * rest. Zero is returned if an error has already occurred.
3535
3536 */
3537
3538 /* Local Variables: */
3539 double deca; /* Apparent DEC */
3540 double raa; /* Apparent RA */
3541 double vobs; /* Velocity of observer relative to earth */
3542 double vearth; /* Velocity of earth realtive to sun */
3543
3544 /* Check the global error status. */
3545 if ( !astOK ) return 0.0;
3546
3547 /* If not already done so, get the parameters defining the transformation
3548 of mean ra and dec to apparent ra and dec, and store in the supplied frame
3549 definition structure. */
3550 if( def->amprms[ 0 ] == AST__BAD ) palMappa( 2000.0, def->epoch,
3551 def->amprms );
3552
3553 /* Convert the source position from mean ra and dec to apparent ra and dec. */
3554 palMapqkz( ra, dec, def->amprms, &raa, &deca );
3555
3556 /* If not already done so, get the local apparent siderial time (in radians)
3557 and store in the supplied frame definition structure. */
3558 if( def->last == AST__BAD ) def->last = palGmst( def->epoch ) +
3559 palEqeqx( def->epoch ) +
3560 def->obslon;
3561
3562 /* Get the component away from the source, of the velocity of the observer
3563 relative to the centre of the earth (in m/s). */
3564 vobs = 1000.0*Rverot( def->obslat, def->obsalt, raa, deca, def->last,
3565 status );
3566
3567 /* Get the component away from the source, of the velocity of the earth's
3568 centre relative to the Sun, in m/s. */
3569 vearth = GeoVel( ra, dec, def, status );
3570
3571 /* Return the total velocity of the observer away from the source in the
3572 frame of the sun. */
3573 return vobs + vearth;
3574 }
3575
Transform(AstMapping * this,AstPointSet * in,int forward,AstPointSet * out,int * status)3576 static AstPointSet *Transform( AstMapping *this, AstPointSet *in,
3577 int forward, AstPointSet *out, int *status ) {
3578 /*
3579 * Name:
3580 * Transform
3581
3582 * Purpose:
3583 * Apply a SpecMap to transform a set of points.
3584
3585 * Type:
3586 * Private function.
3587
3588 * Synopsis:
3589 * #include "specmap.h"
3590 * AstPointSet *Transform( AstMapping *this, AstPointSet *in,
3591 * int forward, AstPointSet *out, int *status )
3592
3593 * Class Membership:
3594 * SpecMap member function (over-rides the astTransform method inherited
3595 * from the Mapping class).
3596
3597 * Description:
3598 * This function takes a SpecMap and a set of points encapsulated
3599 * in a PointSet and transforms the points so as to perform the
3600 * sequence of spectral coordinate conversions specified by
3601 * previous invocations of astSpecAdd.
3602
3603 * Parameters:
3604 * this
3605 * Pointer to the SpecMap.
3606 * in
3607 * Pointer to the PointSet holding the input coordinate data.
3608 * forward
3609 * A non-zero value indicates that the forward coordinate transformation
3610 * should be applied, while a zero value requests the inverse
3611 * transformation.
3612 * out
3613 * Pointer to a PointSet which will hold the transformed (output)
3614 * coordinate values. A NULL value may also be given, in which case a
3615 * new PointSet will be created by this function.
3616 * status
3617 * Pointer to the inherited status variable.
3618
3619 * Returned Value:
3620 * Pointer to the output (possibly new) PointSet.
3621
3622 * Notes:
3623 * - A null pointer will be returned if this function is invoked with the
3624 * global error status set, or if it should fail for any reason.
3625 * - The number of coordinate values per point in the input PointSet must
3626 * match the number of coordinates for the SpecMap being applied.
3627 * - If an output PointSet is supplied, it must have space for sufficient
3628 * number of points and coordinate values per point to accommodate the
3629 * result. Any excess space will be ignored.
3630 */
3631
3632 /* Local Variables: */
3633 AstPointSet *result; /* Pointer to output PointSet */
3634 AstSpecMap *map; /* Pointer to SpecMap to be applied */
3635 double **ptr_in; /* Pointer to input coordinate data */
3636 double **ptr_out; /* Pointer to output coordinate data */
3637 double *spec; /* Pointer to output spectral axis value array */
3638 double *alpha; /* Pointer to output RA axis value array */
3639 double *beta; /* Pointer to output DEC axis value array */
3640 int cvt; /* Loop counter for conversions */
3641 int end; /* Termination index for conversion loop */
3642 int inc; /* Increment for conversion loop */
3643 int map3d; /* Is the SpecMap 3-dimensional? */
3644 int ncoord_in; /* Number of coordinates per input point */
3645 int npoint; /* Number of points */
3646 int start; /* Starting index for conversion loop */
3647
3648 /* Check the global error status. */
3649 if ( !astOK ) return NULL;
3650
3651 /* Obtain a pointer to the SpecMap. */
3652 map = (AstSpecMap *) this;
3653
3654 /* Apply the parent mapping using the stored pointer to the Transform member
3655 function inherited from the parent Mapping class. This function validates
3656 all arguments and generates an output PointSet if necessary, but does not
3657 actually transform any coordinate values. */
3658 result = (*parent_transform)( this, in, forward, out, status );
3659
3660 /* We will now extend the parent astTransform method by performing the
3661 coordinate conversions needed to generate the output coordinate values. */
3662
3663 /* Determine the numbers of points and coordinates per point from the input
3664 PointSet and obtain pointers for accessing the input and output coordinate
3665 values. */
3666 ncoord_in = astGetNcoord( in );
3667 npoint = astGetNpoint( in );
3668 ptr_in = astGetPoints( in );
3669 ptr_out = astGetPoints( result );
3670
3671 /* Determine whether to apply the forward or inverse transformation, according
3672 to the direction specified and whether the mapping has been inverted. */
3673 if ( astGetInvert( this ) ) forward = !forward;
3674
3675 /* Transform the coordinate values. */
3676 /* -------------------------------- */
3677 /* Use "spec" as a synonym for the array of spectral axis values stored in
3678 the output PointSet. */
3679 if ( astOK ) {
3680 spec = ptr_out[ 0 ];
3681
3682 /* If this is a 3D SpecMap use "alpha" as a synonym for the array of RA axis
3683 values and "beta" as a synonym for the array of DEC axis values stored
3684 in the output PointSet. */
3685 map3d = ( ncoord_in == 3 );
3686 if( map3d ) {
3687 alpha = ptr_out[ 1 ];
3688 beta = ptr_out[ 2 ];
3689 } else {
3690 alpha = NULL;
3691 beta = NULL;
3692 }
3693
3694 /* Initialise the output coordinate values by copying the input ones. */
3695 (void) memcpy( spec, ptr_in[ 0 ], sizeof( double ) * (size_t) npoint );
3696 if( map3d ) {
3697 (void) memcpy( alpha, ptr_in[ 1 ], sizeof( double ) * (size_t) npoint );
3698 (void) memcpy( beta, ptr_in[ 2 ], sizeof( double ) * (size_t) npoint );
3699 }
3700
3701 /* We will loop to apply each spectral coordinate conversion in turn to the
3702 (spec) array. However, if the inverse transformation was requested,
3703 we must loop through these transformations in reverse order, so set up
3704 appropriate limits and an increment to control this loop. */
3705 start = forward ? 0 : map->ncvt - 1;
3706 end = forward ? map->ncvt : -1;
3707 inc = forward ? 1 : -1;
3708
3709 /* Loop through the coordinate conversions in the required order. */
3710 for ( cvt = start; cvt != end; cvt += inc ) {
3711
3712 /* Process conversions which correspond to changes of reference frames. */
3713 if( !FrameChange( map->cvttype[ cvt ], npoint, alpha, beta, spec,
3714 map->cvtargs[ cvt ], forward, status ) ) {
3715
3716 /* If this conversion was not a change of reference frame, it must be a
3717 change of system. */
3718 SystemChange( map->cvttype[ cvt ], npoint, spec,
3719 map->cvtargs[ cvt ], forward, status );
3720 }
3721 }
3722 }
3723
3724 /* If an error has occurred and a new PointSet may have been created, then
3725 clean up by annulling it. In any case, ensure that a NULL result is
3726 returned.*/
3727 if ( !astOK ) {
3728 if ( !out ) result = astAnnul( result );
3729 result = NULL;
3730 }
3731
3732 /* Return a pointer to the output PointSet. */
3733 return result;
3734
3735 }
3736
UserVel(double ra,double dec,FrameDef * def,int * status)3737 static double UserVel( double ra, double dec, FrameDef *def, int *status ) {
3738 /*
3739 * Name:
3740 * UserVel
3741
3742 * Purpose:
3743 * Find the component of the velocity of the user-defined rest-frame
3744 * away from the source.
3745
3746 * Type:
3747 * Private function.
3748
3749 * Synopsis:
3750 * #include "specmap.h"
3751 * double UserVel( double ra, double dec, FrameDef *def, int *status )
3752
3753 * Class Membership:
3754 * SpecMap method.
3755
3756 * Description:
3757 * This function finds the component of the velocity of the user-defined
3758 * rest-frame away from a specified position. The magnitude and direction
3759 * of the rest-frames velocity are defined within the supplied "def"
3760 * structure. The user-defined rest-frame is typically used to represent
3761 * the velocity of the source within the heliocentric rest-frame.
3762
3763 * Parameters:
3764 * ra
3765 * The RA (rads, FK5 J2000) of the source.
3766 * dec
3767 * The Dec (rads, FK5 J2000) of the source.
3768 * def
3769 * Pointer to a FrameDef structure which holds the parameters which
3770 * define the frame, together with cached intermediate results.
3771 * status
3772 * Pointer to the inherited status variable.
3773
3774 * Returns:
3775 * The component of the frame's velocity away from the position given by
3776 * "ra" and "dec", in m/s, measured within the Heliographic frame of
3777 * rest. Zero is returned if an error has already occurred.
3778
3779 * Notes:
3780 * - The direction of the user velocity is given by def->refra and
3781 * def->refdec (an FK5 J2000 position). The maginitude of the velocity
3782 * is given by def->veluser, in m/s, positive when the source is moving
3783 * away from the observer towards def->refra, def->refdec, and given
3784 * with respect to the heliocentric rest-frame.
3785
3786 */
3787
3788 /* Local Variables: */
3789 double vb[ 3 ]; /* Source position vector */
3790
3791 /* Check the global error status. */
3792 if ( !astOK ) return 0.0;
3793
3794 /* If not already done so, express the user velocity in the form of a
3795 J2000.0 x,y,z vector. */
3796 if( def->vuser[ 0 ] == AST__BAD ) {
3797 def->vuser[ 0 ] = def->veluser*cos( def->refra )*cos( def->refdec );
3798 def->vuser[ 1 ] = def->veluser*sin( def->refra )*cos( def->refdec );
3799 def->vuser[ 2 ] = def->veluser*sin( def->refdec );
3800 }
3801
3802 /* Convert given J2000 RA,Dec to x,y,z. */
3803 palDcs2c( ra, dec, vb );
3804
3805 /* Return the dot product with the user velocity. Invert it to get the
3806 velocity towards the observer (the def->veluser value is supposed to be
3807 positive if the source is moving away from the observer). */
3808 return -palDvdv( def->vuser, vb );
3809 }
3810
3811 /* Copy constructor. */
3812 /* ----------------- */
Copy(const AstObject * objin,AstObject * objout,int * status)3813 static void Copy( const AstObject *objin, AstObject *objout, int *status ) {
3814 /*
3815 * Name:
3816 * Copy
3817
3818 * Purpose:
3819 * Copy constructor for SpecMap objects.
3820
3821 * Type:
3822 * Private function.
3823
3824 * Synopsis:
3825 * void Copy( const AstObject *objin, AstObject *objout, int *status )
3826
3827 * Description:
3828 * This function implements the copy constructor for SpecMap objects.
3829
3830 * Parameters:
3831 * objin
3832 * Pointer to the object to be copied.
3833 * objout
3834 * Pointer to the object being constructed.
3835 * status
3836 * Pointer to the inherited status variable.
3837
3838 * Returned Value:
3839 * void
3840
3841 * Notes:
3842 * - This constructor makes a deep copy.
3843 */
3844
3845 /* Local Variables: */
3846 AstSpecMap *in; /* Pointer to input SpecMap */
3847 AstSpecMap *out; /* Pointer to output SpecMap */
3848 int cvt; /* Loop counter for coordinate conversions */
3849
3850 /* Check the global error status. */
3851 if ( !astOK ) return;
3852
3853 /* Obtain pointers to the input and output SpecMap structures. */
3854 in = (AstSpecMap *) objin;
3855 out = (AstSpecMap *) objout;
3856
3857 /* For safety, first clear any references to the input memory from the output
3858 SpecMap. */
3859 out->cvtargs = NULL;
3860 out->cvttype = NULL;
3861
3862 /* Allocate memory for the output array of argument list pointers. */
3863 out->cvtargs = astMalloc( sizeof( double * ) * (size_t) in->ncvt );
3864
3865 /* If necessary, allocate memory and make a copy of the input array of
3866 coordinate conversion codes. */
3867 if ( in->cvttype ) out->cvttype = astStore( NULL, in->cvttype,
3868 sizeof( int )
3869 * (size_t) in->ncvt );
3870
3871 /* If OK, loop through each conversion in the input SpecMap and make a copy of
3872 its argument list, storing the new pointer in the output argument list
3873 array. */
3874 if ( astOK ) {
3875 for ( cvt = 0; cvt < in->ncvt; cvt++ ) {
3876 out->cvtargs[ cvt ] = astStore( NULL, in->cvtargs[ cvt ],
3877 astSizeOf( in->cvtargs[ cvt ] ) );
3878 }
3879
3880 /* If an error occurred while copying the argument lists, loop through the
3881 conversions again and clean up by ensuring that the new memory allocated for
3882 each argument list is freed. */
3883 if ( !astOK ) {
3884 for ( cvt = 0; cvt < in->ncvt; cvt++ ) {
3885 out->cvtargs[ cvt ] = astFree( out->cvtargs[ cvt ] );
3886 }
3887 }
3888 }
3889
3890 /* If an error occurred, free all other memory allocated above. */
3891 if ( !astOK ) {
3892 out->cvtargs = astFree( out->cvtargs );
3893 out->cvttype = astFree( out->cvttype );
3894 }
3895 }
3896
3897 /* Destructor. */
3898 /* ----------- */
Delete(AstObject * obj,int * status)3899 static void Delete( AstObject *obj, int *status ) {
3900 /*
3901 * Name:
3902 * Delete
3903
3904 * Purpose:
3905 * Destructor for SpecMap objects.
3906
3907 * Type:
3908 * Private function.
3909
3910 * Synopsis:
3911 * void Delete( AstObject *obj, int *status )
3912
3913 * Description:
3914 * This function implements the destructor for SpecMap objects.
3915
3916 * Parameters:
3917 * obj
3918 * Pointer to the object to be deleted.
3919 * status
3920 * Pointer to the inherited status variable.
3921
3922 * Returned Value:
3923 * void
3924
3925 * Notes:
3926 * This function attempts to execute even if the global error status is
3927 * set.
3928 */
3929
3930 /* Local Variables: */
3931 AstSpecMap *this; /* Pointer to SpecMap */
3932 int cvt; /* Loop counter for coordinate conversions */
3933
3934 /* Obtain a pointer to the SpecMap structure. */
3935 this = (AstSpecMap *) obj;
3936
3937 /* Loop to free the memory containing the argument list for each coordinate
3938 conversion. */
3939 for ( cvt = 0; cvt < this->ncvt; cvt++ ) {
3940 this->cvtargs[ cvt ] = astFree( this->cvtargs[ cvt ] );
3941 }
3942
3943 /* Free the memory holding the array of conversion types and the array of
3944 argument list pointers. */
3945 this->cvtargs = astFree( this->cvtargs );
3946 this->cvttype = astFree( this->cvttype );
3947 }
3948
3949 /* Dump function. */
3950 /* -------------- */
Dump(AstObject * this_object,AstChannel * channel,int * status)3951 static void Dump( AstObject *this_object, AstChannel *channel, int *status ) {
3952 /*
3953 * Name:
3954 * Dump
3955
3956 * Purpose:
3957 * Dump function for SpecMap objects.
3958
3959 * Type:
3960 * Private function.
3961
3962 * Synopsis:
3963 * #include "specmap.h"
3964 * void Dump( AstObject *this, AstChannel *channel, int *status )
3965
3966 * Description:
3967 * This function implements the Dump function which writes out data
3968 * for the SpecMap class to an output Channel.
3969
3970 * Parameters:
3971 * this
3972 * Pointer to the SpecMap whose data are being written.
3973 * channel
3974 * Pointer to the Channel to which the data are being written.
3975 * status
3976 * Pointer to the inherited status variable.
3977 */
3978
3979 /* Local Constants: */
3980 #define KEY_LEN 50 /* Maximum length of a keyword */
3981
3982 /* Local Variables: */
3983 AstSpecMap *this; /* Pointer to the SpecMap structure */
3984 char key[ KEY_LEN + 1 ]; /* Buffer for keyword string */
3985 const char *argdesc[ MAX_ARGS ]; /* Pointers to argument descriptions */
3986 const char *comment; /* Pointer to comment string */
3987 const char *sval; /* Pointer to string value */
3988 int argdec; /* Index of DEC argument */
3989 int argra; /* Index of RA argument */
3990 int iarg; /* Loop counter for arguments */
3991 int icvt; /* Loop counter for conversion steps */
3992 int ival; /* Integer value */
3993 int nargs; /* Number of user-supplied arguments */
3994 int szargs; /* Number of stored arguments */
3995 int set; /* Attribute value set? */
3996
3997 /* Check the global error status. */
3998 if ( !astOK ) return;
3999
4000 /* Obtain a pointer to the SpecMap structure. */
4001 this = (AstSpecMap *) this_object;
4002
4003 /* Write out values representing the instance variables for the SpecMap
4004 class. Accompany these with appropriate comment strings, possibly
4005 depending on the values being written.*/
4006
4007 /* In the case of attributes, we first use the appropriate (private)
4008 Test... member function to see if they are set. If so, we then use
4009 the (private) Get... function to obtain the value to be written
4010 out.
4011
4012 For attributes which are not set, we use the astGet... method to
4013 obtain the value instead. This will supply a default value
4014 (possibly provided by a derived class which over-rides this method)
4015 which is more useful to a human reader as it corresponds to the
4016 actual default attribute value. Since "set" will be zero, these
4017 values are for information only and will not be read back. */
4018
4019 /* Number of conversion steps. */
4020 /* --------------------------- */
4021 /* Regard this as "set" if it is non-zero. */
4022 ival = this->ncvt;
4023 set = ( ival != 0 );
4024 astWriteInt( channel, "Nspec", set, 0, ival, "Number of conversion steps" );
4025
4026 /* Write out data for each conversion step... */
4027 for ( icvt = 0; icvt < this->ncvt; icvt++ ) {
4028
4029 /* Conversion type. */
4030 /* ---------------- */
4031 /* Change each conversion type code into an equivalent string and
4032 obtain associated descriptive information. If the conversion code
4033 was not recognised, report an error and give up. */
4034 if ( astOK ) {
4035 sval = CvtString( this->cvttype[ icvt ], &comment, &argra, &argdec,
4036 &nargs, &szargs, argdesc, status );
4037 if ( astOK && !sval ) {
4038 astError( AST__SPCIN,
4039 "astWrite(%s): Corrupt %s contains invalid SpecMap "
4040 "spectral coordinate conversion code (%d).", status,
4041 astGetClass( channel ), astGetClass( this ),
4042 (int) this->cvttype[ icvt ] );
4043 break;
4044 }
4045
4046 /* Create an appropriate keyword and write out the conversion code
4047 information. */
4048 (void) sprintf( key, "Spec%d", icvt + 1 );
4049 astWriteString( channel, key, 1, 1, sval, comment );
4050
4051 /* Write out data for each conversion argument... */
4052 for ( iarg = 0; iarg < szargs; iarg++ ) {
4053
4054 /* Arguments. */
4055 /* ---------- */
4056 /* Create an appropriate keyword and write out the argument value,
4057 accompanied by the descriptive comment obtained above. */
4058 if( this->cvtargs[ icvt ][ iarg ] != AST__BAD ) {
4059 (void) sprintf( key, "Spec%d%c", icvt + 1, ALPHABET[ iarg ] );
4060 astWriteDouble( channel, key, 1, 1, this->cvtargs[ icvt ][ iarg ],
4061 argdesc[ iarg ] );
4062 }
4063 }
4064
4065 /* Quit looping if an error occurs. */
4066 if ( !astOK ) break;
4067 }
4068 }
4069
4070 /* Undefine macros local to this function. */
4071 #undef KEY_LEN
4072 }
4073
4074 /* Standard class functions. */
4075 /* ========================= */
4076 /* Implement the astIsASpecMap and astCheckSpecMap functions using the macros
4077 defined for this purpose in the "object.h" header file. */
astMAKE_ISA(SpecMap,Mapping)4078 astMAKE_ISA(SpecMap,Mapping)
4079 astMAKE_CHECK(SpecMap)
4080
4081 AstSpecMap *astSpecMap_( int nin, int flags, const char *options, int *status, ...) {
4082 /*
4083 *++
4084 * Name:
4085 c astSpecMap
4086 f AST_SPECMAP
4087
4088 * Purpose:
4089 * Create a SpecMap.
4090
4091 * Type:
4092 * Public function.
4093
4094 * Synopsis:
4095 c #include "specmap.h"
4096 c AstSpecMap *astSpecMap( int nin, int flags, const char *options, ... )
4097 f RESULT = AST_SPECMAP( NIN, FLAGS, OPTIONS, STATUS )
4098
4099 * Class Membership:
4100 * SpecMap constructor.
4101
4102 * Description:
4103 * This function creates a new SpecMap and optionally initialises
4104 * its attributes.
4105 *
4106 * An SpecMap is a specialised form of Mapping which can be used to
4107 * represent a sequence of conversions between standard spectral
4108 * coordinate systems. This includes conversions between frequency,
4109 * wavelength, and various forms of velocity, as well as conversions
4110 * between different standards of rest.
4111 *
4112 * When a SpecMap is first created, it simply performs a unit
4113 c (null) Mapping. Using the astSpecAdd function,
4114 f (null) Mapping. Using the AST_SPECADD routine,
4115 * a series of coordinate conversion steps may then be added, selected
4116 * from the list of supported conversions. This allows multi-step
4117 * conversions between a variety of spectral coordinate systems to
4118 * be assembled out of the building blocks provided by this class.
4119 *
4120 * For details of the individual coordinate conversions available,
4121 c see the description of the astSpecAdd function.
4122 f see the description of the AST_SPECADD routine.
4123 *
4124 * Conversions are available to transform between standards of rest.
4125 * Such conversions need to know the source position as an RA and DEC.
4126 * This information can be supplied in the form of parameters for
4127 * the relevant conversions, in which case the SpecMap is 1-dimensional,
4128 * simply transforming the spectral axis values. This means that the
4129 * same source position will always be used by the SpecMap. However, this
4130 * may not be appropriate for an accurate description of a 3-D spectral
4131 * cube, where changes of spatial position can produce significant
4132 * changes in the Doppler shift introduced when transforming between
4133 * standards of rest. For this situation, a 3-dimensional SpecMap can
4134 * be created in which axes 2 and 3 correspond to the source RA and DEC
4135 * The SpecMap simply copies values for axes 2 and 3 from input to
4136 * output).
4137
4138 * Parameters:
4139 c nin
4140 f NIN = INTEGER (Given)
4141 * The number of inputs to the Mapping (this will also equal the
4142 * number of outputs). This value must be either 1 or 3. In either
4143 * case, the first input and output correspoindis the spectral axis.
4144 * For a 3-axis SpecMap, the second and third axes give the RA and
4145 * DEC (J2000 FK5) of the source. This positional information is
4146 * used by conversions which transform between standards of rest,
4147 * and replaces the "RA" and "DEC" arguments for the individual
4148 * conversions listed in description of the "SpecAdd"
4149 c function.
4150 f routine.
4151 c flags
4152 f FLAGS = INTEGER (Given)
4153 c This parameter is reserved for future use and should currently
4154 f This argument is reserved for future use and should currently
4155 * always be set to zero.
4156 c options
4157 f OPTIONS = CHARACTER * ( * ) (Given)
4158 c Pointer to a null-terminated string containing an optional
4159 c comma-separated list of attribute assignments to be used for
4160 c initialising the new SpecMap. The syntax used is identical to
4161 c that for the astSet function and may include "printf" format
4162 c specifiers identified by "%" symbols in the normal way.
4163 c If no initialisation is required, a zero-length string may be
4164 c supplied.
4165 f A character string containing an optional comma-separated
4166 f list of attribute assignments to be used for initialising the
4167 f new SpecMap. The syntax used is identical to that for the
4168 f AST_SET routine. If no initialisation is required, a blank
4169 f value may be supplied.
4170 c ...
4171 c If the "options" string contains "%" format specifiers, then
4172 c an optional list of additional arguments may follow it in
4173 c order to supply values to be substituted for these
4174 c specifiers. The rules for supplying these are identical to
4175 c those for the astSet function (and for the C "printf"
4176 c function).
4177 f STATUS = INTEGER (Given and Returned)
4178 f The global status.
4179
4180 * Returned Value:
4181 c astSpecMap()
4182 f AST_SPECMAP = INTEGER
4183 * A pointer to the new SpecMap.
4184
4185 * Notes:
4186 * - The nature and units of the coordinate values supplied for the
4187 * first input (i.e. the spectral input) of a SpecMap must be appropriate
4188 * to the first conversion step applied by the SpecMap. For instance, if
4189 * the first conversion step is "FRTOVL" (frequency to relativistic
4190 * velocity), then the coordinate values for the first input should
4191 * be frequency in units of Hz. Similarly, the nature and units of the
4192 * coordinate values returned by a SpecMap will be determined by the
4193 * last conversion step applied by the SpecMap. For instance, if the
4194 * last conversion step is "VLTOVO" (relativistic velocity to optical
4195 * velocity), then the coordinate values for the first output will be optical
4196 * velocity in units of metres per second. See the description of the
4197 c astSpecAdd function for the units expected and returned by each
4198 f AST_SPECADD routine for the units expected and returned by each
4199 * conversion.
4200 * - A null Object pointer (AST__NULL) will be returned if this
4201 c function is invoked with the AST error status set, or if it
4202 f function is invoked with STATUS set to an error value, or if it
4203 * should fail for any reason.
4204 *--
4205 */
4206
4207 /* Local Variables: */
4208 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
4209 AstSpecMap *new; /* Pointer to the new SpecMap */
4210 va_list args; /* Variable argument list */
4211
4212 /* Get a pointer to the thread specific global data structure. */
4213 astGET_GLOBALS(NULL);
4214
4215 /* Check the global status. */
4216 if ( !astOK ) return NULL;
4217
4218 /* Initialise the SpecMap, allocating memory and initialising the virtual
4219 function table as well if necessary. */
4220 new = astInitSpecMap( NULL, sizeof( AstSpecMap ), !class_init, &class_vtab,
4221 "SpecMap", nin, flags );
4222
4223 /* If successful, note that the virtual function table has been initialised. */
4224 if ( astOK ) {
4225 class_init = 1;
4226
4227 /* Obtain the variable argument list and pass it along with the options string
4228 to the astVSet method to initialise the new SpecMap's attributes. */
4229 va_start( args, status );
4230 astVSet( new, options, NULL, args );
4231 va_end( args );
4232
4233 /* If an error occurred, clean up by deleting the new object. */
4234 if ( !astOK ) new = astDelete( new );
4235 }
4236
4237 /* Return a pointer to the new SpecMap. */
4238 return new;
4239 }
4240
astSpecMapId_(int nin,int flags,const char * options,...)4241 AstSpecMap *astSpecMapId_( int nin, int flags, const char *options, ... ) {
4242 /*
4243 * Name:
4244 * astSpecMapId_
4245
4246 * Purpose:
4247 * Create a SpecMap.
4248
4249 * Type:
4250 * Private function.
4251
4252 * Synopsis:
4253 * #include "specmap.h"
4254 * AstSpecMap *astSpecMapId_( int nin, int flags, const char *options, ... )
4255
4256 * Class Membership:
4257 * SpecMap constructor.
4258
4259 * Description:
4260 * This function implements the external (public) interface to the
4261 * astSpecMap constructor function. It returns an ID value (instead
4262 * of a true C pointer) to external users, and must be provided
4263 * because astSpecMap_ has a variable argument list which cannot be
4264 * encapsulated in a macro (where this conversion would otherwise
4265 * occur).
4266 *
4267 * The variable argument list also prevents this function from
4268 * invoking astSpecMap_ directly, so it must be a re-implementation
4269 * of it in all respects, except for the final conversion of the
4270 * result to an ID value.
4271
4272 * Parameters:
4273 * As for astSpecMap_.
4274
4275 * Returned Value:
4276 * The ID value associated with the new SpecMap.
4277 */
4278
4279 /* Local Variables: */
4280 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
4281 AstSpecMap *new; /* Pointer to the new SpecMap */
4282 va_list args; /* Variable argument list */
4283
4284 int *status; /* Pointer to inherited status value */
4285
4286 /* Get a pointer to the inherited status value. */
4287 status = astGetStatusPtr;
4288
4289 /* Get a pointer to the thread specific global data structure. */
4290 astGET_GLOBALS(NULL);
4291
4292 /* Check the global status. */
4293 if ( !astOK ) return NULL;
4294
4295 /* Initialise the SpecMap, allocating memory and initialising the virtual
4296 function table as well if necessary. */
4297 new = astInitSpecMap( NULL, sizeof( AstSpecMap ), !class_init, &class_vtab,
4298 "SpecMap", nin, flags );
4299
4300 /* If successful, note that the virtual function table has been initialised. */
4301 if ( astOK ) {
4302 class_init = 1;
4303
4304 /* Obtain the variable argument list and pass it along with the options string
4305 to the astVSet method to initialise the new SpecMap's attributes. */
4306 va_start( args, options );
4307 astVSet( new, options, NULL, args );
4308 va_end( args );
4309
4310 /* If an error occurred, clean up by deleting the new object. */
4311 if ( !astOK ) new = astDelete( new );
4312 }
4313
4314 /* Return an ID value for the new SpecMap. */
4315 return astMakeId( new );
4316 }
4317
astInitSpecMap_(void * mem,size_t size,int init,AstSpecMapVtab * vtab,const char * name,int nin,int flags,int * status)4318 AstSpecMap *astInitSpecMap_( void *mem, size_t size, int init,
4319 AstSpecMapVtab *vtab, const char *name,
4320 int nin, int flags, int *status ) {
4321 /*
4322 *+
4323 * Name:
4324 * astInitSpecMap
4325
4326 * Purpose:
4327 * Initialise a SpecMap.
4328
4329 * Type:
4330 * Protected function.
4331
4332 * Synopsis:
4333 * #include "specmap.h"
4334 * AstSpecMap *astInitSpecMap( void *mem, size_t size, int init,
4335 * AstSpecMapVtab *vtab, const char *name,
4336 * int nin, int flags )
4337
4338 * Class Membership:
4339 * SpecMap initialiser.
4340
4341 * Description:
4342 * This function is provided for use by class implementations to initialise
4343 * a new SpecMap object. It allocates memory (if necessary) to accommodate
4344 * the SpecMap plus any additional data associated with the derived class.
4345 * It then initialises a SpecMap structure at the start of this memory. If
4346 * the "init" flag is set, it also initialises the contents of a virtual
4347 * function table for a SpecMap at the start of the memory passed via the
4348 * "vtab" parameter.
4349
4350 * Parameters:
4351 * mem
4352 * A pointer to the memory in which the SpecMap is to be initialised.
4353 * This must be of sufficient size to accommodate the SpecMap data
4354 * (sizeof(SpecMap)) plus any data used by the derived class. If a value
4355 * of NULL is given, this function will allocate the memory itself using
4356 * the "size" parameter to determine its size.
4357 * size
4358 * The amount of memory used by the SpecMap (plus derived class data).
4359 * This will be used to allocate memory if a value of NULL is given for
4360 * the "mem" parameter. This value is also stored in the SpecMap
4361 * structure, so a valid value must be supplied even if not required for
4362 * allocating memory.
4363 * init
4364 * A logical flag indicating if the SpecMap's virtual function table is
4365 * to be initialised. If this value is non-zero, the virtual function
4366 * table will be initialised by this function.
4367 * vtab
4368 * Pointer to the start of the virtual function table to be associated
4369 * with the new SpecMap.
4370 * name
4371 * Pointer to a constant null-terminated character string which contains
4372 * the name of the class to which the new object belongs (it is this
4373 * pointer value that will subsequently be returned by the astClass
4374 * method).
4375 * nin
4376 * The number of inputs and outputs for the SpecMap (either 1 or 3).
4377 * flags
4378 * This parameter is reserved for future use. It is currently ignored.
4379
4380 * Returned Value:
4381 * A pointer to the new SpecMap.
4382
4383 * Notes:
4384 * - A null pointer will be returned if this function is invoked with the
4385 * global error status set, or if it should fail for any reason.
4386 *-
4387 */
4388
4389 /* Local Variables: */
4390 AstSpecMap *new; /* Pointer to the new SpecMap */
4391
4392 /* Check the global status. */
4393 if ( !astOK ) return NULL;
4394
4395 /* Check nin is OK (1 or 3). */
4396 if( nin != 1 && nin != 3 ) {
4397 astError( AST__BADNI, "astInitSpecMap(SpecMap): Supplied number of "
4398 "SpecMap axes (%d) is illegal; it should be 1 or 2. ", status,
4399 nin );
4400 }
4401
4402 /* If necessary, initialise the virtual function table. */
4403 if ( init ) astInitSpecMapVtab( vtab, name );
4404
4405 /* Initialise a 1D Mapping structure (the parent class) as the first component
4406 within the SpecMap structure, allocating memory if necessary. Specify that
4407 the Mapping should be defined in both the forward and inverse directions. */
4408 new = (AstSpecMap *) astInitMapping( mem, size, 0,
4409 (AstMappingVtab *) vtab, name,
4410 nin, nin, 1, 1 );
4411
4412 if ( astOK ) {
4413
4414 /* Initialise the SpecMap data. */
4415 /* --------------------------- */
4416 /* The initial state is with no conversions set, in which condition the
4417 SpecMap simply implements a unit mapping. */
4418 new->ncvt = 0;
4419 new->cvtargs = NULL;
4420 new->cvttype = NULL;
4421
4422 /* If an error occurred, clean up by deleting the new object. */
4423 if ( !astOK ) new = astDelete( new );
4424 }
4425
4426 /* Return a pointer to the new object. */
4427 return new;
4428 }
4429
astLoadSpecMap_(void * mem,size_t size,AstSpecMapVtab * vtab,const char * name,AstChannel * channel,int * status)4430 AstSpecMap *astLoadSpecMap_( void *mem, size_t size,
4431 AstSpecMapVtab *vtab, const char *name,
4432 AstChannel *channel, int *status ) {
4433 /*
4434 *+
4435 * Name:
4436 * astLoadSpecMap
4437
4438 * Purpose:
4439 * Load a SpecMap.
4440
4441 * Type:
4442 * Protected function.
4443
4444 * Synopsis:
4445 * #include "specmap.h"
4446 * AstSpecMap *astLoadSpecMap( void *mem, size_t size,
4447 * AstSpecMapVtab *vtab, const char *name,
4448 * AstChannel *channel )
4449
4450 * Class Membership:
4451 * SpecMap loader.
4452
4453 * Description:
4454 * This function is provided to load a new SpecMap using data read
4455 * from a Channel. It first loads the data used by the parent class
4456 * (which allocates memory if necessary) and then initialises a
4457 * SpecMap structure in this memory, using data read from the input
4458 * Channel.
4459 *
4460 * If the "init" flag is set, it also initialises the contents of a
4461 * virtual function table for a SpecMap at the start of the memory
4462 * passed via the "vtab" parameter.
4463
4464
4465 * Parameters:
4466 * mem
4467 * A pointer to the memory into which the SpecMap is to be
4468 * loaded. This must be of sufficient size to accommodate the
4469 * SpecMap data (sizeof(SpecMap)) plus any data used by derived
4470 * classes. If a value of NULL is given, this function will
4471 * allocate the memory itself using the "size" parameter to
4472 * determine its size.
4473 * size
4474 * The amount of memory used by the SpecMap (plus derived class
4475 * data). This will be used to allocate memory if a value of
4476 * NULL is given for the "mem" parameter. This value is also
4477 * stored in the SpecMap structure, so a valid value must be
4478 * supplied even if not required for allocating memory.
4479 *
4480 * If the "vtab" parameter is NULL, the "size" value is ignored
4481 * and sizeof(AstSpecMap) is used instead.
4482 * vtab
4483 * Pointer to the start of the virtual function table to be
4484 * associated with the new SpecMap. If this is NULL, a pointer to
4485 * the (static) virtual function table for the SpecMap class is
4486 * used instead.
4487 * name
4488 * Pointer to a constant null-terminated character string which
4489 * contains the name of the class to which the new object
4490 * belongs (it is this pointer value that will subsequently be
4491 * returned by the astGetClass method).
4492 *
4493 * If the "vtab" parameter is NULL, the "name" value is ignored
4494 * and a pointer to the string "SpecMap" is used instead.
4495
4496 * Returned Value:
4497 * A pointer to the new SpecMap.
4498
4499 * Notes:
4500 * - A null pointer will be returned if this function is invoked
4501 * with the global error status set, or if it should fail for any
4502 * reason.
4503 *-
4504 */
4505
4506 /* Local Constants: */
4507 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
4508 #define KEY_LEN 50 /* Maximum length of a keyword */
4509
4510 /* Local Variables: */
4511 AstSpecMap *new; /* Pointer to the new SpecMap */
4512 char *sval; /* Pointer to string value */
4513 char key[ KEY_LEN + 1 ]; /* Buffer for keyword string */
4514 const char *argdesc[ MAX_ARGS ]; /* Pointers to argument descriptions */
4515 const char *comment; /* Pointer to comment string */
4516 int argdec; /* Index of DEC argument */
4517 int argra; /* Index of RA argument */
4518 int iarg; /* Loop counter for arguments */
4519 int icvt; /* Loop counter for conversion steps */
4520 int nargs; /* Number of user-supplied arguments */
4521 int szargs; /* Number of stored arguments */
4522
4523 /* Get a pointer to the thread specific global data structure. */
4524 astGET_GLOBALS(channel);
4525
4526 /* Initialise. */
4527 new = NULL;
4528
4529 /* Check the global error status. */
4530 if ( !astOK ) return new;
4531
4532 /* If a NULL virtual function table has been supplied, then this is
4533 the first loader to be invoked for this SpecMap. In this case the
4534 SpecMap belongs to this class, so supply appropriate values to be
4535 passed to the parent class loader (and its parent, etc.). */
4536 if ( !vtab ) {
4537 size = sizeof( AstSpecMap );
4538 vtab = &class_vtab;
4539 name = "SpecMap";
4540
4541 /* If required, initialise the virtual function table for this class. */
4542 if ( !class_init ) {
4543 astInitSpecMapVtab( vtab, name );
4544 class_init = 1;
4545 }
4546 }
4547
4548 /* Invoke the parent class loader to load data for all the ancestral
4549 classes of the current one, returning a pointer to the resulting
4550 partly-built SpecMap. */
4551 new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name,
4552 channel );
4553
4554 if ( astOK ) {
4555
4556 /* Read input data. */
4557 /* ================ */
4558 /* Request the input Channel to read all the input data appropriate to
4559 this class into the internal "values list". */
4560 astReadClassData( channel, "SpecMap" );
4561
4562 /* Now read each individual data item from this list and use it to
4563 initialise the appropriate instance variable(s) for this class. */
4564
4565 /* In the case of attributes, we first read the "raw" input value,
4566 supplying the "unset" value as the default. If a "set" value is
4567 obtained, we then use the appropriate (private) Set... member
4568 function to validate and set the value properly. */
4569
4570 /* Number of conversion steps. */
4571 /* --------------------------- */
4572 /* Read the number of conversion steps and allocate memory to hold
4573 data for each step. */
4574 new->ncvt = astReadInt( channel, "nspec", 0 );
4575 if ( new->ncvt < 0 ) new->ncvt = 0;
4576 new->cvttype = astMalloc( sizeof( int ) * (size_t) new->ncvt );
4577 new->cvtargs = astMalloc( sizeof( double * ) * (size_t) new->ncvt );
4578
4579 /* If an error occurred, ensure that all allocated memory is freed. */
4580 if ( !astOK ) {
4581 new->cvttype = astFree( new->cvttype );
4582 new->cvtargs = astFree( new->cvtargs );
4583
4584 /* Otherwise, initialise the argument pointer array. */
4585 } else {
4586 for ( icvt = 0; icvt < new->ncvt; icvt++ ) {
4587 new->cvtargs[ icvt ] = NULL;
4588 }
4589
4590 /* Read in data for each conversion step... */
4591 for ( icvt = 0; icvt < new->ncvt; icvt++ ) {
4592
4593 /* Conversion type. */
4594 /* ---------------- */
4595 /* Create an appropriate keyword and read the string representation of
4596 the conversion type. */
4597 (void) sprintf( key, "spec%d", icvt + 1 );
4598 sval = astReadString( channel, key, NULL );
4599
4600 /* If no value was read, report an error. */
4601 if ( astOK ) {
4602 if ( !sval ) {
4603 astError( AST__BADIN,
4604 "astRead(%s): A spectral coordinate conversion "
4605 "type is missing from the input SpecMap data.", status,
4606 astGetClass( channel ) );
4607
4608 /* Otherwise, convert the string representation into the required
4609 conversion type code. */
4610 } else {
4611 new->cvttype[ icvt ] = CvtCode( sval, status );
4612
4613 /* If the string was not recognised, report an error. */
4614 if ( new->cvttype[ icvt ] == AST__SPEC_NULL ) {
4615 astError( AST__BADIN,
4616 "astRead(%s): Invalid spectral conversion "
4617 "type \"%s\" in SpecMap data.", status,
4618 astGetClass( channel ), sval );
4619 }
4620 }
4621
4622 /* Free the memory holding the string value. */
4623 sval = astFree( sval );
4624 }
4625
4626 /* Obtain the number of arguments associated with the conversion and
4627 allocate memory to hold them. */
4628 (void) CvtString( new->cvttype[ icvt ], &comment, &argra,
4629 &argdec, &nargs, &szargs, argdesc, status );
4630 new->cvtargs[ icvt ] = astMalloc( sizeof( double ) *
4631 (size_t) szargs );
4632
4633 /* Read in data for each argument... */
4634 if ( astOK ) {
4635 for ( iarg = 0; iarg < szargs; iarg++ ) {
4636
4637 /* Arguments. */
4638 /* ---------- */
4639 /* Create an appropriate keyword and read each argument value. */
4640 (void) sprintf( key, "spec%d%c", icvt + 1, ALPHABET[ iarg ] );
4641 new->cvtargs[ icvt ][ iarg ] = astReadDouble( channel, key,
4642 AST__BAD );
4643 }
4644 }
4645
4646 /* Quit looping if an error occurs. */
4647 if ( !astOK ) break;
4648 }
4649 }
4650
4651 /* If an error occurred, clean up by deleting the new SpecMap. */
4652 if ( !astOK ) new = astDelete( new );
4653 }
4654
4655 /* Return the new SpecMap pointer. */
4656 return new;
4657
4658 /* Undefine macros local to this function. */
4659 #undef KEY_LEN
4660 }
4661
4662 /* Virtual function interfaces. */
4663 /* ============================ */
4664 /* These provide the external interface to the virtual functions defined by
4665 this class. Each simply checks the global error status and then locates and
4666 executes the appropriate member function, using the function pointer stored
4667 in the object's virtual function table (this pointer is located using the
4668 astMEMBER macro defined in "object.h").
4669
4670 Note that the member function may not be the one defined here, as it may
4671 have been over-ridden by a derived class. However, it should still have the
4672 same interface. */
astSpecAdd_(AstSpecMap * this,const char * cvt,const double args[],int * status)4673 void astSpecAdd_( AstSpecMap *this, const char *cvt, const double args[], int *status ) {
4674 if ( !astOK ) return;
4675 (**astMEMBER(this,SpecMap,SpecAdd))( this, cvt, args, status );
4676 }
4677
4678
4679
4680
4681