1 /*
2 *class++
3 * Name:
4 * DssMap
5
6 * Purpose:
7 * Map points using a Digitised Sky Survey plate solution.
8
9 * Constructor Function:
10 * The DssMap class does not have a constructor function. A DssMap
11 * is created only as a by-product of reading a FrameSet (using
12 c astRead) from a FitsChan which contains FITS header cards
13 f AST_READ) from a FitsChan which contains FITS header cards
14 * describing a DSS plate solution, and whose Encoding attribute is
15 * set to "DSS". The result of such a read, if successful, is a
16 * FrameSet whose base and current Frames are related by a DssMap.
17
18 * Description:
19 * The DssMap class implements a Mapping which transforms between
20 * 2-dimensional pixel coordinates and an equatorial sky coordinate
21 * system (right ascension and declination) using a Digitised Sky
22 * Survey (DSS) astrometric plate solution.
23 *
24 * The input coordinates are pixel numbers along the first and
25 * second dimensions of an image, where the centre of the first
26 * pixel is located at (1,1) and the spacing between pixel centres
27 * is unity.
28 *
29 * The output coordinates are right ascension and declination in
30 * radians. The celestial coordinate system used (FK4, FK5, etc.)
31 * is unspecified, and will usually be indicated by appropriate
32 * keywords in a FITS header.
33
34 * Inheritance:
35 * The DssMap class inherits from the Mapping class.
36
37 * Attributes:
38 * The DssMap class does not define any new attributes beyond those
39 * which are applicable to all Mappings.
40
41 * Functions:
42 c The DssMap class does not define any new functions beyond those
43 f The DssMap class does not define any new routines beyond those
44 * which are applicable to all Mappings.
45
46 * Copyright:
47 * Copyright (C) 1997-2006 Council for the Central Laboratory of the
48 * Research Councils
49
50 * Licence:
51 * This program is free software: you can redistribute it and/or
52 * modify it under the terms of the GNU Lesser General Public
53 * License as published by the Free Software Foundation, either
54 * version 3 of the License, or (at your option) any later
55 * version.
56 *
57 * This program is distributed in the hope that it will be useful,
58 * but WITHOUT ANY WARRANTY; without even the implied warranty of
59 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
60 * GNU Lesser General Public License for more details.
61 *
62 * You should have received a copy of the GNU Lesser General
63 * License along with this program. If not, see
64 * <http://www.gnu.org/licenses/>.
65 * (except for code supplied by Doug Mink, as noted in this file)
66
67 * Authors:
68 * DSB: D.S. Berry (Starlink)
69 * RFWS: R.F. Warren-Smith (Starlink, RAL)
70
71 * History:
72 * 18-FEB-1997 (DSB):
73 * Original version.
74 * 30-JUN-1997 (DSB):
75 * astDssFits and astDssMap made protected instead of public.
76 * 15-JUL-1997 (RFWS):
77 * Tidied public prologues.
78 * 5-SEP-197 (RFWS):
79 * Added prototypes for platepos and platepix.
80 * 4-NOV-1997 (DSB):
81 * o A copy of the supplied FitsChan is no longer stored inside
82 * the DssMap. The FitsChan returned by DssFits is now derived from
83 * the wcs information stored in the SAOimage "WorldCoor" structure
84 * (stored within the DssMap), and only contains the keywords
85 * necessary to reconstruct the DssMap.
86 * o The external representation of a DssMap is now stored in a set
87 * of scalar values, rather than a FitsChan.
88 * 22-DEC-1997 (DSB):
89 * Bug fixed in MapMerge which caused a core dump when a
90 * DssMap/WinMap combination is succesfully simplified.
91 * 8-JAN-2003 (DSB):
92 * Changed private InitVtab method to protected astInitDssMapVtab
93 * method.
94 * 14-FEB-2006 (DSB):
95 * Override astGetObjSize.
96 * 10-MAY-2006 (DSB):
97 * Override astEqual.
98 *class--
99 */
100
101 /* Module Macros. */
102 /* ============== */
103 /* Set the name of the class we are implementing. This indicates to
104 the header files that define class interfaces that they should make
105 "protected" symbols available. */
106 #define astCLASS DssMap
107
108 /* Macro which returns the nearest integer to a given floating point
109 value. */
110 #define NINT(x) (int)((x)+(((x)>0.0)?0.5:-0.5))
111
112 /* Macros which return the maximum and minimum of two values. */
113 #define MAX(aa,bb) ((aa)>(bb)?(aa):(bb))
114 #define MIN(aa,bb) ((aa)<(bb)?(aa):(bb))
115
116 /* Include files. */
117 /* ============== */
118 /* Interface definitions. */
119 /* ---------------------- */
120 #include "memory.h" /* Memory allocation facilities */
121
122 #include "globals.h" /* Thread-safe global data access */
123 #include "error.h" /* Error reporting facilities */
124 #include "object.h" /* Base Object class */
125 #include "pointset.h" /* Sets of points/coordinates */
126 #include "mapping.h" /* Coordinate mappings (parent class) */
127 #include "channel.h" /* I/O channels */
128 #include "fitschan.h" /* Manipulation of FITS header cards */
129 #include "wcsmap.h" /* Degrees/radians conversion factors */
130 #include "winmap.h" /* Shift and scale mappings */
131 #include "dssmap.h" /* Interface definition for this class */
132
133 /* Error code definitions. */
134 /* ----------------------- */
135 #include "ast_err.h" /* AST error codes */
136
137 /* C header files. */
138 /* --------------- */
139 #include <float.h>
140 #include <stdarg.h>
141 #include <stddef.h>
142 #include <stdio.h>
143 #include <string.h>
144 #include <math.h>
145
146 /* Module Variables. */
147 /* ================= */
148
149 /* Address of this static variable is used as a unique identifier for
150 member of this class. */
151 static int class_check;
152
153 /* Pointers to parent class methods which are extended by this class. */
154 static int (* parent_getobjsize)( AstObject *, int * );
155 static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
156
157
158 #ifdef THREAD_SAFE
159 /* Define how to initialise thread-specific globals. */
160 #define GLOBAL_inits \
161 globals->Class_Init = 0;
162
163 /* Create the function that initialises global data for this module. */
164 astMAKE_INITGLOBALS(DssMap)
165
166 /* Define macros for accessing each item of thread specific global data. */
167 #define class_init astGLOBAL(DssMap,Class_Init)
168 #define class_vtab astGLOBAL(DssMap,Class_Vtab)
169
170
171 #include <pthread.h>
172
173
174 #else
175
176
177 /* Define the class virtual function table and its initialisation flag
178 as static variables. */
179 static AstDssMapVtab class_vtab; /* Virtual function table */
180 static int class_init = 0; /* Virtual function table initialised? */
181
182 #endif
183
184 /* External Interface Function Prototypes. */
185 /* ======================================= */
186 /* The following functions have public prototypes only (i.e. no
187 protected prototypes), so we must provide local prototypes for use
188 within this module. */
189
190 /* Prototypes for Private Member Functions. */
191 /* ======================================== */
192 static AstFitsChan *DssFits( AstDssMap *, int * );
193 static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
194 static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * );
195 static int platepix( double, double, struct WorldCoor *, double *, double * );
196 static int platepos( double, double, struct WorldCoor *, double *, double * );
197 static struct WorldCoor *BuildWcs( AstFitsChan *, const char *, const char *, int * );
198 static void Copy( const AstObject *, AstObject *, int * );
199 static void Delete( AstObject *obj, int * );
200 static void Dump( AstObject *, AstChannel *, int * );
201 static int Equal( AstObject *, AstObject *, int * );
202
203 static int GetObjSize( AstObject *, int * );
204 /* Member functions. */
205 /* ================= */
Equal(AstObject * this_object,AstObject * that_object,int * status)206 static int Equal( AstObject *this_object, AstObject *that_object, int *status ) {
207 /*
208 * Name:
209 * Equal
210
211 * Purpose:
212 * Test if two DssMaps are equivalent.
213
214 * Type:
215 * Private function.
216
217 * Synopsis:
218 * #include "dssmap.h"
219 * int Equal( AstObject *this, AstObject *that, int *status )
220
221 * Class Membership:
222 * DssMap member function (over-rides the astEqual protected
223 * method inherited from the astMapping class).
224
225 * Description:
226 * This function returns a boolean result (0 or 1) to indicate whether
227 * two DssMaps are equivalent.
228
229 * Parameters:
230 * this
231 * Pointer to the first Object (a DssMap).
232 * that
233 * Pointer to the second Object.
234 * status
235 * Pointer to the inherited status variable.
236
237 * Returned Value:
238 * One if the DssMaps are equivalent, zero otherwise.
239
240 * Notes:
241 * - A value of zero will be returned if this function is invoked
242 * with the global status set, or if it should fail for any reason.
243 */
244
245 /* Local Variables: */
246 AstDssMap *that;
247 AstDssMap *this;
248 int i;
249 int nin;
250 int nout;
251 int result;
252 struct WorldCoor *this_wcs;
253 struct WorldCoor *that_wcs;
254
255 /* Initialise. */
256 result = 0;
257
258 /* Check the global error status. */
259 if ( !astOK ) return result;
260
261 /* Obtain pointers to the two DssMap structures. */
262 this = (AstDssMap *) this_object;
263 that = (AstDssMap *) that_object;
264
265 /* Check the second object is a DssMap. We know the first is a
266 DssMap since we have arrived at this implementation of the virtual
267 function. */
268 if( astIsADssMap( that ) ) {
269
270 /* Get the number of inputs and outputs and check they are the same for both. */
271 nin = astGetNin( this );
272 nout = astGetNout( this );
273 if( astGetNin( that ) == nin && astGetNout( that ) == nout ) {
274
275 /* If the Invert flags for the two DssMaps differ, it may still be possible
276 for them to be equivalent. First compare the DssMaps if their Invert
277 flags are the same. In this case all the attributes of the two DssMaps
278 must be identical. */
279 if( astGetInvert( this ) == astGetInvert( that ) ) {
280
281 this_wcs = ( struct WorldCoor *) this->wcs;
282 that_wcs = ( struct WorldCoor *) that->wcs;
283
284 if( this_wcs->x_pixel_offset == that_wcs->x_pixel_offset &&
285 this_wcs->y_pixel_offset == that_wcs->y_pixel_offset &&
286 this_wcs->ppo_coeff[2] == that_wcs->ppo_coeff[2] &&
287 this_wcs->ppo_coeff[5] == that_wcs->ppo_coeff[5] &&
288 this_wcs->x_pixel_size == that_wcs->x_pixel_size &&
289 this_wcs->y_pixel_size == that_wcs->y_pixel_size &&
290 this_wcs->plate_dec == that_wcs->plate_dec &&
291 this_wcs->plate_ra == that_wcs->plate_ra ) {
292
293 result = 1;
294 for( i = 0; i < 13; i++ ) {
295 if( this_wcs->amd_x_coeff[i] != that_wcs->amd_x_coeff[i] ||
296 this_wcs->amd_y_coeff[i] != that_wcs->amd_y_coeff[i] ) {
297 result = 0;
298 break;
299 }
300 }
301
302 }
303
304 /* If the Invert flags for the two DssMaps differ, the attributes of the two
305 DssMaps must be inversely related to each other. */
306 } else {
307
308 /* In the specific case of a DssMap, Invert flags must be equal. */
309 result = 0;
310
311 }
312 }
313 }
314
315 /* If an error occurred, clear the result value. */
316 if ( !astOK ) result = 0;
317
318 /* Return the result, */
319 return result;
320 }
321
GetObjSize(AstObject * this_object,int * status)322 static int GetObjSize( AstObject *this_object, int *status ) {
323 /*
324 * Name:
325 * GetObjSize
326
327 * Purpose:
328 * Return the in-memory size of an Object.
329
330 * Type:
331 * Private function.
332
333 * Synopsis:
334 * #include "dssmap.h"
335 * int GetObjSize( AstObject *this, int *status )
336
337 * Class Membership:
338 * DssMap member function (over-rides the astGetObjSize protected
339 * method inherited from the parent class).
340
341 * Description:
342 * This function returns the in-memory size of the supplied DssMap,
343 * in bytes.
344
345 * Parameters:
346 * this
347 * Pointer to the DssMap.
348 * status
349 * Pointer to the inherited status variable.
350
351 * Returned Value:
352 * The Object size, in bytes.
353
354 * Notes:
355 * - A value of zero will be returned if this function is invoked
356 * with the global status set, or if it should fail for any reason.
357 */
358
359 /* Local Variables: */
360 AstDssMap *this; /* Pointer to DssMap structure */
361 int result; /* Result value to return */
362
363 /* Initialise. */
364 result = 0;
365
366 /* Check the global error status. */
367 if ( !astOK ) return result;
368
369 /* Obtain a pointers to the DssMap structure. */
370 this = (AstDssMap *) this_object;
371
372 /* Invoke the GetObjSize method inherited from the parent class, and then
373 add on any components of the class structure defined by thsi class
374 which are stored in dynamically allocated memory. */
375 result = (*parent_getobjsize)( this_object, status );
376 result += astTSizeOf( this->wcs );
377
378 /* If an error occurred, clear the result value. */
379 if ( !astOK ) result = 0;
380
381 /* Return the result, */
382 return result;
383 }
384
385
BuildWcs(AstFitsChan * fits,const char * method,const char * class,int * status)386 static struct WorldCoor *BuildWcs( AstFitsChan *fits, const char *method,
387 const char *class, int *status ) {
388 /*
389 * Name:
390 * BuildWcs
391
392 * Purpose:
393 * Copy DSS plate fit information from a FitsChan to a SAOimage
394 * WorldCoor structure.
395
396 * Type:
397 * Private function.
398
399 * Synopsis:
400 * #include "dssmap.h"
401 * struct WorldCoor *BuildWcs( AstFitsChan *fits, const char *method,
402 * const char *class )
403
404 * Class Membership:
405 * DssMap member function.
406
407 * Description:
408 * This creates a WorldCoor structure and copies the required data
409 * from the supplied FitsChan into the new WorldCoor structure. Note,
410 * only those components of the WorldCoor structure which are needed to
411 * transform between pixel and sky coordinates are initialised in the
412 * returned structure.
413
414 * Parameters:
415 * fits
416 * Pointer to the FitsChan containing the FITS header describing
417 * the DSS plate fit.
418 * method
419 * The calling method (for error messages).
420 * class
421 * The object class (for error messages).
422
423 * Returned Value:
424 * A pointer to the new WorldCoor structure. This should be freed
425 * using astFree when no longer needed.
426
427 * Notes:
428 * - A NULL pointer is returned if an error has already occurred, or
429 * if this function should fail for any reason.
430
431 */
432
433 /* Local Variables: */
434 char name_buff[ 10 ]; /* Buffer for keyword name */
435 char *name; /* Pointer to jeyword name string */
436 char *ckeyval; /* Pointer to string keyword value */
437 struct WorldCoor *ret; /* Pointer to the returned structure */
438 double rah,ram,ras; /* Centre RA hours, minutes and seconds */
439 double dsign; /* Sign of centre dec */
440 double decd,decm,decs; /* Centre Dec degrees, minutes, seconds */
441 double dec_deg; /* Centre Dec in degrees */
442 double ra_hours; /* Centre RA in hours */
443 int i; /* Coefficient index */
444
445 /* Check the local error status. */
446 if ( !astOK ) return NULL;
447
448 /* Get memory to hold the returned structure. */
449 ret = (struct WorldCoor *) astMalloc( sizeof( struct WorldCoor ) );
450
451 /* Check the memory can be used. */
452 if( astOK ){
453
454 /* The following code is based on the "wcsinit" function in SAOimage file
455 wcs.c. Note, only the values needed in the platepos and platepix
456 functions are set up. The FITS keywords are accessed in the order in
457 which they are usually stored in a FITS file. This will cut down the
458 time spent searching for keywords. Report an error if any required
459 keyword is not found. */
460
461 /* Plate center RA */
462 rah = 0.0;
463 ram = 0.0;
464 ras = 0.0;
465
466 name = "PLTRAH";
467 if( !astGetFitsF( fits, name, &rah ) && astOK ){
468 astError( AST__BDFTS, "%s(%s): No value has been supplied for the "
469 "FITS keyword '%s'.", status, method, class, name );
470 }
471
472 name = "PLTRAM";
473 if( !astGetFitsF( fits, name, &ram ) && astOK ){
474 astError( AST__BDFTS, "%s(%s): No value has been supplied for the "
475 "FITS keyword '%s'.", status, method, class, name );
476 }
477
478 name = "PLTRAS";
479 if( !astGetFitsF( fits, name, &ras ) && astOK ){
480 astError( AST__BDFTS, "%s(%s): No value has been supplied for the "
481 "FITS keyword '%s'.", status, method, class, name );
482 }
483
484 ra_hours = rah + (ram / (double)60.0) + (ras / (double)3600.0);
485 ret->plate_ra = AST__DD2R*15.0*ra_hours;
486
487
488 /* Plate center Dec */
489 name = "PLTDECSN";
490 if( !astGetFitsS( fits, name, &ckeyval ) && astOK ){
491 dsign = 1.0;
492
493 } else {
494 if( *ckeyval == '-' ){
495 dsign = -1.0;
496 } else {
497 dsign = 1.0;
498 }
499
500 }
501
502 decd = 0.0;
503 decm = 0.0;
504 decs = 0.0;
505
506 name = "PLTDECD";
507 if( !astGetFitsF( fits, name, &decd ) && astOK ){
508 astError( AST__BDFTS, "%s(%s): No value has been supplied for the "
509 "FITS keyword '%s'.", status, method, class, name );
510 }
511
512 name = "PLTDECM";
513 if( !astGetFitsF( fits, name, &decm ) && astOK ){
514 astError( AST__BDFTS, "%s(%s): No value has been supplied for the "
515 "FITS keyword '%s'.", status, method, class, name );
516 }
517
518 name = "PLTDECS";
519 if( !astGetFitsF( fits, name, &decs ) && astOK ){
520 astError( AST__BDFTS, "%s(%s): No value has been supplied for the "
521 "FITS keyword '%s'.", status, method, class, name );
522 }
523
524 dec_deg = dsign * (decd+(decm/(double)60.0)+(decs/(double)3600.0));
525 ret->plate_dec = AST__DD2R*dec_deg;
526
527 /* Plate Scale arcsec per mm */
528 name = "PLTSCALE";
529 if( !astGetFitsF( fits, name, &ret->plate_scale ) && astOK ){
530 astError( AST__BDFTS, "%s(%s): No value has been supplied for the "
531 "FITS keyword '%s'.", status, method, class, name );
532 }
533
534 /* X and Y corners (in pixels) */
535 name = "CNPIX1";
536 if( !astGetFitsF( fits, name, &ret->x_pixel_offset ) && astOK ){
537 astError( AST__BDFTS, "%s(%s): No value has been supplied for the "
538 "FITS keyword '%s'.", status, method, class, name );
539 }
540
541 name = "CNPIX2";
542 if( !astGetFitsF( fits, name, &ret->y_pixel_offset ) && astOK ){
543 astError( AST__BDFTS, "%s(%s): No value has been supplied for the "
544 "FITS keyword '%s'.", status, method, class, name );
545 }
546
547 /* X and Y pixel sizes (microns). */
548 name = "XPIXELSZ";
549 if( !astGetFitsF( fits, name, &ret->x_pixel_size ) && astOK ){
550 astError( AST__BDFTS, "%s(%s): No value has been supplied for the "
551 "FITS keyword '%s'.", status, method, class, name );
552 }
553
554 name = "YPIXELSZ";
555 if( !astGetFitsF( fits, name, &ret->y_pixel_size ) && astOK ){
556 astError( AST__BDFTS, "%s(%s): No value has been supplied for the "
557 "FITS keyword '%s'.", status, method, class, name );
558 }
559
560 /* Orientation Coefficients. Only report an error if PPO3 or PPO6 are
561 missing (these are the only two which are actually used). Assume a
562 value of zero for any of the others which are missing. */
563 name = name_buff;
564 for ( i = 0; i < 6; i++ ) {
565 sprintf( name_buff, "PPO%d", i + 1 );
566 if( !astGetFitsF( fits, name, &ret->ppo_coeff[i] ) ) {
567 ret->ppo_coeff[i] = 0.0;
568 if( ( i == 2 || i == 5 ) && astOK ) {
569 astError( AST__BDFTS, "%s(%s): No value has been supplied "
570 "for the FITS keyword '%s'.", status, method, class,
571 name );
572 break;
573 }
574 }
575 }
576
577 /* Plate solution x and y coefficients. Report an error if any of
578 coefficients 1 to 14 are missing. Assume a value of zero for any
579 others which are missing. */
580 name = name_buff;
581 for( i = 0; i < 19; i++ ){
582 sprintf( name_buff, "AMDX%d", i + 1 );
583 if( !astGetFitsF( fits, name, &ret->amd_x_coeff[i] ) ) {
584 ret->amd_x_coeff[i] = 0.0;
585 if( i < 13 && astOK ){
586 astError( AST__BDFTS, "%s(%s): No value has been supplied "
587 "for the FITS keyword '%s'.", status, method, class, name );
588 break;
589 }
590 }
591 }
592
593 for( i = 0; i < 19; i++ ){
594 sprintf( name_buff, "AMDY%d", i + 1 );
595 if( !astGetFitsF( fits, name, &ret->amd_y_coeff[i] ) ){
596 ret->amd_y_coeff[i] = 0.0;
597 if( i < 13 && astOK ){
598 astError( AST__BDFTS, "%s(%s): No value has been supplied "
599 "for the FITS keyword '%s'.", status, method, class, name );
600 break;
601 }
602 }
603 }
604
605 /* If anything went wrong, free the returned structure. */
606 if( !astOK ) ret = (struct WorldCoor *) astFree( (void *) ret );
607 }
608
609 /* Return the pointer. */
610 return ret;
611 }
612
DssFits(AstDssMap * this,int * status)613 static AstFitsChan *DssFits( AstDssMap *this, int *status ) {
614 /*
615 *+
616 * Name:
617 * astDssFits
618
619 * Purpose:
620 * Return a pointer to a FitsChan describing a DssMap.
621
622 * Type:
623 * Protected virtual function.
624
625 * Synopsis:
626 * #include "dssmap.h"
627 * AstFitsChan *DssFits( AstDssMap *this )
628
629 * Class Membership:
630 * DssMap method.
631
632 * Description:
633 * This function returns a pointer to a DSS-encoded FitsChan containing
634 * cards generated from the information stored with the DssMap. The
635 * keywords contained in the FitsChan are those which would ne needed to
636 * re-create the DssMap (see astDSSMap).
637
638 * Parameters:
639 * this
640 * Pointer to the DssMap.
641
642 * Returned Value:
643 * astDssFits()
644 * A pointer to the FitsChan.
645
646 * Notes:
647 * - The returned pointer should be annuled using astAnnul when no longer
648 * needed.
649 * - A value of NULL will be returned if this function is invoked
650 * with the global error status set, or if it should fail for any
651 * reason.
652 *-
653 */
654
655 /* Local Variables: */
656 AstFitsChan *ret; /* Pointer to the returned FitsChan */
657 char *comm; /* Pointer to keyword comment string */
658 char *name; /* Pointer to keyword name string */
659 char name_buff[ 10 ]; /* Buffer for keyword name */
660 double dec; /* Centre Dec in degrees */
661 double decd,decm,decs; /* Centre Dec degrees, minutes, seconds */
662 double ra; /* Centre RA in hours */
663 double rah,ram,ras; /* Centre RA hours, minutes and seconds */
664 int i; /* Coefficient index */
665 struct WorldCoor *wcs; /* WCS information from the DssMap */
666
667 /* Check the global error status. */
668 if ( !astOK ) return NULL;
669
670 /* Store a pointer to the WCS information stored in the DSSMap. */
671 wcs = (struct WorldCoor *) this->wcs,
672
673 /* Create a new empty FitsChan, using DSS encoding. */
674 ret = astFitsChan( NULL, NULL, "Encoding=DSS", status );
675
676 /* Create the keyword values and stored them in the returned FitsChan... */
677
678 /* Plate centre RA. */
679 ra = wcs->plate_ra/( AST__DD2R*15.0 );
680 ra = modf( ra, &rah );
681 ra = modf( 60.0*ra, &ram );
682 ras = 60.0*ra;
683
684 astSetFitsI( ret, "PLTRAH", NINT( rah ), "Plate centre RA", 0 );
685 astSetFitsI( ret, "PLTRAM", NINT( ram ), " ", 0 );
686 astSetFitsF( ret, "PLTRAS", ras, " ", 0 );
687
688 /* Plate centre DEC. */
689 dec = wcs->plate_dec/AST__DD2R;
690 if( dec < 0.0 ) {
691 dec = -dec;
692 astSetFitsS( ret, "PLTDECSN", "-", "Plate centre DEC", 0 );
693 } else {
694 astSetFitsS( ret, "PLTDECSN", "+", "Plate centre DEC", 0 );
695 }
696
697 dec = modf( dec, &decd );
698 dec = modf( 60.0*dec, &decm );
699 decs = 60.0*dec;
700
701 astSetFitsI( ret, "PLTDECD", NINT( decd ), " ", 0 );
702 astSetFitsI( ret, "PLTDECM", NINT( decm ), " ", 0 );
703 astSetFitsF( ret, "PLTDECS", decs, " ", 0 );
704
705 /* Plate Scale arcsec per mm */
706 astSetFitsF( ret, "PLTSCALE", wcs->plate_scale, "Plate Scale arcsec per mm",
707 0 );
708
709 /* X and Y corners (in pixels) */
710 astSetFitsI( ret, "CNPIX1", NINT( wcs->x_pixel_offset ),
711 "X corner (pixels)", 0 );
712 astSetFitsI( ret, "CNPIX2", NINT( wcs->y_pixel_offset ),
713 "Y corner", 0 );
714
715 /* X and Y pixel sizes (microns). */
716 astSetFitsF( ret, "XPIXELSZ", wcs->x_pixel_size,
717 "X pixel size (microns)", 0 );
718 astSetFitsF( ret, "YPIXELSZ", wcs->y_pixel_size,
719 "Y pixel size (microns)", 0 );
720
721 /* Orientation Coefficients. */
722 name = name_buff;
723 comm = "Orientation Coefficients";
724 for ( i = 0; i < 6; i++ ) {
725 sprintf( name_buff, "PPO%d", i + 1 );
726 astSetFitsF( ret, name, wcs->ppo_coeff[i], comm, 0 );
727 comm = " ";
728 }
729
730 /* Plate solution x and y coefficients. */
731 comm = "Plate solution x coefficients";
732 for( i = 0; i < 19; i++ ){
733 sprintf( name_buff, "AMDX%d", i + 1 );
734 astSetFitsF( ret, name, wcs->amd_x_coeff[i], comm, 0 );
735 comm = " ";
736 }
737
738 comm = "Plate solution y coefficients";
739 for( i = 0; i < 19; i++ ){
740 sprintf( name_buff, "AMDY%d", i + 1 );
741 astSetFitsF( ret, name, wcs->amd_y_coeff[i], comm, 0 );
742 comm = " ";
743 }
744
745 /* Return a pointer to the FitsChan. */
746 return ret;
747 }
748
astInitDssMapVtab_(AstDssMapVtab * vtab,const char * name,int * status)749 void astInitDssMapVtab_( AstDssMapVtab *vtab, const char *name, int *status ) {
750 /*
751 *+
752 * Name:
753 * astInitDssMapVtab
754
755 * Purpose:
756 * Initialise a virtual function table for a DssMap.
757
758 * Type:
759 * Protected function.
760
761 * Synopsis:
762 * #include "dssmap.h"
763 * void astInitDssMapVtab( AstDssMapVtab *vtab, const char *name )
764
765 * Class Membership:
766 * DssMap vtab initialiser.
767
768 * Description:
769 * This function initialises the component of a virtual function
770 * table which is used by the DssMap class.
771
772 * Parameters:
773 * vtab
774 * Pointer to the virtual function table. The components used by
775 * all ancestral classes will be initialised if they have not already
776 * been initialised.
777 * name
778 * Pointer to a constant null-terminated character string which contains
779 * the name of the class to which the virtual function table belongs (it
780 * is this pointer value that will subsequently be returned by the Object
781 * astClass function).
782 *-
783 */
784
785 /* Local Variables: */
786 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
787 AstObjectVtab *object; /* Pointer to Object component of Vtab */
788 AstMappingVtab *mapping; /* Pointer to Mapping component of Vtab */
789
790 /* Check the local error status. */
791 if ( !astOK ) return;
792
793 /* Get a pointer to the thread specific global data structure. */
794 astGET_GLOBALS(NULL);
795
796 /* Initialize the component of the virtual function table used by the
797 parent class. */
798 astInitMappingVtab( (AstMappingVtab *) vtab, name );
799
800 /* Store a unique "magic" value in the virtual function table. This
801 will be used (by astIsADssMap) to determine if an object belongs
802 to this class. We can conveniently use the address of the (static)
803 class_check variable to generate this unique value. */
804 vtab->id.check = &class_check;
805 vtab->id.parent = &(((AstMappingVtab *) vtab)->id);
806
807 /* Initialise member function pointers. */
808 /* ------------------------------------ */
809 /* Store pointers to the member functions (implemented here) that provide
810 virtual methods for this class. */
811 vtab->DssFits = DssFits;
812
813 /* Save the inherited pointers to methods that will be extended, and
814 replace them with pointers to the new member functions. */
815 mapping = (AstMappingVtab *) vtab;
816 object = (AstObjectVtab *) vtab;
817
818 parent_transform = mapping->Transform;
819 parent_getobjsize = object->GetObjSize;
820 object->GetObjSize = GetObjSize;
821 mapping->Transform = Transform;
822
823 /* Store replacement pointers for methods which will be over-ridden by
824 new member functions implemented here. */
825 object->Equal = Equal;
826 mapping->MapMerge = MapMerge;
827
828 /* Declare the class dump, copy and delete function. */
829 astSetDump( object, Dump, "DssMap", "DSS plate fit mapping" );
830 astSetCopy( object, Copy );
831 astSetDelete( object, Delete );
832
833 /* If we have just initialised the vtab for the current class, indicate
834 that the vtab is now initialised, and store a pointer to the class
835 identifier in the base "object" level of the vtab. */
836 if( vtab == &class_vtab ) {
837 class_init = 1;
838 astSetVtabClassIdentifier( vtab, &(vtab->id) );
839 }
840 }
841
MapMerge(AstMapping * this,int where,int series,int * nmap,AstMapping *** map_list,int ** invert_list,int * status)842 static int MapMerge( AstMapping *this, int where, int series, int *nmap,
843 AstMapping ***map_list, int **invert_list, int *status ) {
844 /*
845 * Name:
846 * MapMerge
847
848 * Purpose:
849 * Simplify a sequence of Mappings containing a DssMap.
850
851 * Type:
852 * Private function.
853
854 * Synopsis:
855 * #include "mapping.h"
856 * int MapMerge( AstMapping *this, int where, int series, int *nmap,
857 * AstMapping ***map_list, int **invert_list, int *status )
858
859 * Class Membership:
860 * DssMap method (over-rides the protected astMapMerge method
861 * inherited from the Mapping class).
862
863 * Description:
864 * This function attempts to simplify a sequence of Mappings by
865 * merging a nominated DssMap in the sequence with its neighbours,
866 * so as to shorten the sequence if possible.
867 *
868 * In many cases, simplification will not be possible and the
869 * function will return -1 to indicate this, without further
870 * action.
871 *
872 * In most cases of interest, however, this function will either
873 * attempt to replace the nominated DssMap with a Mapping which it
874 * considers simpler, or to merge it with the Mappings which
875 * immediately precede it or follow it in the sequence (both will
876 * normally be considered). This is sufficient to ensure the
877 * eventual simplification of most Mapping sequences by repeated
878 * application of this function.
879 *
880 * In some cases, the function may attempt more elaborate
881 * simplification, involving any number of other Mappings in the
882 * sequence. It is not restricted in the type or scope of
883 * simplification it may perform, but will normally only attempt
884 * elaborate simplification in cases where a more straightforward
885 * approach is not adequate.
886
887 * Parameters:
888 * this
889 * Pointer to the nominated DssMap which is to be merged with
890 * its neighbours. This should be a cloned copy of the DssMap
891 * pointer contained in the array element "(*map_list)[where]"
892 * (see below). This pointer will not be annulled, and the
893 * DssMap it identifies will not be modified by this function.
894 * where
895 * Index in the "*map_list" array (below) at which the pointer
896 * to the nominated DssMap resides.
897 * series
898 * A non-zero value indicates that the sequence of Mappings to
899 * be simplified will be applied in series (i.e. one after the
900 * other), whereas a zero value indicates that they will be
901 * applied in parallel (i.e. on successive sub-sets of the
902 * input/output coordinates).
903 * nmap
904 * Address of an int which counts the number of Mappings in the
905 * sequence. On entry this should be set to the initial number
906 * of Mappings. On exit it will be updated to record the number
907 * of Mappings remaining after simplification.
908 * map_list
909 * Address of a pointer to a dynamically allocated array of
910 * Mapping pointers (produced, for example, by the astMapList
911 * method) which identifies the sequence of Mappings. On entry,
912 * the initial sequence of Mappings to be simplified should be
913 * supplied.
914 *
915 * On exit, the contents of this array will be modified to
916 * reflect any simplification carried out. Any form of
917 * simplification may be performed. This may involve any of: (a)
918 * removing Mappings by annulling any of the pointers supplied,
919 * (b) replacing them with pointers to new Mappings, (c)
920 * inserting additional Mappings and (d) changing their order.
921 *
922 * The intention is to reduce the number of Mappings in the
923 * sequence, if possible, and any reduction will be reflected in
924 * the value of "*nmap" returned. However, simplifications which
925 * do not reduce the length of the sequence (but improve its
926 * execution time, for example) may also be performed, and the
927 * sequence might conceivably increase in length (but normally
928 * only in order to split up a Mapping into pieces that can be
929 * more easily merged with their neighbours on subsequent
930 * invocations of this function).
931 *
932 * If Mappings are removed from the sequence, any gaps that
933 * remain will be closed up, by moving subsequent Mapping
934 * pointers along in the array, so that vacated elements occur
935 * at the end. If the sequence increases in length, the array
936 * will be extended (and its pointer updated) if necessary to
937 * accommodate any new elements.
938 *
939 * Note that any (or all) of the Mapping pointers supplied in
940 * this array may be annulled by this function, but the Mappings
941 * to which they refer are not modified in any way (although
942 * they may, of course, be deleted if the annulled pointer is
943 * the final one).
944 * invert_list
945 * Address of a pointer to a dynamically allocated array which,
946 * on entry, should contain values to be assigned to the Invert
947 * attributes of the Mappings identified in the "*map_list"
948 * array before they are applied (this array might have been
949 * produced, for example, by the astMapList method). These
950 * values will be used by this function instead of the actual
951 * Invert attributes of the Mappings supplied, which are
952 * ignored.
953 *
954 * On exit, the contents of this array will be updated to
955 * correspond with the possibly modified contents of the
956 * "*map_list" array. If the Mapping sequence increases in
957 * length, the "*invert_list" array will be extended (and its
958 * pointer updated) if necessary to accommodate any new
959 * elements.
960 * status
961 * Pointer to the inherited status variable.
962
963 * Returned Value:
964 * If simplification was possible, the function returns the index
965 * in the "map_list" array of the first element which was
966 * modified. Otherwise, it returns -1 (and makes no changes to the
967 * arrays supplied).
968
969 * Notes:
970 * - A value of -1 will be returned if this function is invoked
971 * with the global error status set, or if it should fail for any
972 * reason.
973 */
974
975 /* Local Variables: */
976 AstDssMap *dm; /* Pointer to supplied DssMap */
977 AstDssMap *dmnew; /* Pointer to replacement DssMap */
978 AstFitsChan *fits; /* FITS headers for replacement DssMap */
979 AstFitsChan *fits_dss;/* FITS headers for supplied DssMap */
980 AstWinMap *wm; /* Pointer to the adjacent WinMap */
981 double *a; /* Pointer to shift terms */
982 double *b; /* Pointer to scale terms */
983 double cnpix1; /* X pixel origin */
984 double cnpix2; /* Y pixel origin */
985 double xpixelsz; /* X pixel size */
986 double ypixelsz; /* Y pixel size */
987 int i; /* Loop counter */
988 int ok; /* All FITS keywords found? */
989 int old_winv; /* original Invert value for supplied WinMap */
990 int result; /* Result value to return */
991 int wmi; /* Index of adjacent WinMap in map list */
992 struct WorldCoor *wcs;/* Pointer to SAOimage wcs structure */
993
994 /* Initialise. */
995 result = -1;
996
997 /* Check the global error status. */
998 if ( !astOK ) return result;
999
1000 /* The only simplification easily possible is if a WinMap maps the pixel
1001 coordinates prior to a DssMap. If the DssMap has not been inverted, the
1002 WinMap must be applied before the DssMap, otherwise the WinMap must be
1003 applied after the DssMap. */
1004 if( series ){
1005
1006 if( !( *invert_list )[ where ] ){
1007 wmi = where - 1;
1008 } else {
1009 wmi = where + 1;
1010 }
1011
1012 if( wmi >= 0 && wmi < *nmap ){
1013 if( !strcmp( astGetClass( ( *map_list )[ wmi ] ), "WinMap" ) ){
1014
1015 /* Temporarily set the Invert attribute of the WinMap to the supplied value. */
1016 wm = (AstWinMap *) ( *map_list )[ wmi ];
1017 old_winv = astGetInvert( wm );
1018 astSetInvert( wm, ( *invert_list )[ wmi ] );
1019
1020 /* Get a copy of the scale and shift terms from the WinMap. */
1021 astWinTerms( wm, &a, &b );
1022
1023 /* Check that the scale and shift terms are usable. */
1024 if( astOK &&
1025 a[ 0 ] != AST__BAD && b[ 0 ] != AST__BAD && b[ 0 ] != 0.0 &&
1026 a[ 1 ] != AST__BAD && b[ 1 ] != AST__BAD && b[ 1 ] != 0.0 ){
1027
1028 /* Get the values of the keywords which define the origin and scales of
1029 the DssMap pixel coordinate system. */
1030 dm = (AstDssMap *) ( *map_list )[ where ];
1031 wcs = (struct WorldCoor *) ( dm->wcs );
1032
1033 cnpix1 = wcs->x_pixel_offset;
1034 cnpix2 = wcs->y_pixel_offset;
1035 xpixelsz = wcs->x_pixel_size;
1036 ypixelsz = wcs->y_pixel_size;
1037
1038 /* Calculate new values which take into account the WinMap. */
1039 if( wmi == where - 1 ){
1040 xpixelsz *= b[ 0 ];
1041 ypixelsz *= b[ 1 ];
1042 cnpix1 = 0.5 + ( cnpix1 + a[ 0 ] - 0.5 )/b[ 0 ];
1043 cnpix2 = 0.5 + ( cnpix2 + a[ 1 ] - 0.5 )/b[ 1 ];
1044
1045 } else {
1046 xpixelsz /= b[ 0 ];
1047 ypixelsz /= b[ 1 ];
1048 cnpix1 = b[ 0 ]*( cnpix1 - 0.5 ) - a[ 0 ] + 0.5;
1049 cnpix2 = b[ 1 ]*( cnpix2 - 0.5 ) - a[ 1 ] + 0.5;
1050 }
1051
1052 /* The CNPIX1 and CNPIX2 keywords are integer keywords. Therefore, we can
1053 only do the simplification if the new values are integer to a good
1054 approximation. We use one hundredth of a pixel. */
1055 if( fabs( cnpix1 - NINT( cnpix1 ) ) < 0.01 &&
1056 fabs( cnpix2 - NINT( cnpix2 ) ) < 0.01 ){
1057
1058 /* Get a copy of the FitsChan holding the header cards which define the
1059 DssMap. */
1060 fits_dss = astDssFits( dm );
1061 fits = astCopy( fits_dss );
1062 fits_dss = astAnnul( fits_dss );
1063
1064 /* Update the value of each of the changed keywords. */
1065 ok = 1;
1066
1067 astClearCard( fits );
1068 if( astFindFits( fits, "CNPIX1", NULL, 0 ) ){
1069 astSetFitsI( fits, "CNPIX1", NINT( cnpix1 ), NULL, 1 );
1070 } else {
1071 ok = 0;
1072 }
1073
1074 astClearCard( fits );
1075 if( astFindFits( fits, "CNPIX2", NULL, 0 ) ){
1076 astSetFitsI( fits, "CNPIX2", NINT( cnpix2 ), NULL, 1 );
1077 } else {
1078 ok = 0;
1079 }
1080
1081 astClearCard( fits );
1082 if( astFindFits( fits, "XPIXELSZ", NULL, 0 ) ){
1083 astSetFitsF( fits, "XPIXELSZ", xpixelsz, NULL, 1 );
1084 } else {
1085 ok = 0;
1086 }
1087
1088 astClearCard( fits );
1089 if( astFindFits( fits, "YPIXELSZ", NULL, 0 ) ){
1090 astSetFitsF( fits, "YPIXELSZ", ypixelsz, NULL, 1 );
1091 } else {
1092 ok = 0;
1093 }
1094
1095 /* If all the keywords were updated succesfully, create the new DssMap
1096 based on the modified FITS header cards. */
1097 if( ok ){
1098 dmnew = astDssMap( fits, "", status );
1099
1100 /* Anull the DssMap pointer in the list and replace it with the new one.
1101 The invert flag is left unchanged. */
1102 dm = astAnnul( dm );
1103 ( *map_list )[ where ] = (AstMapping *) dmnew;
1104
1105 /* Annul the WinMap pointer in the list, and shuffle any remaining
1106 Mappings down to fill the gap. */
1107 wm = astAnnul( wm );
1108 for ( i = wmi + 1; i < *nmap; i++ ) {
1109 ( *map_list )[ i - 1 ] = ( *map_list )[ i ];
1110 ( *invert_list )[ i - 1 ] = ( *invert_list )[ i ];
1111 }
1112
1113 /* Clear the vacated element at the end. */
1114 ( *map_list )[ *nmap - 1 ] = NULL;
1115 ( *invert_list )[ *nmap - 1 ] = 0;
1116
1117 /* Decrement the Mapping count and return the index of the first
1118 modified element. */
1119 ( *nmap )--;
1120 result = MIN( wmi, where );
1121
1122 }
1123
1124 /* Annul the FitsChan holding the modified header cards. */
1125 fits = astAnnul( fits );
1126 }
1127 }
1128
1129 /* Free the arrays holding scale and shift terms from the WinMap. */
1130 a = (double *) astFree( (void *) a );
1131 b = (double *) astFree( (void *) b );
1132
1133 /* Reinstate the original setting of the Invert attribute of the WinMap (if
1134 it still exists). */
1135 if( wm ) astSetInvert( wm, old_winv );
1136
1137 }
1138 }
1139 }
1140
1141 /* Return the result. */
1142 return result;
1143 }
1144
Transform(AstMapping * this,AstPointSet * in,int forward,AstPointSet * out,int * status)1145 static AstPointSet *Transform( AstMapping *this, AstPointSet *in,
1146 int forward, AstPointSet *out, int *status ) {
1147 /*
1148 * Name:
1149 * Transform
1150
1151 * Purpose:
1152 * Apply a DssMap to transform a set of points.
1153
1154 * Type:
1155 * Private function.
1156
1157 * Synopsis:
1158 * #include "dssmap.h"
1159 * AstPointSet *Transform( AstMapping *this, AstPointSet *in,
1160 * int forward, AstPointSet *out, int *status )
1161
1162 * Class Membership:
1163 * DssMap member function (over-rides the astTransform protected
1164 * method inherited from the Mapping class).
1165
1166 * Description:
1167 * This function takes a DssMap and a set of points encapsulated in a
1168 * PointSet and transforms the points so as to apply the required DSS
1169 * plate fit.
1170
1171 * Parameters:
1172 * this
1173 * Pointer to the DssMap.
1174 * in
1175 * Pointer to the PointSet holding the input coordinate data.
1176 * forward
1177 * A non-zero value indicates that the forward coordinate transformation
1178 * should be applied, while a zero value requests the inverse
1179 * transformation.
1180 * out
1181 * Pointer to a PointSet which will hold the transformed (output)
1182 * coordinate values. A NULL value may also be given, in which case a
1183 * new PointSet will be created by this function.
1184 * status
1185 * Pointer to the inherited status variable.
1186
1187 * Returned Value:
1188 * Pointer to the output (possibly new) PointSet.
1189
1190 * Notes:
1191 * - A null pointer will be returned if this function is invoked with the
1192 * global error status set, or if it should fail for any reason.
1193 * - The number of coordinate values per point in the input PointSet must
1194 * be 2.
1195 * - If an output PointSet is supplied, it must have space for sufficient
1196 * number of points and coordinate values per point to accommodate the
1197 * result. Any excess space will be ignored.
1198 */
1199
1200 /* Local Variables: */
1201 AstPointSet *result; /* Pointer to output PointSet */
1202 AstDssMap *map; /* Pointer to DssMap to be applied */
1203 double **ptr_in; /* Pointer to input coordinate data */
1204 double **ptr_out; /* Pointer to output coordinate data */
1205 double *aa; /* Pointer to next longitude value */
1206 double *bb; /* Pointer to next latitude value */
1207 double *xx; /* Pointer to next pixel X value */
1208 double *yy; /* Pointer to next pixel Y value */
1209 int npoint; /* Number of points */
1210 int point; /* Loop counter for points */
1211
1212 /* Check the global error status. */
1213 if ( !astOK ) return NULL;
1214
1215 /* Obtain a pointer to the DssMap. */
1216 map = (AstDssMap *) this;
1217
1218 /* Apply the parent mapping using the stored pointer to the Transform member
1219 function inherited from the parent Mapping class. This function validates
1220 all arguments and generates an output PointSet if necessary, but does not
1221 actually transform any coordinate values. */
1222 result = (*parent_transform)( this, in, forward, out, status );
1223
1224 /* We will now extend the parent astTransform method by performing the
1225 calculations needed to generate the output coordinate values. */
1226
1227 /* Determine the numbers of points from the input PointSet and obtain
1228 pointers for accessing the input and output coordinate values. */
1229 npoint = astGetNpoint( in );
1230 ptr_in = astGetPoints( in );
1231 ptr_out = astGetPoints( result );
1232
1233 /* Determine whether to apply the forward or inverse mapping, according to the
1234 direction specified and whether the mapping has been inverted. */
1235 if ( astGetInvert( map ) ) forward = !forward;
1236
1237 /* Perform coordinate arithmetic. */
1238 /* ------------------------------ */
1239 if ( astOK ) {
1240
1241 /* First deal with forward transformations. */
1242 if( forward ){
1243
1244 /* Store pointers to the next value on each axis. */
1245 xx = ptr_in[ 0 ];
1246 yy = ptr_in[ 1 ];
1247 aa = ptr_out[ 0 ];
1248 bb = ptr_out[ 1 ];
1249
1250 /* Loop to apply the plate fit to all the points, checking for (and
1251 propagating) bad values in the process. */
1252 for ( point = 0; point < npoint; point++ ) {
1253 if( *xx != AST__BAD && *yy != AST__BAD ){
1254
1255 /* If the pixel position is transformed succesfully, convert the returned
1256 RA/DEC from degrees to radians. Otherwise, store bad values. NB,
1257 platepos returns zero for success. */
1258 if( !platepos( *xx, *yy, (struct WorldCoor *) map->wcs,
1259 aa, bb ) ){
1260 (*aa) *= AST__DD2R;
1261 (*bb) *= AST__DD2R;
1262
1263 } else {
1264 *aa = AST__BAD;
1265 *bb = AST__BAD;
1266 }
1267
1268 } else {
1269 *aa = AST__BAD;
1270 *bb = AST__BAD;
1271 }
1272
1273 /* Move on to the next point. */
1274 xx++;
1275 yy++;
1276 aa++;
1277 bb++;
1278 }
1279
1280 /* Now deal with inverse transformations in the same way. */
1281 } else {
1282 aa = ptr_in[ 0 ];
1283 bb = ptr_in[ 1 ];
1284 xx = ptr_out[ 0 ];
1285 yy = ptr_out[ 1 ];
1286
1287 for ( point = 0; point < npoint; point++ ) {
1288 if( *aa != AST__BAD && *bb != AST__BAD ){
1289
1290 if( platepix( AST__DR2D*(*aa), AST__DR2D*(*bb),
1291 (struct WorldCoor *) map->wcs, xx, yy ) ){
1292 *xx = AST__BAD;
1293 *yy = AST__BAD;
1294 }
1295
1296 } else {
1297 *xx = AST__BAD;
1298 *yy = AST__BAD;
1299 }
1300
1301 xx++;
1302 yy++;
1303 aa++;
1304 bb++;
1305 }
1306 }
1307 }
1308
1309 /* Return a pointer to the output PointSet. */
1310 return result;
1311 }
1312
1313 /* Functions which access class attributes. */
1314 /* ---------------------------------------- */
1315 /* Implement member functions to access the attributes associated with
1316 this class using the macros defined for this purpose in the
1317 "object.h" file. For a description of each attribute, see the class
1318 interface (in the associated .h file). */
1319
1320 /* Copy constructor. */
1321 /* ----------------- */
Copy(const AstObject * objin,AstObject * objout,int * status)1322 static void Copy( const AstObject *objin, AstObject *objout, int *status ) {
1323 /*
1324 * Name:
1325 * Copy
1326
1327 * Purpose:
1328 * Copy constructor for DssMap objects.
1329
1330 * Type:
1331 * Private function.
1332
1333 * Synopsis:
1334 * void Copy( const AstObject *objin, AstObject *objout, int *status )
1335
1336 * Description:
1337 * This function implements the copy constructor for DssMap objects.
1338
1339 * Parameters:
1340 * objin
1341 * Pointer to the object to be copied.
1342 * objout
1343 * Pointer to the object being constructed.
1344 * status
1345 * Pointer to the inherited status variable.
1346
1347 * Returned Value:
1348 * void
1349
1350 * Notes:
1351 * - This constructor makes a deep copy.
1352 */
1353
1354
1355 /* Local Variables: */
1356 AstDssMap *in; /* Pointer to input DssMap */
1357 AstDssMap *out; /* Pointer to output DssMap */
1358
1359 /* Check the global error status. */
1360 if ( !astOK ) return;
1361
1362 /* Obtain pointers to the input and output DssMaps. */
1363 in = (AstDssMap *) objin;
1364 out = (AstDssMap *) objout;
1365
1366 /* Store a copy of the input SAOIMAGE WorldCoor structure in the output. */
1367 out->wcs = astStore( NULL, in->wcs, sizeof( struct WorldCoor ) );
1368
1369 return;
1370
1371 }
1372
1373 /* Destructor. */
1374 /* ----------- */
Delete(AstObject * obj,int * status)1375 static void Delete( AstObject *obj, int *status ) {
1376 /*
1377 * Name:
1378 * Delete
1379
1380 * Purpose:
1381 * Destructor for DssMap objects.
1382
1383 * Type:
1384 * Private function.
1385
1386 * Synopsis:
1387 * void Delete( AstObject *obj, int *status )
1388
1389 * Description:
1390 * This function implements the destructor for DssMap objects.
1391
1392 * Parameters:
1393 * obj
1394 * Pointer to the object to be deleted.
1395 * status
1396 * Pointer to the inherited status variable.
1397
1398 * Returned Value:
1399 * void
1400
1401 * Notes:
1402 * This function attempts to execute even if the global error status is
1403 * set.
1404 */
1405
1406 /* Local Variables: */
1407 AstDssMap *this; /* Pointer to DssMap */
1408
1409 /* Obtain a pointer to the DssMap structure. */
1410 this = (AstDssMap *) obj;
1411
1412 /* Free the SAOIMAGE WorldCoor structure. */
1413 this->wcs = astFree( this->wcs );
1414
1415 }
1416
1417 /* Dump function. */
1418 /* -------------- */
Dump(AstObject * this_object,AstChannel * channel,int * status)1419 static void Dump( AstObject *this_object, AstChannel *channel, int *status ) {
1420 /*
1421 * Name:
1422 * Dump
1423
1424 * Purpose:
1425 * Dump function for DssMap objects.
1426
1427 * Type:
1428 * Private function.
1429
1430 * Synopsis:
1431 * void Dump( AstObject *this, AstChannel *channel, int *status )
1432
1433 * Description:
1434 * This function implements the Dump function which writes out data
1435 * for the DssMap class to an output Channel.
1436
1437 * Parameters:
1438 * this
1439 * Pointer to the DssMap whose data are being written.
1440 * channel
1441 * Pointer to the Channel to which the data are being written.
1442 * status
1443 * Pointer to the inherited status variable.
1444 */
1445
1446 AstDssMap *this; /* Pointer to the DssMap structure */
1447 struct WorldCoor *wcs; /* Pointer to SAOimage wcs structure */
1448 char name_buff[ 11 ]; /* Buffer for keyword string */
1449 int i; /* Coefficient index */
1450
1451 /* Check the global error status. */
1452 if ( !astOK ) return;
1453
1454 /* Obtain a pointer to the DssMap structure. */
1455 this = (AstDssMap *) this_object;
1456
1457 /* Get a pointer to the WorldCoor structure holding the description of the
1458 DssMap. */
1459 wcs = (struct WorldCoor *) ( this->wcs );
1460
1461 /* Write out values representing the contents of the WorldCoor structure.
1462 Only the components which are required to re-create the DssMap are
1463 written out. */
1464 astWriteDouble( channel, "PltRA", 1, 1, wcs->plate_ra, "Plate centre RA (radians)" );
1465 astWriteDouble( channel, "PltDec", 1, 1, wcs->plate_dec, "Plate centre Dec (radians)" );
1466 astWriteDouble( channel, "PltScl", 1, 1, wcs->plate_scale, "Plate scale (arcsec/mm)" );
1467 astWriteDouble( channel, "CNPix1", 1, 1, wcs->x_pixel_offset, "X Pixel offset (pixels)" );
1468 astWriteDouble( channel, "CNPix2", 1, 1, wcs->y_pixel_offset, "Y Pixel offset (pixels)" );
1469 astWriteDouble( channel, "XPixSz", 1, 1, wcs->x_pixel_size, "X Pixel size (microns)" );
1470 astWriteDouble( channel, "YPixSz", 1, 1, wcs->y_pixel_size, "Y Pixel size (microns)" );
1471
1472 for( i = 0; i < 6; i++ ) {
1473 sprintf( name_buff, "PPO%d", i + 1 );
1474 astWriteDouble( channel, name_buff, 1, 1, wcs->ppo_coeff[i],
1475 "Orientation coefficients" );
1476 }
1477
1478 for( i = 0; i < 19; i++ ) {
1479 sprintf( name_buff, "AMDX%d", i + 1 );
1480 astWriteDouble( channel, name_buff, 1, 1, wcs->amd_x_coeff[i],
1481 "Plate solution X coefficients" );
1482 }
1483
1484 for( i = 0; i < 19; i++ ) {
1485 sprintf( name_buff, "AMDY%d", i + 1 );
1486 astWriteDouble( channel, name_buff, 1, 1, wcs->amd_y_coeff[i],
1487 "Plate solution Y coefficients" );
1488 }
1489
1490 }
1491
1492 /* Standard class functions. */
1493 /* ========================= */
1494 /* Implement the astIsADssMap and astCheckDssMap functions using the macros
1495 defined for this purpose in the "object.h" header file. */
astMAKE_ISA(DssMap,Mapping)1496 astMAKE_ISA(DssMap,Mapping)
1497 astMAKE_CHECK(DssMap)
1498
1499 AstDssMap *astDssMap_( void *fits_void, const char *options, int *status, ...) {
1500 /*
1501 *+
1502 * Name:
1503 * astDssMap
1504
1505 * Purpose:
1506 * Create a DssMap.
1507
1508 * Type:
1509 * Protected function.
1510
1511 * Synopsis:
1512 * #include "dssmap.h"
1513 * AstDssMap *astDssMap( AstFitsChan *fits, const char *options, int *status, ... )
1514
1515 * Class Membership:
1516 * DssMap constructor.
1517
1518 * Description:
1519 * This function creates a new DssMap and optionally initialises its
1520 * attributes.
1521 *
1522 * A DssMap is a Mapping which uses a Digitised Sky Survey plate fit to
1523 * transform a set of points from pixel coordinates to equatorial
1524 * coordinates (i.e. Right Ascension and Declination).
1525
1526 * Parameters:
1527 * fits
1528 * A pointer to a FitsChan holding a set of FITS header cards
1529 * describing the plate fit to be used. The FitsChan may contain
1530 * other header cards which will be ignored, and it is unchanged on
1531 * exit. The required information is copied from the FitsChan, and
1532 * so the supplied FitsChan may subsequently be changed or deleted
1533 * without changing the DssMap.
1534 * options
1535 * Pointer to a null-terminated string containing an optional
1536 * comma-separated list of attribute assignments to be used for
1537 * initialising the new DssMap. The syntax used is identical to
1538 * that for the astSet function and may include "printf" format
1539 * specifiers identified by "%" symbols in the normal way.
1540 * status
1541 * Pointer to the inherited status variable.
1542 * ...
1543 * If the "options" string contains "%" format specifiers, then
1544 * an optional list of additional arguments may follow it in
1545 * order to supply values to be substituted for these
1546 * specifiers. The rules for supplying these are identical to
1547 * those for the astSet function (and for the C "printf"
1548 * function).
1549
1550 * Returned Value:
1551 * astDssMap()
1552 * A pointer to the new DssMap.
1553
1554 * Attributes:
1555 * The DssMap class has no additional attributes over and above those
1556 * common to all Mappings.
1557
1558 * Notes:
1559 * - The supplied FitsChan must contain values for the following FITS
1560 * keywords: CNPIX1, CNPIX2, PPO3, PPO6, XPIXELSZ, YPIXELSZ, PLTRAH,
1561 * PLTRAM, PLTRAS, PLTDECD, PLTDECM, PLTDECS, PLTDECSN, PLTSCALE,
1562 * AMDX1, AMDX2, ..., AMDX13, AMDY1, AMDY2, ..., AMDY13.
1563 * - A null Object pointer (AST__NULL) will be returned if this
1564 * function is invoked with the AST error status set, or if it
1565 * should fail for any reason.
1566 *-
1567 */
1568
1569 /* Local Variables: */
1570 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
1571 AstFitsChan *fits; /* Pointer to supplied FitsChan */
1572 AstDssMap *new; /* Pointer to new DssMap */
1573 va_list args; /* Variable argument list */
1574
1575 /* Get a pointer to the thread specific global data structure. */
1576 astGET_GLOBALS(NULL);
1577
1578 /* Check the global status. */
1579 new = NULL;
1580 if ( !astOK ) return new;
1581
1582 /* Obtain and validate a pointer to the FitsChan structure provided. */
1583 fits = astCheckFitsChan( fits_void );
1584 if ( astOK ) {
1585
1586 /* Initialise the DssMap, allocating memory and initialising the
1587 virtual function table as well if necessary. */
1588 new = astInitDssMap( NULL, sizeof( AstDssMap ), !class_init, &class_vtab,
1589 "DssMap", fits );
1590
1591 /* If successful, note that the virtual function table has been
1592 initialised. */
1593 if ( astOK ) {
1594 class_init = 1;
1595
1596 /* Obtain the variable argument list and pass it along with the options string
1597 to the astVSet method to initialise the new DssMap's attributes. */
1598 va_start( args, status );
1599 astVSet( new, options, NULL, args );
1600 va_end( args );
1601
1602 /* If an error occurred, clean up by deleting the new object. */
1603 if ( !astOK ) new = astDelete( new );
1604 }
1605 }
1606
1607 /* Return a pointer to the new DssMap. */
1608 return new;
1609 }
1610
astInitDssMap_(void * mem,size_t size,int init,AstDssMapVtab * vtab,const char * name,AstFitsChan * fits,int * status)1611 AstDssMap *astInitDssMap_( void *mem, size_t size, int init,
1612 AstDssMapVtab *vtab, const char *name,
1613 AstFitsChan *fits, int *status ) {
1614 /*
1615 *+
1616 * Name:
1617 * astInitDssMap
1618
1619 * Purpose:
1620 * Initialise a DssMap.
1621
1622 * Type:
1623 * Protected function.
1624
1625 * Synopsis:
1626 * #include "dssmap.h"
1627 * AstDssMap *astInitDssMap( void *mem, size_t size, int init,
1628 * AstDssMapVtab *vtab, const char *name,
1629 * AstFitsChan *fits )
1630
1631 * Class Membership:
1632 * DssMap initialiser.
1633
1634 * Description:
1635 * This function is provided for use by class implementations to initialise
1636 * a new DssMap object. It allocates memory (if necessary) to accommodate
1637 * the DssMap plus any additional data associated with the derived class.
1638 * It then initialises a DssMap structure at the start of this memory. If
1639 * the "init" flag is set, it also initialises the contents of a virtual
1640 * function table for a DssMap at the start of the memory passed via the
1641 * "vtab" parameter.
1642
1643 * Parameters:
1644 * mem
1645 * A pointer to the memory in which the DssMap is to be initialised.
1646 * This must be of sufficient size to accommodate the DssMap data
1647 * (sizeof(DssMap)) plus any data used by the derived class. If a value
1648 * of NULL is given, this function will allocate the memory itself using
1649 * the "size" parameter to determine its size.
1650 * size
1651 * The amount of memory used by the DssMap (plus derived class data).
1652 * This will be used to allocate memory if a value of NULL is given for
1653 * the "mem" parameter. This value is also stored in the DssMap
1654 * structure, so a valid value must be supplied even if not required for
1655 * allocating memory.
1656 * init
1657 * A logical flag indicating if the DssMap's virtual function table is
1658 * to be initialised. If this value is non-zero, the virtual function
1659 * table will be initialised by this function.
1660 * vtab
1661 * Pointer to the start of the virtual function table to be associated
1662 * with the new DssMap.
1663 * name
1664 * Pointer to a constant null-terminated character string which contains
1665 * the name of the class to which the new object belongs (it is this
1666 * pointer value that will subsequently be returned by the astGetClass
1667 * method).
1668 * fits
1669 * Pointer to a FitsChan containing the DSS FITS Header.
1670
1671 * Returned Value:
1672 * A pointer to the new DssMap.
1673
1674 * Notes:
1675 * - A null pointer will be returned if this function is invoked with the
1676 * global error status set, or if it should fail for any reason.
1677 *-
1678 */
1679
1680 /* Local Variables: */
1681 AstDssMap *new; /* Pointer to new DssMap */
1682 struct WorldCoor *wcs; /* Pointer to SAOIMAGE wcs structure */
1683
1684 /* Check the global status. */
1685 if ( !astOK ) return NULL;
1686
1687 /* If necessary, initialise the virtual function table. */
1688 if ( init ) astInitDssMapVtab( vtab, name );
1689
1690 /* Initialise. */
1691 new = NULL;
1692
1693 /* Create a structure holding the information required by the SAOIMAGE
1694 "platepos" function. The required values are extracted from the
1695 supplied FitsChan. An error is reported and NULL returned if any required
1696 keywords are missing or unusable. */
1697 if ( ( wcs = BuildWcs( fits, "astInitDssMap", name, status ) ) ) {
1698
1699 /* Initialise a 2-D Mapping structure (the parent class) as the first component
1700 within the DssMap structure, allocating memory if necessary. Specify that
1701 the Mapping should be defined in both the forward and inverse directions. */
1702 new = (AstDssMap *) astInitMapping( mem, size, 0,
1703 (AstMappingVtab *) vtab, name,
1704 2, 2, 1, 1 );
1705
1706 if ( astOK ) {
1707
1708 /* Initialise the DssMap data. */
1709 /* --------------------------- */
1710 /* Store a copy of the SAOIMAGE wcs structure. */
1711 new->wcs = astStore( NULL, (void *) wcs, sizeof( struct WorldCoor ) );
1712
1713 /* If an error occurred, clean up by deleting the new DssMap. */
1714 if ( !astOK ) new = astDelete( new );
1715 }
1716
1717 /* Free the SAOIMAGE wcs structure. */
1718 wcs = (struct WorldCoor *) astFree( (void *) wcs );
1719
1720 }
1721
1722 /* Return a pointer to the new DssMap. */
1723 return new;
1724 }
1725
astLoadDssMap_(void * mem,size_t size,AstDssMapVtab * vtab,const char * name,AstChannel * channel,int * status)1726 AstDssMap *astLoadDssMap_( void *mem, size_t size,
1727 AstDssMapVtab *vtab, const char *name,
1728 AstChannel *channel, int *status ) {
1729 /*
1730 *+
1731 * Name:
1732 * astLoadDssMap
1733
1734 * Purpose:
1735 * Load a DssMap.
1736
1737 * Type:
1738 * Protected function.
1739
1740 * Synopsis:
1741 * #include "dssmap.h"
1742 * AstDssMap *astLoadDssMap( void *mem, size_t size,
1743 * AstDssMapVtab *vtab, const char *name,
1744 * AstChannel *channel )
1745
1746 * Class Membership:
1747 * DssMap loader.
1748
1749 * Description:
1750 * This function is provided to load a new DssMap using data read
1751 * from a Channel. It first loads the data used by the parent class
1752 * (which allocates memory if necessary) and then initialises a
1753 * DssMap structure in this memory, using data read from the input
1754 * Channel.
1755 *
1756 * If the "init" flag is set, it also initialises the contents of a
1757 * virtual function table for a DssMap at the start of the memory
1758 * passed via the "vtab" parameter.
1759
1760
1761 * Parameters:
1762 * mem
1763 * A pointer to the memory into which the DssMap is to be
1764 * loaded. This must be of sufficient size to accommodate the
1765 * DssMap data (sizeof(DssMap)) plus any data used by derived
1766 * classes. If a value of NULL is given, this function will
1767 * allocate the memory itself using the "size" parameter to
1768 * determine its size.
1769 * size
1770 * The amount of memory used by the DssMap (plus derived class
1771 * data). This will be used to allocate memory if a value of
1772 * NULL is given for the "mem" parameter. This value is also
1773 * stored in the DssMap structure, so a valid value must be
1774 * supplied even if not required for allocating memory.
1775 *
1776 * If the "vtab" parameter is NULL, the "size" value is ignored
1777 * and sizeof(AstDssMap) is used instead.
1778 * vtab
1779 * Pointer to the start of the virtual function table to be
1780 * associated with the new DssMap. If this is NULL, a pointer
1781 * to the (static) virtual function table for the DssMap class
1782 * is used instead.
1783 * name
1784 * Pointer to a constant null-terminated character string which
1785 * contains the name of the class to which the new object
1786 * belongs (it is this pointer value that will subsequently be
1787 * returned by the astGetClass method).
1788 *
1789 * If the "vtab" parameter is NULL, the "name" value is ignored
1790 * and a pointer to the string "DssMap" is used instead.
1791
1792 * Returned Value:
1793 * A pointer to the new DssMap.
1794
1795 * Notes:
1796 * - A null pointer will be returned if this function is invoked
1797 * with the global error status set, or if it should fail for any
1798 * reason.
1799 *-
1800 */
1801
1802 /* Local Variables: */
1803 astDECLARE_GLOBALS /* Pointer to thread-specific global data */
1804 AstDssMap *new; /* Pointer to the new DssMap */
1805 char name_buff[ 11 ]; /* Buffer for item name */
1806 int i; /* Coefficient index */
1807 struct WorldCoor *wcs; /* Pointer to Wcs information */
1808
1809 /* Get a pointer to the thread specific global data structure. */
1810 astGET_GLOBALS(channel);
1811
1812 /* Initialise. */
1813 new = NULL;
1814
1815 /* Check the global error status. */
1816 if ( !astOK ) return new;
1817
1818 /* If a NULL virtual function table has been supplied, then this is
1819 the first loader to be invoked for this DssMap. In this case the
1820 DssMap belongs to this class, so supply appropriate values to be
1821 passed to the parent class loader (and its parent, etc.). */
1822 if ( !vtab ) {
1823 size = sizeof( AstDssMap );
1824 vtab = &class_vtab;
1825 name = "DssMap";
1826
1827 /* If required, initialise the virtual function table for this class. */
1828 if ( !class_init ) {
1829 astInitDssMapVtab( vtab, name );
1830 class_init = 1;
1831 }
1832 }
1833
1834 /* Invoke the parent class loader to load data for all the ancestral
1835 classes of the current one, returning a pointer to the resulting
1836 partly-built DssMap. */
1837 new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name,
1838 channel );
1839
1840 if ( astOK ) {
1841
1842
1843 /* Read input data. */
1844 /* ================ */
1845 /* Request the input Channel to read all the input data appropriate to
1846 this class into the internal "values list". */
1847 astReadClassData( channel, "DssMap" );
1848
1849 /* Allocate memory to hold the WorldCoor structure which holds the wcs
1850 information in a form usable by the SAOimage projection functions. */
1851 new->wcs = astMalloc( sizeof(struct WorldCoor) );
1852 if( astOK ) {
1853
1854 /* Get a pointer to the WorldCoor structure holding the description of the
1855 DssMap. */
1856 wcs = (struct WorldCoor *) ( new->wcs );
1857
1858 /* Read the values representing the contents of the WorldCoor structure.
1859 Only the components which are required to re-create the DssMap are
1860 read. */
1861 wcs->plate_ra = astReadDouble( channel, "pltra", AST__BAD );
1862 if( wcs->plate_ra == AST__BAD && astOK ){
1863 astError( AST__RDERR, "astRead(DssMap): 'PltRA' object (Plate "
1864 "centre RA) missing from input." , status);
1865 }
1866
1867 wcs->plate_dec = astReadDouble( channel, "pltdec", AST__BAD );
1868 if( wcs->plate_dec == AST__BAD && astOK ){
1869 astError( AST__RDERR, "astRead(DssMap): 'PltDec' object (Plate "
1870 "centre Dec) missing from input." , status);
1871 }
1872
1873 wcs->plate_scale = astReadDouble( channel, "pltscl", AST__BAD );
1874 if( wcs->plate_scale == AST__BAD && astOK ){
1875 astError( AST__RDERR, "astRead(DssMap): 'PltScl' object (Plate "
1876 "scale) missing from input." , status);
1877 }
1878
1879 wcs->x_pixel_offset = astReadDouble( channel, "cnpix1", AST__BAD );
1880 if( wcs->x_pixel_offset == AST__BAD && astOK ){
1881 astError( AST__RDERR, "astRead(DssMap): 'CNPix1' object (X pixel "
1882 "offset) missing from input." , status);
1883 }
1884
1885 wcs->y_pixel_offset = astReadDouble( channel, "cnpix2", AST__BAD );
1886 if( wcs->y_pixel_offset == AST__BAD && astOK ){
1887 astError( AST__RDERR, "astRead(DssMap): 'CNPix2' object (Y pixel "
1888 "offset) missing from input." , status);
1889 }
1890
1891 wcs->x_pixel_size = astReadDouble( channel, "xpixsz", AST__BAD );
1892 if( wcs->x_pixel_size == AST__BAD && astOK ){
1893 astError( AST__RDERR, "astRead(DssMap): 'XPixSz' object (X pixel "
1894 "size) missing from input." , status);
1895 }
1896
1897 wcs->y_pixel_size = astReadDouble( channel, "ypixsz", AST__BAD );
1898 if( wcs->y_pixel_size == AST__BAD && astOK ){
1899 astError( AST__RDERR, "astRead(DssMap): 'YPixSz' object (Y pixel "
1900 "size) missing from input." , status);
1901 }
1902
1903 for( i = 0; i < 6 && astOK; i++ ) {
1904 sprintf( name_buff, "ppo%d", i + 1 );
1905 wcs->ppo_coeff[i] = astReadDouble( channel, name_buff, AST__BAD );
1906 if( wcs->ppo_coeff[i] == AST__BAD ){
1907 if( i == 2 || i == 5 ) {
1908 if( astOK ) astError( AST__RDERR, "astRead(DssMap): 'PPO%d' "
1909 "object (orientation coefficient %d) "
1910 "missing from input.", status, i + 1, i + 1 );
1911 } else {
1912 wcs->ppo_coeff[i] = 0.0;
1913 }
1914 }
1915 }
1916
1917 for( i = 0; i < 19 && astOK; i++ ) {
1918 sprintf( name_buff, "amdx%d", i + 1 );
1919 wcs->amd_x_coeff[i] = astReadDouble( channel, name_buff, AST__BAD );
1920 if( wcs->amd_x_coeff[i] == AST__BAD ){
1921 if( i < 13 ){
1922 if( astOK ) astError( AST__RDERR, "astRead(DssMap): 'AMDX%d' "
1923 "object (plate solution X coefficient "
1924 "%d) missing from input.", status, i + 1, i + 1 );
1925 } else {
1926 wcs->amd_x_coeff[i] = 0.0;
1927 }
1928 }
1929 }
1930
1931 for( i = 0; i < 19 && astOK; i++ ) {
1932 sprintf( name_buff, "amdy%d", i + 1 );
1933 wcs->amd_y_coeff[i] = astReadDouble( channel, name_buff, AST__BAD );
1934 if( wcs->amd_y_coeff[i] == AST__BAD ){
1935 if( i < 13 ){
1936 if( astOK ) astError( AST__RDERR, "astRead(DssMap): 'AMDY%d' "
1937 "object (plate solution Y coefficient "
1938 "%d) missing from input.", status, i + 1, i + 1 );
1939 } else {
1940 wcs->amd_y_coeff[i] = 0.0;
1941 }
1942 }
1943 }
1944 }
1945
1946 /* If an error occurred, clean up by deleting the new DssMap. */
1947 if ( !astOK ) new = astDelete( new );
1948 }
1949
1950 /* Return the new DssMap pointer. */
1951 return new;
1952 }
1953
1954 /* Virtual function interfaces. */
1955 /* ============================ */
1956 /* These provide the external interface to the virtual functions defined by
1957 this class. Each simply checks the global error status and then locates and
1958 executes the appropriate member function, using the function pointer stored
1959 in the object's virtual function table (this pointer is located using the
1960 astMEMBER macro defined in "object.h").
1961
1962 Note that the member function may not be the one defined here, as it may
1963 have been over-ridden by a derived class. However, it should still have the
1964 same interface. */
astDssFits_(AstDssMap * this,int * status)1965 AstFitsChan *astDssFits_( AstDssMap *this, int *status ){
1966 if( !astOK ) return NULL;
1967 return (**astMEMBER(this,DssMap,DssFits))( this, status );
1968 }
1969
1970 /* The code which follows in this file is covered by the following
1971 statement of terms and conditions, which differ from the terms and
1972 conditions which apply above.
1973
1974 ***************************************************************************
1975 *
1976 * Copyright: 1988 Smithsonian Astrophysical Observatory
1977 * You may do anything you like with these files except remove
1978 * this copyright. The Smithsonian Astrophysical Observatory
1979 * makes no representations about the suitability of this
1980 * software for any purpose. It is provided "as is" without
1981 * express or implied warranty.
1982 *
1983 *****************************************************************************
1984 */
1985
1986 /* >>>>>>>>>>>>>>>>>>>>>> platepos.c <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
1987
1988 /* File saoimage/wcslib/platepos.c
1989 * February 25, 1996
1990 * By Doug Mink, Harvard-Smithsonian Center for Astrophysics
1991
1992 * Module: platepos.c (Plate solution WCS conversion
1993 * Purpose: Compute WCS from Digital Sky Survey plate fit
1994 * Subroutine: platepos() converts from pixel location to RA,Dec
1995 * Subroutine: platepix() converts from RA,Dec to pixel location
1996
1997 These functions are based on the astrmcal.c portion of GETIMAGE by
1998 J. Doggett and the documentation distributed with the Digital Sky Survey.
1999
2000 >>>>>>> STARLINK VERSION <<<<<<<<
2001
2002 */
2003
2004 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
2005
2006 Changed by R.F. Warren-Smith (Starlink) to make the function static. */
2007
2008 static int
platepos(xpix,ypix,wcs,xpos,ypos)2009 platepos (xpix, ypix, wcs, xpos, ypos)
2010
2011 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> */
2012
2013 /* Routine to determine accurate position for pixel coordinates */
2014 /* returns 0 if successful otherwise 1 = angle too large for projection; */
2015 /* based on amdpos() from getimage */
2016
2017 /* Input: */
2018 double xpix; /* x pixel number (RA or long without rotation) */
2019 double ypix; /* y pixel number (dec or lat without rotation) */
2020 struct WorldCoor *wcs; /* WCS parameter structure */
2021
2022 /* Output: */
2023 double *xpos; /* Right ascension or longitude in degrees */
2024 double *ypos; /* Declination or latitude in degrees */
2025
2026 {
2027 double x, y, xmm, ymm, xmm2, ymm2, xmm3, ymm3, x2y2;
2028 double xi, xir, eta, etar, raoff, ra, dec;
2029 double cond2r = 1.745329252e-2;
2030 double cons2r = 206264.8062470964;
2031 double twopi = 6.28318530717959;
2032 double ctan, ccos;
2033
2034 /* Ignore magnitude and color terms
2035 double mag = 0.0;
2036 double color = 0.0; */
2037
2038 /* Convert from image pixels to plate pixels */
2039 x = xpix + wcs->x_pixel_offset - 1.0 + 0.5;
2040 y = ypix + wcs->y_pixel_offset - 1.0 + 0.5;
2041
2042 /* Convert from pixels to millimeters */
2043 xmm = (wcs->ppo_coeff[2] - x * wcs->x_pixel_size) / 1000.0;
2044 ymm = (y * wcs->y_pixel_size - wcs->ppo_coeff[5]) / 1000.0;
2045 xmm2 = xmm * xmm;
2046 ymm2 = ymm * ymm;
2047 xmm3 = xmm * xmm2;
2048 ymm3 = ymm * ymm2;
2049 x2y2 = xmm2 + ymm2;
2050
2051 /* Compute coordinates from x,y and plate model */
2052
2053 xi = wcs->amd_x_coeff[ 0]*xmm + wcs->amd_x_coeff[ 1]*ymm +
2054 wcs->amd_x_coeff[ 2] + wcs->amd_x_coeff[ 3]*xmm2 +
2055 wcs->amd_x_coeff[ 4]*xmm*ymm + wcs->amd_x_coeff[ 5]*ymm2 +
2056 wcs->amd_x_coeff[ 6]*(x2y2) + wcs->amd_x_coeff[ 7]*xmm3 +
2057 wcs->amd_x_coeff[ 8]*xmm2*ymm + wcs->amd_x_coeff[ 9]*xmm*ymm2 +
2058 wcs->amd_x_coeff[10]*ymm3 + wcs->amd_x_coeff[11]*xmm*(x2y2) +
2059 wcs->amd_x_coeff[12]*xmm*x2y2*x2y2;
2060
2061 /* Ignore magnitude and color terms
2062 + wcs->amd_x_coeff[13]*mag + wcs->amd_x_coeff[14]*mag*mag +
2063 wcs->amd_x_coeff[15]*mag*mag*mag + wcs->amd_x_coeff[16]*mag*xmm +
2064 wcs->amd_x_coeff[17]*mag*x2y2 + wcs->amd_x_coeff[18]*mag*xmm*x2y2 +
2065 wcs->amd_x_coeff[19]*color; */
2066
2067 eta = wcs->amd_y_coeff[ 0]*ymm + wcs->amd_y_coeff[ 1]*xmm +
2068 wcs->amd_y_coeff[ 2] + wcs->amd_y_coeff[ 3]*ymm2 +
2069 wcs->amd_y_coeff[ 4]*xmm*ymm + wcs->amd_y_coeff[ 5]*xmm2 +
2070 wcs->amd_y_coeff[ 6]*(x2y2) + wcs->amd_y_coeff[ 7]*ymm3 +
2071 wcs->amd_y_coeff[ 8]*ymm2*xmm + wcs->amd_y_coeff[ 9]*ymm*xmm2 +
2072 wcs->amd_y_coeff[10]*xmm3 + wcs->amd_y_coeff[11]*ymm*(x2y2) +
2073 wcs->amd_y_coeff[12]*ymm*x2y2*x2y2;
2074
2075 /* Ignore magnitude and color terms
2076 + wcs->amd_y_coeff[13]*mag + wcs->amd_y_coeff[14]*mag*mag +
2077 wcs->amd_y_coeff[15]*mag*mag*mag + wcs->amd_y_coeff[16]*mag*ymm +
2078 wcs->amd_y_coeff[17]*mag*x2y2) + wcs->amd_y_coeff[18]*mag*ymm*x2y2 +
2079 wcs->amd_y_coeff[19]*color; */
2080
2081 /* Convert to radians */
2082
2083 xir = xi / cons2r;
2084 etar = eta / cons2r;
2085
2086 /* Convert to RA and Dec */
2087
2088 ctan = tan (wcs->plate_dec);
2089 ccos = cos (wcs->plate_dec);
2090 raoff = atan2 (xir / ccos, 1.0 - etar * ctan);
2091 ra = raoff + wcs->plate_ra;
2092 if (ra < 0.0) ra = ra + twopi;
2093 *xpos = ra / cond2r;
2094
2095 dec = atan (cos (raoff) / ((1.0 - (etar * ctan)) / (etar + ctan)));
2096 *ypos = dec / cond2r;
2097 return 0;
2098 }
2099
2100
2101 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
2102
2103 Changed by R.F. Warren-Smith (Starlink) to make the function static. */
2104
2105 static int
platepix(xpos,ypos,wcs,xpix,ypix)2106 platepix (xpos, ypos, wcs, xpix, ypix)
2107
2108 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> */
2109
2110 /* Routine to determine pixel coordinates for sky position */
2111 /* returns 0 if successful otherwise 1 = angle too large for projection; */
2112 /* based on amdinv() from getimage */
2113
2114 /* Input: */
2115 double xpos; /* Right ascension or longitude in degrees */
2116 double ypos; /* Declination or latitude in degrees */
2117 struct WorldCoor *wcs; /* WCS parameter structure */
2118
2119 /* Output: */
2120 double *xpix; /* x pixel number (RA or long without rotation) */
2121 double *ypix; /* y pixel number (dec or lat without rotation) */
2122
2123 {
2124 double div,xi,eta,x,y,xy,x2,y2,x2y,y2x,x3,y3,x4,y4,x2y2,cjunk,dx,dy;
2125 double sypos,cypos,syplate,cyplate,sxdiff,cxdiff;
2126 double f,fx,fy,g,gx,gy, xmm, ymm;
2127 double conr2s = 206264.8062470964;
2128 double tolerance = 0.0000005;
2129 int max_iterations = 50;
2130 int i;
2131 double xr, yr; /* position in radians */
2132 double cond2r = 1.745329252e-2;
2133
2134 /* Convert RA and Dec in radians to standard coordinates on a plate */
2135 xr = xpos * cond2r;
2136 yr = ypos * cond2r;
2137 sypos = sin (yr);
2138 cypos = cos (yr);
2139 syplate = sin (wcs->plate_dec);
2140 cyplate = cos (wcs->plate_dec);
2141 sxdiff = sin (xr - wcs->plate_ra);
2142 cxdiff = cos (xr - wcs->plate_ra);
2143 div = (sypos * syplate) + (cypos * cyplate * cxdiff);
2144 xi = cypos * sxdiff * conr2s / div;
2145 eta = ((sypos * cyplate) - (cypos * syplate * cxdiff)) * conr2s / div;
2146
2147 /* Set initial value for x,y */
2148 xmm = xi / wcs->plate_scale;
2149 ymm = eta / wcs->plate_scale;
2150
2151 /* Iterate by Newton's method */
2152 for (i = 0; i < max_iterations; i++) {
2153
2154 /* X plate model */
2155 xy = xmm * ymm;
2156 x2 = xmm * xmm;
2157 y2 = ymm * ymm;
2158 x2y = x2 * ymm;
2159 y2x = y2 * xmm;
2160 x2y2 = x2 + y2;
2161 cjunk = x2y2 * x2y2;
2162 x3 = x2 * xmm;
2163 y3 = y2 * ymm;
2164 x4 = x2 * x2;
2165 y4 = y2 * y2;
2166 f = wcs->amd_x_coeff[0]*xmm + wcs->amd_x_coeff[1]*ymm +
2167 wcs->amd_x_coeff[2] + wcs->amd_x_coeff[3]*x2 +
2168 wcs->amd_x_coeff[4]*xy + wcs->amd_x_coeff[5]*y2 +
2169 wcs->amd_x_coeff[6]*x2y2 + wcs->amd_x_coeff[7]*x3 +
2170 wcs->amd_x_coeff[8]*x2y + wcs->amd_x_coeff[9]*y2x +
2171 wcs->amd_x_coeff[10]*y3 + wcs->amd_x_coeff[11]*xmm*x2y2 +
2172 wcs->amd_x_coeff[12]*xmm*cjunk;
2173 /* magnitude and color terms ignored
2174 + wcs->amd_x_coeff[13]*mag +
2175 wcs->amd_x_coeff[14]*mag*mag + wcs->amd_x_coeff[15]*mag*mag*mag +
2176 wcs->amd_x_coeff[16]*mag*xmm + wcs->amd_x_coeff[17]*mag*(x2+y2) +
2177 wcs->amd_x_coeff[18]*mag*xmm*(x2+y2) + wcs->amd_x_coeff[19]*color;
2178 */
2179
2180 /* Derivative of X model wrt x */
2181 fx = wcs->amd_x_coeff[0] + wcs->amd_x_coeff[3]*2.0*xmm +
2182 wcs->amd_x_coeff[4]*ymm + wcs->amd_x_coeff[6]*2.0*xmm +
2183 wcs->amd_x_coeff[7]*3.0*x2 + wcs->amd_x_coeff[8]*2.0*xy +
2184 wcs->amd_x_coeff[9]*y2 + wcs->amd_x_coeff[11]*(3.0*x2+y2) +
2185 wcs->amd_x_coeff[12]*(5.0*x4 +6.0*x2*y2+y4);
2186 /* magnitude and color terms ignored
2187 wcs->amd_x_coeff[16]*mag + wcs->amd_x_coeff[17]*mag*2.0*xmm +
2188 wcs->amd_x_coeff[18]*mag*(3.0*x2+y2);
2189 */
2190
2191 /* Derivative of X model wrt y */
2192 fy = wcs->amd_x_coeff[1] + wcs->amd_x_coeff[4]*xmm +
2193 wcs->amd_x_coeff[5]*2.0*ymm + wcs->amd_x_coeff[6]*2.0*ymm +
2194 wcs->amd_x_coeff[8]*x2 + wcs->amd_x_coeff[9]*2.0*xy +
2195 wcs->amd_x_coeff[10]*3.0*y2 + wcs->amd_x_coeff[11]*2.0*xy +
2196 wcs->amd_x_coeff[12]*4.0*xy*x2y2;
2197 /* magnitude and color terms ignored
2198 wcs->amd_x_coeff[17]*mag*2.0*ymm +
2199 wcs->amd_x_coeff[18]*mag*2.0*xy;
2200 */
2201
2202 /* Y plate model */
2203 g = wcs->amd_y_coeff[0]*ymm + wcs->amd_y_coeff[1]*xmm +
2204 wcs->amd_y_coeff[2] + wcs->amd_y_coeff[3]*y2 +
2205 wcs->amd_y_coeff[4]*xy + wcs->amd_y_coeff[5]*x2 +
2206 wcs->amd_y_coeff[6]*x2y2 + wcs->amd_y_coeff[7]*y3 +
2207 wcs->amd_y_coeff[8]*y2x + wcs->amd_y_coeff[9]*x2y +
2208 wcs->amd_y_coeff[10]*x3 + wcs->amd_y_coeff[11]*ymm*x2y2 +
2209 wcs->amd_y_coeff[12]*ymm*cjunk;
2210 /* magnitude and color terms ignored
2211 wcs->amd_y_coeff[13]*mag + wcs->amd_y_coeff[14]*mag*mag +
2212 wcs->amd_y_coeff[15]*mag*mag*mag + wcs->amd_y_coeff[16]*mag*ymm +
2213 wcs->amd_y_coeff[17]*mag*x2y2 +
2214 wcs->amd_y_coeff[18]*mag*ymm*x2y2 + wcs->amd_y_coeff[19]*color;
2215 */
2216
2217 /* Derivative of Y model wrt x */
2218 gx = wcs->amd_y_coeff[1] + wcs->amd_y_coeff[4]*ymm +
2219 wcs->amd_y_coeff[5]*2.0*xmm + wcs->amd_y_coeff[6]*2.0*xmm +
2220 wcs->amd_y_coeff[8]*y2 + wcs->amd_y_coeff[9]*2.0*xy +
2221 wcs->amd_y_coeff[10]*3.0*x2 + wcs->amd_y_coeff[11]*2.0*xy +
2222 wcs->amd_y_coeff[12]*4.0*xy*x2y2;
2223 /* magnitude and color terms ignored
2224 wcs->amd_y_coeff[17]*mag*2.0*xmm +
2225 wcs->amd_y_coeff[18]*mag*ymm*2.0*xmm;
2226 */
2227
2228 /* Derivative of Y model wrt y */
2229 gy = wcs->amd_y_coeff[0] + wcs->amd_y_coeff[3]*2.0*ymm +
2230 wcs->amd_y_coeff[4]*xmm + wcs->amd_y_coeff[6]*2.0*ymm +
2231 wcs->amd_y_coeff[7]*3.0*y2 + wcs->amd_y_coeff[8]*2.0*xy +
2232 wcs->amd_y_coeff[9]*x2 + wcs->amd_y_coeff[11]*(x2+3.0*y2) +
2233 wcs->amd_y_coeff[12]*(5.0*y4 + 6.0*x2*y2 + x4);
2234 /* magnitude and color terms ignored
2235 wcs->amd_y_coeff[16]*mag + wcs->amd_y_coeff[17]*mag*2.0*ymm +
2236 wcs->amd_y_coeff[18]*mag*(x2+3.0*y2);
2237 */
2238
2239 f = f - xi;
2240 g = g - eta;
2241 dx = ((-f * gy) + (g * fy)) / ((fx * gy) - (fy * gx));
2242 dy = ((-g * fx) + (f * gx)) / ((fx * gy) - (fy * gx));
2243 xmm = xmm + dx;
2244 ymm = ymm + dy;
2245 if ((fabs(dx) < tolerance) && (fabs(dy) < tolerance)) break;
2246 }
2247
2248 /* Convert mm from plate center to plate pixels */
2249 x = (wcs->ppo_coeff[2] - xmm*1000.0) / wcs->x_pixel_size;
2250 y = (wcs->ppo_coeff[5] + ymm*1000.0) / wcs->y_pixel_size;
2251
2252 /* Convert from plate pixels to image pixels */
2253 *xpix = x - wcs->x_pixel_offset + 1.0 - 0.5;
2254 *ypix = y - wcs->y_pixel_offset + 1.0 - 0.5;
2255
2256 /* If position is off of the image, return offscale code */
2257
2258 /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
2259
2260 Commented out by D.Berry (Starlink) in order to remove dependancy
2261 on NAXIS1/NAXIS2 keywords >>>>>>>>
2262
2263 if (*xpix < 0.5 || *xpix > wcs->nxpix+0.5)
2264 return -1;
2265 if (*ypix < 0.5 || *ypix > wcs->nypix+0.5)
2266 return -1;
2267
2268 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> */
2269
2270 return 0;
2271 }
2272 /* Mar 6 1995 Original version of this code
2273 May 4 1995 Fix eta cross terms which were all in y
2274 Jun 21 1995 Add inverse routine
2275 Oct 17 1995 Fix inverse routine (degrees -> radians)
2276 Nov 7 1995 Add half pixel to image coordinates to get astrometric
2277 plate coordinates
2278 Feb 26 1996 Fix plate to image pixel conversion error
2279 Feb 18 1997 Modified by D.S. Berry (Starlink) to avoid use of the image
2280 dimensions stored in wcs->nxpix and wcs->nypix.
2281 Sep 5 1997 Modified by R.F. Warren-Smith (Starlink) to make the
2282 platepos and platepix functions static.
2283 */
2284
2285
2286
2287
2288