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