1 /*
2 *class++
3 *  Name:
4 *     WinMap
5 
6 *  Purpose:
7 *     Map one window on to another by scaling and shifting each axis.
8 
9 *  Constructor Function:
10 c     astWinMap
11 f     AST_WINMAP
12 
13 *  Description:
14 *     A Winmap is a linear Mapping which transforms a rectangular
15 *     window in one coordinate system into a similar window in another
16 *     coordinate system by scaling and shifting each axis (the window
17 *     edges being parallel to the coordinate axes).
18 *
19 *     A WinMap is specified by giving the coordinates of two opposite
20 *     corners (A and B) of the window in both the input and output
21 *     coordinate systems.
22 
23 *  Inheritance:
24 *     The WinMap class inherits from the Mapping class.
25 
26 *  Attributes:
27 *     The WinMap class does not define any new attributes beyond those
28 *     which are applicable to all Mappings.
29 
30 *  Functions:
31 c     The WinMap class does not define any new functions beyond those
32 f     The WinMap class does not define any new routines beyond those
33 *     which are applicable to all Mappings.
34 
35 *  Copyright:
36 *     Copyright (C) 1997-2006 Council for the Central Laboratory of the
37 *     Research Councils
38 
39 *  Licence:
40 *     This program is free software: you can redistribute it and/or
41 *     modify it under the terms of the GNU Lesser General Public
42 *     License as published by the Free Software Foundation, either
43 *     version 3 of the License, or (at your option) any later
44 *     version.
45 *
46 *     This program is distributed in the hope that it will be useful,
47 *     but WITHOUT ANY WARRANTY; without even the implied warranty of
48 *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
49 *     GNU Lesser General Public License for more details.
50 *
51 *     You should have received a copy of the GNU Lesser General
52 *     License along with this program.  If not, see
53 *     <http://www.gnu.org/licenses/>.
54 
55 *  Authors:
56 *     DSB: David Berry (Starlink)
57 *     RFWS: R.F. Warren-Smith (Starlink)
58 
59 *  History:
60 *     23-OCT-1996 (DSB):
61 *        Original version.
62 *     4-MAR-1997 (RFWS):
63 *        Tidied public prologues.
64 *     11-MAR-1997 (DSB):
65 *        Added MapMerge method and associated bits.
66 *     30-JUN-1997 (DSB):
67 *        Bug fixed which caused the MapMerge method to generate a
68 *        segmentation violation.
69 *     24-MAR-1998 (RFWS):
70 *        Improved output format from Dump.
71 *     9-APR-1998 (DSB):
72 *        MapMerge modified to allow merging of WinMaps with ZoomMaps and
73 *        and UnitMaps in parallel.
74 *     4-SEP-1998 (DSB):
75 *        Improved MapMerge so that WinMaps can change places with a wider
76 *        range of PermMaps, allowing them to approach closer to a Mapping
77 *        with which they can merge.
78 *     22-FEB-1999 (DSB):
79 *        Corrected logic of MapMerge method to avoid infinite looping.
80 *     5-MAY-1999 (DSB):
81 *        More corrections to MapMerge: Cleared up errors in the use of the
82 *        supplied invert flags, and corrected logic for deciding which
83 *        neighbouring Mapping to swap with.
84 *     16-JUL-1999 (DSB):
85 *        Fixed memory leaks in WinMat and MapMerge.
86 *     8-JAN-2003 (DSB):
87 *        Changed private InitVtab method to protected astInitWinMapVtab
88 *        method.
89 *     8-SEP-2003 (DSB):
90 *        Allow WinMaps to swap with WcsMaps if possible.
91 *     10-NOV-2003 (DSB):
92 *        Modified functions which swap a WinMap with another Mapping
93 *        (e.g. WinPerm, etc), to simplify the returned Mappings.
94 *     23-APR-2004 (DSB):
95 *        Changes to simplification algorithm.
96 *     1-SEP-2004 (DSB):
97 *        Ensure do1 and do2 are initialised before use in MapMerge.
98 *     7-SEP-2005 (DSB):
99 *        Take account of the Invert flag when using the soom factor from
100 *        a ZoomMap.
101 *     14-FEB-2006 (DSB):
102 *        Override astGetObjSize.
103 *     15-MAR-2006 (DSB):
104 *        Override astEqual.
105 *     23-AUG-2006 (DSB):
106 *        Correct initialisation of "result" in the Equal function.
107 *     19-JAN-2007 (DSB):
108 *        Fix memory leak.
109 *     3-MAY-2013 (DSB):
110 *        Improve simplification by adding check for inverse pairs of
111 *        WinMaps in function WinWin.
112 *class--
113 */
114 
115 /* Module Macros. */
116 /* ============== */
117 /* Set the name of the class we are implementing. This indicates to
118    the header files that define class interfaces that they should make
119    "protected" symbols available. */
120 #define astCLASS WinMap
121 
122 /* Include files. */
123 /* ============== */
124 /* Interface definitions. */
125 /* ---------------------- */
126 
127 #include "globals.h"             /* Thread-safe global data access */
128 #include "error.h"               /* Error reporting facilities */
129 #include "memory.h"              /* Memory management facilities */
130 #include "object.h"              /* Base Object class */
131 #include "pointset.h"            /* Sets of points/coordinates */
132 #include "matrixmap.h"           /* Linear mappings */
133 #include "unitmap.h"             /* Unit mappings */
134 #include "zoommap.h"             /* Zoom mappings */
135 #include "permmap.h"             /* Axis permutations */
136 #include "cmpmap.h"              /* Compound mappings */
137 #include "wcsmap.h"              /* Celestial projections */
138 #include "mapping.h"             /* Coordinate mappings (parent class) */
139 #include "channel.h"             /* I/O channels */
140 #include "winmap.h"              /* Interface definition for this class */
141 
142 /* Error code definitions. */
143 /* ----------------------- */
144 #include "ast_err.h"             /* AST error codes */
145 
146 /* C header files. */
147 /* --------------- */
148 #include <float.h>
149 #include <math.h>
150 #include <stdarg.h>
151 #include <stddef.h>
152 #include <stdio.h>
153 #include <string.h>
154 
155 /* Module Variables. */
156 /* ================= */
157 
158 /* Address of this static variable is used as a unique identifier for
159    member of this class. */
160 static int class_check;
161 
162 /* Pointers to parent class methods which are extended by this class. */
163 static int (* parent_getobjsize)( AstObject *, int * );
164 static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
165 static const char *(* parent_getattrib)( AstObject *, const char *, int * );
166 static int (* parent_testattrib)( AstObject *, const char *, int * );
167 static void (* parent_clearattrib)( AstObject *, const char *, int * );
168 static void (* parent_setattrib)( AstObject *, const char *, int * );
169 
170 
171 #ifdef THREAD_SAFE
172 /* Define how to initialise thread-specific globals. */
173 #define GLOBAL_inits \
174    globals->Class_Init = 0;
175 
176 /* Create the function that initialises global data for this module. */
177 astMAKE_INITGLOBALS(WinMap)
178 
179 /* Define macros for accessing each item of thread specific global data. */
180 #define class_init astGLOBAL(WinMap,Class_Init)
181 #define class_vtab astGLOBAL(WinMap,Class_Vtab)
182 
183 
184 #include <pthread.h>
185 
186 
187 #else
188 
189 
190 /* Define the class virtual function table and its initialisation flag
191    as static variables. */
192 static AstWinMapVtab class_vtab;   /* Virtual function table */
193 static int class_init = 0;       /* Virtual function table initialised? */
194 
195 #endif
196 
197 /* External Interface Function Prototypes. */
198 /* ======================================= */
199 /* The following functions have public prototypes only (i.e. no
200    protected prototypes), so we must provide local prototypes for use
201    within this module. */
202 AstWinMap *astWinMapId_( int, const double [], const double [],
203                          const double [], const double [], const char *, ... );
204 
205 /* Prototypes for Private Member Functions. */
206 /* ======================================== */
207 
208 static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * );
209 static AstWinMap *WinUnit( AstWinMap *, AstUnitMap *, int, int, int * );
210 static AstWinMap *WinWin( AstMapping *, AstMapping *, int, int, int, int * );
211 static AstWinMap *WinZoom( AstWinMap *, AstZoomMap *, int, int, int, int, int * );
212 static int GetObjSize( AstObject *, int * );
213 static const char *GetAttrib( AstObject *, const char *, int * );
214 static double Rate( AstMapping *, double *, int, int, int * );
215 static int CanSwap( AstMapping *, AstMapping *, int, int, int *, int * );
216 static int Equal( AstObject *, AstObject *, int * );
217 static int GetIsLinear( AstMapping *, int * );
218 static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * );
219 static int TestAttrib( AstObject *, const char *, int * );
220 static int WinTerms( AstWinMap *, double **, double **, int * );
221 static void ClearAttrib( AstObject *, const char *, int * );
222 static void Copy( const AstObject *, AstObject *, int * );
223 static void Delete( AstObject *, int * );
224 static void Dump( AstObject *, AstChannel *, int * );
225 static void PermGet( AstPermMap *, int **, int **, double **, int * );
226 static void SetAttrib( AstObject *, const char *, int * );
227 static void WinMat( AstMapping **, int *, int, int * );
228 static void WinPerm( AstMapping **, int *, int, int * );
229 static void WinWcs( AstMapping **, int *, int, int * );
230 static int *MapSplit( AstMapping *, int, const int *, AstMapping **, int * );
231 
232 /* Function Macros */
233 /* =============== */
234 /* Macros which return the maximum and minimum of two values. */
235 #define MAX(aa,bb) ((aa)>(bb)?(aa):(bb))
236 #define MIN(aa,bb) ((aa)<(bb)?(aa):(bb))
237 
238 /* Macro to check for equality of floating point values. We cannot
239 compare bad values directory because of the danger of floating point
240 exceptions, so bad values are dealt with explicitly. */
241 #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))))
242 
243 /* Member functions. */
244 /* ================= */
CanSwap(AstMapping * map1,AstMapping * map2,int inv1,int inv2,int * simpler,int * status)245 static int CanSwap( AstMapping *map1, AstMapping *map2, int inv1, int inv2,
246                     int *simpler, int *status ){
247 /*
248 *  Name:
249 *     CanSwap
250 
251 *  Purpose:
252 *     Determine if two Mappings could be swapped.
253 
254 *  Type:
255 *     Private function.
256 
257 *  Synopsis:
258 *     #include "winmap.h"
259 *     int CanSwap( AstMapping *map1, AstMapping *map2, int inv1, int inv2,
260 *                  int *simpler, int *status )
261 
262 *  Class Membership:
263 *     WinMap member function
264 
265 *  Description:
266 *     This function returns a flag indicating if the pair of supplied
267 *     Mappings could be replaced by an equivalent pair of Mappings from the
268 *     same classes as the supplied pair, but in reversed order. Each pair
269 *     of Mappings is considered to be compunded in series. The supplied
270 *     Mapings are not changed in any way.
271 
272 *  Parameters:
273 *     map1
274 *        The Mapping to be applied first.
275 *     map2
276 *        The Mapping to be applied second.
277 *     inv1
278 *        The invert flag to use with map1. A value of zero causes the forward
279 *        mapping to be used, and a non-zero value causes the inverse
280 *        mapping to be used.
281 *     inv2
282 *        The invert flag to use with map2.
283 *     simpler
284 *        Addresss of a location at which to return a flag indicating if
285 *        the swapped Mappings would be intrinsically simpler than the
286 *        original Mappings.
287 *     status
288 *        Pointer to the inherited status variable.
289 
290 *  Returned Value:
291 *     1 if the Mappings could be swapped, 0 otherwise.
292 
293 *  Notes:
294 *     -  One of the supplied pair of Mappings must be a WinMap.
295 *     -  A value of 0 is returned if the two Mappings could be merged into
296 *     a single Mapping.
297 *     -  A value of 0 is returned if an error has already occurred, or if
298 *     this function should fail for any reason.
299 */
300 
301 /* Local Variables: */
302    AstMapping *nowin;        /* Pointer to non-WinMap Mapping */
303    AstWinMap *win;           /* Pointer to the WinMap */
304    const char *class1;       /* Pointer to map1 class string */
305    const char *class2;       /* Pointer to map2 class string */
306    const char *nowin_class;  /* Pointer to non-WinMap class string */
307    double *consts;           /* Pointer to constants array */
308    int *inperm;              /* Pointer to input axis permutation array */
309    int *outperm;             /* Pointer to output axis permutation array */
310    int axlat;                /* Latitude axis in WcsMap */
311    int axlon;                /* Longitude axis in WcsMap */
312    int i;                    /* Loop count */
313    int invert[ 2 ];          /* Original invert flags */
314    int nin;                  /* No. of input coordinates for the PermMap */
315    int nout;                 /* No. of output coordinates for the PermMap */
316    int ret;                  /* Returned flag */
317 
318 /* Check the global error status. */
319    if ( !astOK ) return 0;
320 
321 /* Initialise */
322    ret = 0;
323    *simpler = 0;
324 
325 /* Temporarily set the Invert attributes of both Mappings to the supplied
326    values. */
327    invert[ 0 ] = astGetInvert( map1 );
328    astSetInvert( map1, inv1 );
329 
330    invert[ 1 ] = astGetInvert( map2 );
331    astSetInvert( map2, inv2 );
332 
333 /* Get the classes of the two mappings. */
334    class1 = astGetClass( map1 );
335    class2 = astGetClass( map2 );
336    if( astOK ){
337 
338 /* Get a pointer to the non-WinMap Mapping. */
339       if( !strcmp( class1, "WinMap" ) ){
340          nowin = map2;
341          nowin_class = class2;
342          win = (AstWinMap *) map1;
343       } else {
344          nowin = map1;
345          nowin_class = class1;
346          win = (AstWinMap *) map2;
347       }
348 
349 /* If it is a MatrixMap, the Mappings can be swapped. */
350       if( !strcmp( nowin_class, "MatrixMap" ) ){
351          ret = 1;
352 
353 /* If it is a WcsMap, the Mappings can be swapped if the WinMap is
354    equivalent to a unit transformation on the celestial axes of the
355    WcsMap. */
356       } else if( !strcmp( nowin_class, "WcsMap" ) ){
357 
358 /* Get the indices of the celestial coordinates inthe WcsMap. */
359          axlat = astGetWcsAxis( (AstWcsMap *) nowin, 1 );
360          axlon = astGetWcsAxis( (AstWcsMap *) nowin, 0 );
361 
362 /* Check the shift and scale for these axes. */
363          ret = ( win->a[ axlon ] == 0.0 && win->b[ axlon ] == 1.0 &&
364                  win->a[ axlat ] == 0.0 && win->b[ axlat ] == 1.0 );
365 
366 /* If it is a PermMap, the Mappings can be swapped so long as all links
367    between input and output axes in the PermMap are bi-directional. This
368    does not preclude the existence of unconnected axes, which do not
369    have links (bi-directional or otherwise). */
370       } else if( !strcmp( nowin_class, "PermMap" ) ){
371 
372 /* Get the number of input and output coordinates. */
373          nin = astGetNin( nowin );
374          nout = astGetNout( nowin );
375 
376 /* We need to know the axis permutation arrays and constants array for
377    the PermMap. */
378          PermGet( (AstPermMap *) nowin, &outperm, &inperm, &consts, status );
379          if( astOK ) {
380 
381 /* Indicate we can swap with the PermMap. */
382             ret = 1;
383 
384 /* Check each output axis. If any links between axes are found which are
385    not bi-directional, indicate that we cannot swap with the PermMap. */
386             for( i = 0; i < nout; i++ ){
387                if( outperm[ i ] >= 0 && outperm[ i ] < nin ) {
388                   if( inperm[ outperm[ i ] ] != i ) {
389                      ret = 0;
390                      break;
391                   }
392                }
393             }
394 
395 /* Check each input axis. If any links between axes are found which are
396    not bi-directional, indicate that we cannot swap with the PermMap. */
397             for( i = 0; i < nin; i++ ){
398                if( inperm[ i ] >= 0 && inperm[ i ] < nout ) {
399                   if( outperm[ inperm[ i ] ] != i ) {
400                      ret = 0;
401                      break;
402                   }
403                }
404             }
405 
406 /* If we can swap with the PermMap, the swapped Mappings may be
407    intrinsically simpler than the original mappings. */
408             if( ret ) {
409 
410 /* If the PermMap precedes the WinMap, this will be the case if the PermMap
411    has more outputs than inputs. If the WinMap precedes the PermMap, this
412    will be the case if the PermMap has more inputs than outputs. */
413                *simpler = ( nowin == map1 ) ? nout > nin : nin > nout;
414             }
415 
416 /* Free the axis permutation and constants arrays. */
417             outperm = (int *) astFree( (void *) outperm );
418             inperm = (int *) astFree( (void *) inperm );
419             consts = (double *) astFree( (void *) consts );
420          }
421       }
422    }
423 
424 /* Re-instate the original settings of the Invert attributes for the
425    supplied MatrixMaps. */
426    astSetInvert( map1, invert[ 0 ] );
427    astSetInvert( map2, invert[ 1 ] );
428 
429 /* Return the answer. */
430    return astOK ? ret : 0;
431 }
432 
ClearAttrib(AstObject * this_object,const char * attrib,int * status)433 static void ClearAttrib( AstObject *this_object, const char *attrib, int *status ) {
434 /*
435 *  Name:
436 *     ClearAttrib
437 
438 *  Purpose:
439 *     Clear an attribute value for a WinMap.
440 
441 *  Type:
442 *     Private function.
443 
444 *  Synopsis:
445 *     #include "winmap.h"
446 *     void ClearAttrib( AstObject *this, const char *attrib, int *status )
447 
448 *  Class Membership:
449 *     WinMap member function (over-rides the astClearAttrib protected
450 *     method inherited from the Mapping class).
451 
452 *  Description:
453 *     This function clears the value of a specified attribute for a
454 *     WinMap, so that the default value will subsequently be used.
455 
456 *  Parameters:
457 *     this
458 *        Pointer to the WinMap.
459 *     attrib
460 *        Pointer to a null-terminated string specifying the attribute
461 *        name.  This should be in lower case with no surrounding white
462 *        space.
463 *     status
464 *        Pointer to the inherited status variable.
465 */
466 
467 /* Local Variables: */
468    AstWinMap *this;             /* Pointer to the WinMap structure */
469 
470 /* Check the global error status. */
471    if ( !astOK ) return;
472 
473 /* Obtain a pointer to the WinMap structure. */
474    this = (AstWinMap *) this_object;
475 
476 /* At the moment the WinMap class has no attributes, so pass it on to the
477    parent method for further interpretation. */
478    (*parent_clearattrib)( this_object, attrib, status );
479 
480 }
481 
Equal(AstObject * this_object,AstObject * that_object,int * status)482 static int Equal( AstObject *this_object, AstObject *that_object, int *status ) {
483 /*
484 *  Name:
485 *     Equal
486 
487 *  Purpose:
488 *     Test if two WinMaps are equivalent.
489 
490 *  Type:
491 *     Private function.
492 
493 *  Synopsis:
494 *     #include "winmap.h"
495 *     int Equal( AstObject *this, AstObject *that, int *status )
496 
497 *  Class Membership:
498 *     WinMap member function (over-rides the astEqual protected
499 *     method inherited from the astMapping class).
500 
501 *  Description:
502 *     This function returns a boolean result (0 or 1) to indicate whether
503 *     two WinMaps are equivalent.
504 
505 *  Parameters:
506 *     this
507 *        Pointer to the first Object (a WinMap).
508 *     that
509 *        Pointer to the second Object.
510 *     status
511 *        Pointer to the inherited status variable.
512 
513 *  Returned Value:
514 *     One if the WinMaps are equivalent, zero otherwise.
515 
516 *  Notes:
517 *     - A value of zero will be returned if this function is invoked
518 *     with the global status set, or if it should fail for any reason.
519 */
520 
521 /* Local Variables: */
522    AstWinMap *that;
523    AstWinMap *this;
524    double *a_that;
525    double *a_this;
526    double *b_that;
527    double *b_this;
528    int i;
529    int nin;
530    int result;
531 
532 /* Initialise. */
533    result = 0;
534 
535 /* Check the global error status. */
536    if ( !astOK ) return result;
537 
538 /* Obtain pointers to the two WinMap structures. */
539    this = (AstWinMap *) this_object;
540    that = (AstWinMap *) that_object;
541 
542 /* Check the second object is a WinMap. We know the first is a
543    WinMap since we have arrived at this implementation of the virtual
544    function. */
545    if( astIsAWinMap( that ) ) {
546 
547 /* Get the number of inputs and outputs and check they are the same for both. */
548       nin = astGetNin( this );
549       if( astGetNin( that ) == nin ) {
550 
551 /* Assume the WinMaps are equivalent. */
552          result = 1;
553 
554 /* Compare the shift and scale terms from both WinMaps ignoring the
555    setting of the Invert flag for the moment. */
556          for( i = 0; i < nin; i++ ) {
557             if( !EQUAL( this->a[ i ], that->a[ i ] ) ||
558                 !EQUAL( this->b[ i ], that->b[ i ] ) ) {
559                result = 0;
560                break;
561             }
562          }
563 
564 /* If the scale and shifts are equal, check the Invert flags are equal. */
565          if( result ) {
566             result= ( astGetInvert( this ) == astGetInvert( that ) );
567 
568 /* If the scale and shifts differ, there is still a chance that the
569    WinMaps may be equivalent if their Invert flags differ. */
570          } else if( astGetInvert( this ) != astGetInvert( that ) ) {
571 
572 /* Create copies of the scale and shift terms from the two WinMaps, taking
573    into account the setting of the Invert attribute. Finding the inverted
574    terms involves arithmetic which introduces rounding errors, so this
575    test is not as reliable as the above direct comparison of terms. */
576             astWinTerms( this, &a_this, &b_this );
577             astWinTerms( that, &a_that, &b_that );
578             result = 1;
579 
580             for( i = 0; i < nin; i++ ) {
581                if( !EQUAL( a_this[ i ], a_that[ i ] ) ||
582                    !EQUAL( b_this[ i ], b_that[ i ] ) ) {
583                   result = 0;
584                   break;
585                }
586             }
587 
588 /* Free resources */
589             a_this = astFree( a_this );
590             a_that = astFree( a_that );
591             b_this = astFree( b_this );
592             b_that = astFree( b_that );
593          }
594       }
595    }
596 
597 /* If an error occurred, clear the result value. */
598    if ( !astOK ) result = 0;
599 
600 /* Return the result, */
601    return result;
602 }
603 
GetIsLinear(AstMapping * this_mapping,int * status)604 static int GetIsLinear( AstMapping *this_mapping, int *status ){
605 /*
606 *  Name:
607 *     GetIsLinear
608 
609 *  Purpose:
610 *     Return the value of the IsLinear attribute for a WinMap.
611 
612 *  Type:
613 *     Private function.
614 
615 *  Synopsis:
616 *     #include "mapping.h"
617 *     void GetIsLinear( AstMapping *this, int *status )
618 
619 *  Class Membership:
620 *     WinMap member function (over-rides the protected astGetIsLinear
621 *     method inherited from the Mapping class).
622 
623 *  Description:
624 *     This function returns the value of the IsLinear attribute for a
625 *     Frame, which is always one.
626 
627 *  Parameters:
628 *     this
629 *        Pointer to the WinMap.
630 *     status
631 *        Pointer to the inherited status variable.
632 */
633    return 1;
634 }
635 
GetObjSize(AstObject * this_object,int * status)636 static int GetObjSize( AstObject *this_object, int *status ) {
637 /*
638 *  Name:
639 *     GetObjSize
640 
641 *  Purpose:
642 *     Return the in-memory size of an Object.
643 
644 *  Type:
645 *     Private function.
646 
647 *  Synopsis:
648 *     #include "winmap.h"
649 *     int GetObjSize( AstObject *this, int *status )
650 
651 *  Class Membership:
652 *     WinMap member function (over-rides the astGetObjSize protected
653 *     method inherited from the parent class).
654 
655 *  Description:
656 *     This function returns the in-memory size of the supplied WinMap,
657 *     in bytes.
658 
659 *  Parameters:
660 *     this
661 *        Pointer to the WinMap.
662 *     status
663 *        Pointer to the inherited status variable.
664 
665 *  Returned Value:
666 *     The Object size, in bytes.
667 
668 *  Notes:
669 *     - A value of zero will be returned if this function is invoked
670 *     with the global status set, or if it should fail for any reason.
671 */
672 
673 /* Local Variables: */
674    AstWinMap *this;         /* Pointer to WinMap structure */
675    int result;                /* Result value to return */
676 
677 /* Initialise. */
678    result = 0;
679 
680 /* Check the global error status. */
681    if ( !astOK ) return result;
682 
683 /* Obtain a pointers to the WinMap structure. */
684    this = (AstWinMap *) this_object;
685 
686 /* Invoke the GetObjSize method inherited from the parent class, and then
687    add on any components of the class structure defined by thsi class
688    which are stored in dynamically allocated memory. */
689    result = (*parent_getobjsize)( this_object, status );
690    result += astTSizeOf( this->a );
691    result += astTSizeOf( this->b );
692 
693 /* If an error occurred, clear the result value. */
694    if ( !astOK ) result = 0;
695 
696 /* Return the result, */
697    return result;
698 }
699 
GetAttrib(AstObject * this_object,const char * attrib,int * status)700 static const char *GetAttrib( AstObject *this_object, const char *attrib, int *status ) {
701 /*
702 *  Name:
703 *     GetAttrib
704 
705 *  Purpose:
706 *     Get the value of a specified attribute for a WinMap.
707 
708 *  Type:
709 *     Private function.
710 
711 *  Synopsis:
712 *     #include "winmap.h"
713 *     const char *GetAttrib( AstObject *this, const char *attrib, int *status )
714 
715 *  Class Membership:
716 *     WinMap member function (over-rides the protected astGetAttrib
717 *     method inherited from the Mapping class).
718 
719 *  Description:
720 *     This function returns a pointer to the value of a specified
721 *     attribute for a WinMap, formatted as a character string.
722 
723 *  Parameters:
724 *     this
725 *        Pointer to the WinMap.
726 *     attrib
727 *        Pointer to a null-terminated string containing the name of
728 *        the attribute whose value is required. This name should be in
729 *        lower case, with all white space removed.
730 *     status
731 *        Pointer to the inherited status variable.
732 
733 *  Returned Value:
734 *     - Pointer to a null-terminated string containing the attribute
735 *     value.
736 
737 *  Notes:
738 *     - The returned string pointer may point at memory allocated
739 *     within the WinMap, or at static memory. The contents of the
740 *     string may be over-written or the pointer may become invalid
741 *     following a further invocation of the same function or any
742 *     modification of the WinMap. A copy of the string should
743 *     therefore be made if necessary.
744 *     - A NULL pointer will be returned if this function is invoked
745 *     with the global error status set, or if it should fail for any
746 *     reason.
747 */
748 
749 /* Local Constants: */
750 #define BUFF_LEN 50              /* Max. characters in result buffer */
751 
752 /* Local Variables: */
753    AstWinMap *this;             /* Pointer to the WinMap structure */
754    const char *result;           /* Pointer value to return */
755 
756 /* Initialise. */
757    result = NULL;
758 
759 /* Check the global error status. */
760    if ( !astOK ) return result;
761 
762 /* Obtain a pointer to the WinMap structure. */
763    this = (AstWinMap *) this_object;
764 
765 /* At the moment the WinMap class has no attributes, so pass it on to the
766    parent method for further interpretation. */
767    result = (*parent_getattrib)( this_object, attrib, status );
768 
769 /* Return the result. */
770    return result;
771 
772 /* Undefine macros local to this function. */
773 #undef BUFF_LEN
774 }
775 
astInitWinMapVtab_(AstWinMapVtab * vtab,const char * name,int * status)776 void astInitWinMapVtab_(  AstWinMapVtab *vtab, const char *name, int *status ) {
777 /*
778 *+
779 *  Name:
780 *     astInitWinMapVtab
781 
782 *  Purpose:
783 *     Initialise a virtual function table for a WinMap.
784 
785 *  Type:
786 *     Protected function.
787 
788 *  Synopsis:
789 *     #include "winmap.h"
790 *     void astInitWinMapVtab( AstWinMapVtab *vtab, const char *name )
791 
792 *  Class Membership:
793 *     WinMap vtab initialiser.
794 
795 *  Description:
796 *     This function initialises the component of a virtual function
797 *     table which is used by the WinMap class.
798 
799 *  Parameters:
800 *     vtab
801 *        Pointer to the virtual function table. The components used by
802 *        all ancestral classes will be initialised if they have not already
803 *        been initialised.
804 *     name
805 *        Pointer to a constant null-terminated character string which contains
806 *        the name of the class to which the virtual function table belongs (it
807 *        is this pointer value that will subsequently be returned by the Object
808 *        astClass function).
809 *-
810 */
811 
812 /* Local Variables: */
813    astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
814    AstObjectVtab *object;        /* Pointer to Object component of Vtab */
815    AstMappingVtab *mapping;      /* Pointer to Mapping component of Vtab */
816 
817 /* Check the local error status. */
818    if ( !astOK ) return;
819 
820 /* Get a pointer to the thread specific global data structure. */
821    astGET_GLOBALS(NULL);
822 
823 /* Initialize the component of the virtual function table used by the
824    parent class. */
825    astInitMappingVtab( (AstMappingVtab *) vtab, name );
826 
827 /* Store a unique "magic" value in the virtual function table. This
828    will be used (by astIsAWinMap) to determine if an object belongs
829    to this class.  We can conveniently use the address of the (static)
830    class_check variable to generate this unique value. */
831    vtab->id.check = &class_check;
832    vtab->id.parent = &(((AstMappingVtab *) vtab)->id);
833 
834 /* Initialise member function pointers. */
835 /* ------------------------------------ */
836 /* Store pointers to the member functions (implemented here) that provide
837    virtual methods for this class. */
838    vtab->WinTerms = WinTerms;
839 
840 /* Save the inherited pointers to methods that will be extended, and
841    replace them with pointers to the new member functions. */
842    object = (AstObjectVtab *) vtab;
843    mapping = (AstMappingVtab *) vtab;
844    parent_getobjsize = object->GetObjSize;
845    object->GetObjSize = GetObjSize;
846 
847    parent_clearattrib = object->ClearAttrib;
848    object->ClearAttrib = ClearAttrib;
849    parent_getattrib = object->GetAttrib;
850    object->GetAttrib = GetAttrib;
851    parent_setattrib = object->SetAttrib;
852    object->SetAttrib = SetAttrib;
853    parent_testattrib = object->TestAttrib;
854    object->TestAttrib = TestAttrib;
855 
856    parent_transform = mapping->Transform;
857    mapping->Transform = Transform;
858 
859 /* Store replacement pointers for methods which will be over-ridden by
860    new member functions implemented here. */
861    object->Equal = Equal;
862    mapping->MapMerge = MapMerge;
863    mapping->MapSplit = MapSplit;
864    mapping->Rate = Rate;
865    mapping->GetIsLinear = GetIsLinear;
866 
867 /* Declare the class dump, copy and delete functions.*/
868    astSetDump( vtab, Dump, "WinMap", "Map one window on to another" );
869    astSetCopy( (AstObjectVtab *) vtab, Copy );
870    astSetDelete( (AstObjectVtab *) vtab, Delete );
871 
872 /* If we have just initialised the vtab for the current class, indicate
873    that the vtab is now initialised, and store a pointer to the class
874    identifier in the base "object" level of the vtab. */
875    if( vtab == &class_vtab ) {
876       class_init = 1;
877       astSetVtabClassIdentifier( vtab, &(vtab->id) );
878    }
879 }
880 
MapMerge(AstMapping * this,int where,int series,int * nmap,AstMapping *** map_list,int ** invert_list,int * status)881 static int MapMerge( AstMapping *this, int where, int series, int *nmap,
882                      AstMapping ***map_list, int **invert_list, int *status ) {
883 /*
884 *  Name:
885 *     MapMerge
886 
887 *  Purpose:
888 *     Simplify a sequence of Mappings containing a WinMap.
889 
890 *  Type:
891 *     Private function.
892 
893 *  Synopsis:
894 *     #include "mapping.h"
895 *     int MapMerge( AstMapping *this, int where, int series, int *nmap,
896 *                   AstMapping ***map_list, int **invert_list, int *status )
897 
898 *  Class Membership:
899 *     WinMap method (over-rides the protected astMapMerge method
900 *     inherited from the Mapping class).
901 
902 *  Description:
903 *     This function attempts to simplify a sequence of Mappings by
904 *     merging a nominated WinMap in the sequence with its neighbours,
905 *     so as to shorten the sequence if possible.
906 *
907 *     In many cases, simplification will not be possible and the
908 *     function will return -1 to indicate this, without further
909 *     action.
910 *
911 *     In most cases of interest, however, this function will either
912 *     attempt to replace the nominated WinMap with a Mapping which it
913 *     considers simpler, or to merge it with the Mappings which
914 *     immediately precede it or follow it in the sequence (both will
915 *     normally be considered). This is sufficient to ensure the
916 *     eventual simplification of most Mapping sequences by repeated
917 *     application of this function.
918 *
919 *     In some cases, the function may attempt more elaborate
920 *     simplification, involving any number of other Mappings in the
921 *     sequence. It is not restricted in the type or scope of
922 *     simplification it may perform, but will normally only attempt
923 *     elaborate simplification in cases where a more straightforward
924 *     approach is not adequate.
925 
926 *  Parameters:
927 *     this
928 *        Pointer to the nominated WinMap which is to be merged with
929 *        its neighbours. This should be a cloned copy of the WinMap
930 *        pointer contained in the array element "(*map_list)[where]"
931 *        (see below). This pointer will not be annulled, and the
932 *        WinMap it identifies will not be modified by this function.
933 *     where
934 *        Index in the "*map_list" array (below) at which the pointer
935 *        to the nominated WinMap resides.
936 *     series
937 *        A non-zero value indicates that the sequence of Mappings to
938 *        be simplified will be applied in series (i.e. one after the
939 *        other), whereas a zero value indicates that they will be
940 *        applied in parallel (i.e. on successive sub-sets of the
941 *        input/output coordinates).
942 *     nmap
943 *        Address of an int which counts the number of Mappings in the
944 *        sequence. On entry this should be set to the initial number
945 *        of Mappings. On exit it will be updated to record the number
946 *        of Mappings remaining after simplification.
947 *     map_list
948 *        Address of a pointer to a dynamically allocated array of
949 *        Mapping pointers (produced, for example, by the astMapList
950 *        method) which identifies the sequence of Mappings. On entry,
951 *        the initial sequence of Mappings to be simplified should be
952 *        supplied.
953 *
954 *        On exit, the contents of this array will be modified to
955 *        reflect any simplification carried out. Any form of
956 *        simplification may be performed. This may involve any of: (a)
957 *        removing Mappings by annulling any of the pointers supplied,
958 *        (b) replacing them with pointers to new Mappings, (c)
959 *        inserting additional Mappings and (d) changing their order.
960 *
961 *        The intention is to reduce the number of Mappings in the
962 *        sequence, if possible, and any reduction will be reflected in
963 *        the value of "*nmap" returned. However, simplifications which
964 *        do not reduce the length of the sequence (but improve its
965 *        execution time, for example) may also be performed, and the
966 *        sequence might conceivably increase in length (but normally
967 *        only in order to split up a Mapping into pieces that can be
968 *        more easily merged with their neighbours on subsequent
969 *        invocations of this function).
970 *
971 *        If Mappings are removed from the sequence, any gaps that
972 *        remain will be closed up, by moving subsequent Mapping
973 *        pointers along in the array, so that vacated elements occur
974 *        at the end. If the sequence increases in length, the array
975 *        will be extended (and its pointer updated) if necessary to
976 *        accommodate any new elements.
977 *
978 *        Note that any (or all) of the Mapping pointers supplied in
979 *        this array may be annulled by this function, but the Mappings
980 *        to which they refer are not modified in any way (although
981 *        they may, of course, be deleted if the annulled pointer is
982 *        the final one).
983 *     invert_list
984 *        Address of a pointer to a dynamically allocated array which,
985 *        on entry, should contain values to be assigned to the Invert
986 *        attributes of the Mappings identified in the "*map_list"
987 *        array before they are applied (this array might have been
988 *        produced, for example, by the astMapList method). These
989 *        values will be used by this function instead of the actual
990 *        Invert attributes of the Mappings supplied, which are
991 *        ignored.
992 *
993 *        On exit, the contents of this array will be updated to
994 *        correspond with the possibly modified contents of the
995 *        "*map_list" array.  If the Mapping sequence increases in
996 *        length, the "*invert_list" array will be extended (and its
997 *        pointer updated) if necessary to accommodate any new
998 *        elements.
999 *     status
1000 *        Pointer to the inherited status variable.
1001 
1002 *  Returned Value:
1003 *     If simplification was possible, the function returns the index
1004 *     in the "map_list" array of the first element which was
1005 *     modified. Otherwise, it returns -1 (and makes no changes to the
1006 *     arrays supplied).
1007 
1008 *  Notes:
1009 *     - A value of -1 will be returned if this function is invoked
1010 *     with the global error status set, or if it should fail for any
1011 *     reason.
1012 */
1013 
1014 /* Local Variables: */
1015    AstCmpMap *cm;        /* Pointer to neighbouring CmpMap */
1016    AstMapping **maplt;   /* New mappings list pointer */
1017    AstMapping *map2;     /* Pointer to replacement Mapping */
1018    AstMapping *mc[2];    /* Copies of supplied Mappings to swap */
1019    AstMapping *nc[2];    /* Copies of neighbouring Mappings to merge */
1020    AstMapping *smc0;     /* Simplified Mapping */
1021    AstMapping *smc1;     /* Simplified Mapping */
1022    AstMapping *simp1;    /* Simplified Mapping */
1023    AstMapping *simp2;    /* Simplified Mapping */
1024    AstMatrixMap *mtr;    /* Pointer to replacement MatrixMap */
1025    AstWinMap *newwm2;    /* Second component WinMap */
1026    AstWinMap *newwm;     /* Pointer to replacement WinMap */
1027    AstWinMap *oldwm;     /* Pointer to supplied WinMap */
1028    const char *class1;   /* Pointer to first Mapping class string */
1029    const char *class2;   /* Pointer to second Mapping class string */
1030    const char *nclass;   /* Pointer to neighbouring Mapping class */
1031    double *a;            /* Pointer to zero terms */
1032    double *b;            /* Pointer to scale terms */
1033    int *invlt;           /* New invert flags list pointer */
1034    int cmlow;            /* Is lower neighbour a CmpMap? */
1035    int diag;             /* Is WinMap equivalent to a diagonal matrix? */
1036    int do1;              /* Would a backward swap make a simplification? */
1037    int do2;              /* Would a forward swap make a simplification? */
1038    int i1;               /* Index of first WinMap to merge */
1039    int i2;               /* Index of last WinMap to merge */
1040    int i;                /* Loop counter */
1041    int ic[2];            /* Copies of supplied invert flags to swap */
1042    int inc[4];           /* Copies of supplied invert flags to merge */
1043    int invert;           /* Should the inverted Mapping be used? */
1044    int neighbour;        /* Index of Mapping with which to swap */
1045    int nin2;             /* No. of inputs for second component WinMap */
1046    int nin;              /* Number of coordinates for WinMap */
1047    int nmapt;            /* No. of Mappings in list */
1048    int nstep1;           /* No. of Mappings backwards to next mergable Mapping */
1049    int nstep2;           /* No. of Mappings forward to next mergable Mapping */
1050    int old_winv;         /* original Invert value for supplied WinMap */
1051    int result;           /* Result value to return */
1052    int ser;              /* Are Mappings applied in series? */
1053    int simpler;          /* Is the resulting Mapping simpler than original? */
1054    int swap;             /* Is there an advantage in swapping mappings? */
1055    int swaphi;           /* Can WinMap be swapped with higher neighbour? */
1056    int swaplo;           /* Can WinMap be swapped with lower neighbour? */
1057 
1058 /* Initialise. */
1059    result = -1;
1060 
1061 /* Check the global error status. */
1062    if ( !astOK ) return result;
1063 
1064 /* Initialise variables to avoid "used of uninitialised variable"
1065    messages from dumb compilers. */
1066    i1 = 0;
1067    i2 = 0;
1068    neighbour = 0;
1069 
1070 /* Get the number of axes for the WinMap. */
1071    nin = astGetNin( ( *map_list )[ where ] );
1072 
1073 /* Get a pointer to the WinMap. */
1074    oldwm = (AstWinMap *) this;
1075 
1076 /* First of all, see if the WinMap can be replaced by a simpler Mapping,
1077    without reference to the neighbouring Mappings in the list.           */
1078 /* ======================================================================*/
1079 /* If the shift terms in the WinMap are all zero, the WinMap can be
1080    replaced by a diagonal MatrixMap (which is faster to compute). Check the
1081    shift terms. */
1082    diag = 1;
1083    newwm = (AstWinMap *) ( *map_list )[ where ];
1084    for( i = 0; i < nin; i++ ){
1085       if( !EQUAL( ( newwm->a )[ i ], 0.0 ) ){
1086          diag = 0;
1087          break;
1088       }
1089    }
1090 
1091 /* If all the shift terms are zero... */
1092    if( diag ){
1093 
1094 /* Temporarily set the Invert attribute of the WinMap to the supplied
1095    value. */
1096       old_winv = astGetInvert( newwm );
1097       astSetInvert( newwm, ( *invert_list )[ where ] );
1098 
1099 /* Get a copy of the scale terms from the WinMap. */
1100       astWinTerms( newwm, NULL, &b );
1101 
1102 /* Create a diagonal MatrixMap holding the scale terms. */
1103       mtr = astMatrixMap( nin, nin, 1, b, "", status );
1104 
1105 /* Restore the Invert attribute of the supplied WinMap. */
1106       astSetInvert( newwm, old_winv );
1107 
1108 /* Free the memory used to hold the scale terms. */
1109       b = (double *) astFree( (void *) b );
1110 
1111 /* Annul the WinMap pointer in the list and replace it with the MatrixMap
1112    pointer, and indicate that the forward transformation of the returned
1113    MatrixMap should be used. */
1114       (void) astAnnul( ( *map_list )[ where ] );
1115       ( *map_list )[ where ] = (AstMapping *) mtr;
1116       ( *invert_list )[ where ] = 0;
1117 
1118 /* Return the index of the first modified element. */
1119       result = where;
1120 
1121 /* If the WinMap itself could not be simplified, see if it can be merged
1122    with the Mappings on either side of it in the list. */
1123    } else {
1124 
1125 /* Store the classes of the neighbouring Mappings in the list. */
1126        class1 = ( where > 0 ) ? astGetClass( ( *map_list )[ where - 1 ] ) : NULL;
1127        class2 = ( where < *nmap - 1 ) ? astGetClass( ( *map_list )[ where + 1 ] ) : NULL;
1128 
1129 /* In series. */
1130 /* ========== */
1131       if ( series ) {
1132 
1133 /* We first look to see if the WinMap can be merged with one of its
1134    neighbours, resulting in a reduction of one in the number of Mappings
1135    in the list. WinMaps can only merge directly with another WinMap, a
1136    ZoomMap, or a UnitMap. */
1137          if( class1 && ( !strcmp( class1, "WinMap" ) ||
1138                          !strcmp( class1, "ZoomMap" ) ||
1139                          !strcmp( class1, "UnitMap" ) ) ){
1140             nclass = class1;
1141             i1 = where - 1;
1142             i2 = where;
1143 
1144          } else if( class2 && ( !strcmp( class2, "WinMap" ) ||
1145                                 !strcmp( class2, "ZoomMap" ) ||
1146                                 !strcmp( class2, "UnitMap" ) ) ){
1147             nclass = class2;
1148             i1 = where;
1149             i2 = where + 1;
1150 
1151          } else {
1152             nclass = NULL;
1153          }
1154 
1155 /* If the WinMap can merge with one of its neighbours, create the merged
1156    Mapping. */
1157          if( nclass ){
1158 
1159             if( !strcmp( nclass, "WinMap" ) ){
1160                newwm = WinWin( ( *map_list )[ i1 ], ( *map_list )[ i2 ],
1161                                ( *invert_list )[ i1 ], ( *invert_list )[ i2 ],
1162                                1, status );
1163                invert = 0;
1164 
1165             } else if( !strcmp( nclass, "ZoomMap" ) ){
1166                if( i1 == where ){
1167                   newwm = WinZoom( (AstWinMap *)( *map_list )[ i1 ],
1168                                    (AstZoomMap *)( *map_list )[ i2 ],
1169                            ( *invert_list )[ i1 ], ( *invert_list )[ i2 ], 1, 1, status );
1170                } else {
1171                   newwm = WinZoom( (AstWinMap *)( *map_list )[ i2 ],
1172                                    (AstZoomMap *)( *map_list )[ i1 ],
1173                            ( *invert_list )[ i2 ], ( *invert_list )[ i1 ], 0, 1, status );
1174                }
1175                invert = 0;
1176 
1177             } else {
1178                newwm = astClone( ( *map_list )[ where ] );
1179                invert = ( *invert_list )[ where ];
1180             }
1181 
1182 /* If succesfull... */
1183             if( astOK ){
1184 
1185 /* Annul the first of the two Mappings, and replace it with the merged
1186    WinMap. Also set the invert flag. */
1187                (void) astAnnul( ( *map_list )[ i1 ] );
1188                ( *map_list )[ i1 ] = (AstMapping *) newwm;
1189                ( *invert_list )[ i1 ] = invert;
1190 
1191 /* Annul the second of the two Mappings, and shuffle down the rest of the
1192    list to fill the gap. */
1193                (void) astAnnul( ( *map_list )[ i2 ] );
1194                for ( i = i2 + 1; i < *nmap; i++ ) {
1195                   ( *map_list )[ i - 1 ] = ( *map_list )[ i ];
1196                   ( *invert_list )[ i - 1 ] = ( *invert_list )[ i ];
1197                }
1198 
1199 /* Clear the vacated element at the end. */
1200                ( *map_list )[ *nmap - 1 ] = NULL;
1201                ( *invert_list )[ *nmap - 1 ] = 0;
1202 
1203 /* Decrement the Mapping count and return the index of the first
1204    modified element. */
1205                ( *nmap )--;
1206                result = i1;
1207 
1208             }
1209 
1210 /* If one of the neighbours is a (parallel) CmpMap, we convert the WinMap
1211    into an equivalent parallel CmpMap, and then merge this parallel
1212    CmpMap with the neighbouring parallel CmpMap to create a parallel CmpMap
1213    containing two series CmpMaps. */
1214          } else if( ( class1 && !strcmp( "CmpMap", class1 ) ) ||
1215                     ( class2 && !strcmp( "CmpMap", class2 ) ) ) {
1216 
1217 /* Identify the WinMap and the CmpMap. */
1218             if( class1 && !strcmp( "CmpMap", class1 ) ) {
1219                i1 = where - 1;
1220                i2 = where;
1221                cm = (AstCmpMap *) ( *map_list )[ where - 1 ];
1222                cmlow = 1;
1223 
1224             } else {
1225                i1 = where;
1226                i2 = where + 1;
1227                cm = (AstCmpMap *) ( *map_list )[ where + 1 ];
1228                cmlow = 0;
1229 
1230             }
1231 
1232 /* Temporarily set the required Invert attributes in the two Mappings. */
1233             inc[ 0 ] = astGetInvert( ( *map_list )[ i1 ] );
1234             astSetInvert( ( *map_list )[ i1 ], ( *invert_list )[ i1 ] );
1235 
1236             inc[ 1 ] = astGetInvert( ( *map_list )[ i2 ] );
1237             astSetInvert( ( *map_list )[ i2 ], ( *invert_list )[ i2 ] );
1238 
1239 /* Now get pointers to the scale and zero terms of the nominated WinMap
1240    (these describe the forward transformation, taking into account the
1241    setting of the Invert flag). */
1242             (void) astWinTerms( oldwm , &a, &b );
1243 
1244 /* Get pointers to the two components of the parallel CmpMap. */
1245             astDecompose( cm, mc, mc + 1, &ser, ic, ic + 1 );
1246 
1247 /* Check component Mappings are combined in parallel. */
1248             map2 = NULL;
1249             if( astOK && !ser ) {
1250 
1251 /* Temporarily set the required Invert attributes in the two component
1252    Mappings to the indicated values. */
1253                inc[ 2 ] = astGetInvert( mc[ 0 ] );
1254                astSetInvert( mc[ 0 ], ic[ 0 ] );
1255 
1256                inc[ 3 ] = astGetInvert( mc[ 1 ] );
1257                astSetInvert( mc[ 1 ], ic[ 1 ] );
1258 
1259 /* Create the first of two corresponding WinMaps, initially with undefined
1260    corners. These could be combined into a parallel CmpMap which would be
1261    equivalent to the nominated WinMap. The number of inputs for each WinMap
1262    is equal to either the number of outputs or inputs of the corresponding
1263    component of the CmpMap, depending on whether the CmpMap is upper or lower
1264    neighbour. */
1265                nin = cmlow ? astGetNout( mc[ 0 ] ):astGetNin( mc[ 0 ] );
1266                newwm = astWinMap( nin, NULL, NULL, NULL, NULL, "", status );
1267                if( astOK ) {
1268 
1269 /* Store the first "nin" scale and zero terms from the nominated WinMap
1270    in the new WinMap. */
1271                   for( i = 0; i < nin; i++ ) {
1272                      (newwm->a)[ i ] = a[ i ];
1273                      (newwm->b)[ i ] = b[ i ];
1274                   }
1275                }
1276 
1277 /* Now create the second WinMap in the same way, which transforms the
1278    remaining outputs of the CmpMap. */
1279                nin2 = cmlow ? astGetNout( mc[ 1 ] ):astGetNin( mc[ 1 ] );
1280                newwm2 = astWinMap( nin2, NULL, NULL, NULL, NULL, "", status );
1281                if( astOK ) {
1282 
1283 /* Store the remaining scale and zero terms from the nominated WinMap
1284    in the new WinMap. */
1285                   for( i = 0; i < nin2; i++ ) {
1286                      (newwm2->a)[ i ] = a[ i + nin ];
1287                      (newwm2->b)[ i ] = b[ i + nin ];
1288                   }
1289                }
1290 
1291 /* Combine the two corresponding lower component Mappings into a series
1292    CmpMap, and likewise combine the two corresponding upper component
1293    Mappings into a series CmpMap. */
1294                if( cmlow ) {
1295                   nc[ 0 ] = (AstMapping *) astCmpMap( mc[ 0 ], newwm, 1, "", status );
1296                   nc[ 1 ] = (AstMapping *) astCmpMap( mc[ 1 ], newwm2, 1, "", status );
1297                } else {
1298                   nc[ 0 ] = (AstMapping *) astCmpMap( newwm, mc[ 0 ], 1, "", status );
1299                   nc[ 1 ] = (AstMapping *) astCmpMap( newwm2, mc[ 1 ], 1, "", status );
1300                }
1301                newwm = astAnnul( newwm );
1302                newwm2 = astAnnul( newwm2 );
1303 
1304 /* Attempt to simplify each of the two new series CmpMaps. If neither of
1305    them simplify then there is no point in doing the current merger. In fact
1306    it would be dangerous to do so since we may end up in an infinite loop
1307    where the resulting parallel CmpMap gets converted back into the
1308    existing series CmpMap by the CmpMap MapMerge method, and then back
1309    again by this method, etc. */
1310                simp1 = astSimplify( nc[ 0 ] );
1311                simp2 = astSimplify( nc[ 1 ] );
1312 
1313 /* Test if either could be simplified by checking if its pointer value
1314    has changed. */
1315                simpler = ( simp1 != nc[ 0 ] ) || ( simp2 != nc[ 1 ] );
1316 
1317 /* If either CmpMap was simplified, then combine the two series CmpMap into
1318    a single parallel CmpMap. */
1319                if( simpler ) {
1320                   map2 = (AstMapping *) astCmpMap( simp1, simp2, 0, "", status );
1321                }
1322 
1323 /* Re-instate the original Invert attributes in the two component Mappings. */
1324                astSetInvert( mc[ 0 ], inc[ 2 ] );
1325                astSetInvert( mc[ 1 ], inc[ 3 ] );
1326 
1327 /* Free resources. */
1328                simp1 = astAnnul( simp1 );
1329                simp2 = astAnnul( simp2 );
1330                nc[ 0 ] = astAnnul( nc[ 0 ] );
1331                nc[ 1 ] = astAnnul( nc[ 1 ] );
1332 
1333             }
1334 
1335 /* Free resources. */
1336             mc[ 0 ] = astAnnul( mc[ 0 ] );
1337             mc[ 1 ] = astAnnul( mc[ 1 ] );
1338             a = astFree( a );
1339             b = astFree( b );
1340 
1341 /* Re-instate the original Invert attributes. */
1342             astSetInvert( ( *map_list )[ i1 ], inc[ 0 ] );
1343             astSetInvert( ( *map_list )[ i2 ], inc[ 1 ] );
1344 
1345 /* If the above produced a new Mapping, annul the supplied pointers for
1346    the two merged Mappings, store the pointer for the new merged Mapping,
1347    and shuffle the remaining Mappings down to fill the space left. Nullify
1348    the end slot which is no longer used, reduce the number of Mappings in
1349    the list by 1, and return the index of the first modified Mapping. */
1350             if( map2 ) {
1351                (void) astAnnul( ( *map_list )[ i1 ] );
1352                (void) astAnnul( ( *map_list )[ i2 ] );
1353                ( *map_list )[ i1 ] = map2;
1354                ( *invert_list )[ i1 ] = 0;
1355                for( i = i2 + 1; i < *nmap; i++ ){
1356                   ( *map_list )[ i - 1 ] = ( *map_list )[ i ];
1357                   ( *invert_list )[ i - 1 ] = ( *invert_list )[ i ];
1358                }
1359                ( *map_list )[ *nmap - 1 ] = NULL;
1360                (*nmap)--;
1361                result = i1;
1362             }
1363 
1364 /* If the WinMap could not merge directly with either of its neighbours,
1365    we consider whether it would be worthwhile to swap the WinMap with
1366    either of its neighbours. This can only be done for certain classes
1367    of Mapping (MatrixMap & some PermMaps & WcsMaps), and will usually require both
1368    Mappings to be modified (unless they are commutative). The advantage of
1369    swapping the order of the Mappings is that it may result in the WinMap
1370    being adjacent to a Mapping with which it can merge directly on the next
1371    invocation of this function, thus reducing the number of Mappings
1372    in the list. */
1373          } else {
1374 
1375 /* Set a flag if we could swap the WinMap with its higher neighbour. "do2"
1376    is returned if swapping the Mappings would simplify either of the
1377    Mappings. */
1378             if( where + 1 < *nmap ){
1379                swaphi = CanSwap(  ( *map_list )[ where ],
1380                                   ( *map_list )[ where + 1 ],
1381                                   ( *invert_list )[ where ],
1382                                   ( *invert_list )[ where + 1 ], &do2, status );
1383             } else {
1384                swaphi = 0;
1385                do2 = 0;
1386             }
1387 
1388 /* If so, step through each of the Mappings which follow the WinMap,
1389    looking for a Mapping with which the WinMap could merge directly. Stop
1390    when such a Mapping is found, or if a Mapping is found with which the
1391    WinMap could definitely not swap. Note the number of Mappings which
1392    separate the WinMap from the Mapping with which it could merge (if
1393    any). */
1394             nstep2 = -1;
1395             if( swaphi ){
1396                for( i2 = where + 1; i2 < *nmap; i2++ ){
1397 
1398 /* See if we can merge with this Mapping. If so, note the number of steps
1399    between the two Mappings and leave the loop. */
1400                   nclass = astGetClass( ( *map_list )[ i2 ] );
1401                   if( !strcmp( nclass, "WinMap" ) ||
1402                       !strcmp( nclass, "ZoomMap" ) ||
1403                       !strcmp( nclass, "UnitMap" ) ) {
1404                      nstep2 = i2 - where - 1;
1405                      break;
1406                   }
1407 
1408 /* If there is no chance that we can swap with this Mapping, leave the loop
1409    with -1 for the number of steps to indicate that no merging is possible.
1410    WinMaps can swap with MatrixMaps and some PermMaps. */
1411                   if( strcmp( nclass, "MatrixMap" ) &&
1412                       strcmp( nclass, "WcsMap" ) &&
1413                       strcmp( nclass, "PermMap" ) ) {
1414                      break;
1415                   }
1416 
1417                }
1418 
1419             }
1420 
1421 /* Do the same working forward from the WinMap towards the start of the map
1422    list. */
1423             if( where > 0 ){
1424                swaplo = CanSwap(  ( *map_list )[ where - 1 ],
1425                                   ( *map_list )[ where ],
1426                                   ( *invert_list )[ where - 1 ],
1427                                   ( *invert_list )[ where ], &do1, status );
1428             } else {
1429                swaplo = 0;
1430                do1 = 0;
1431             }
1432 
1433             nstep1 = -1;
1434             if( swaplo ){
1435                for( i1 = where - 1; i1 >= 0; i1-- ){
1436 
1437                   nclass = astGetClass( ( *map_list )[ i1 ] );
1438                   if( !strcmp( nclass, "WinMap" ) ||
1439                       !strcmp( nclass, "ZoomMap" ) ||
1440                       !strcmp( nclass, "UnitMap" ) ) {
1441                      nstep1 = where - 1 - i1;
1442                      break;
1443                   }
1444 
1445                   if( strcmp( nclass, "MatrixMap" ) &&
1446                       strcmp( nclass, "WcsMap" ) &&
1447                       strcmp( nclass, "PermMap" ) ) {
1448                      break;
1449                   }
1450 
1451                }
1452 
1453             }
1454 
1455 /* Choose which neighbour to swap with so that the WinMap moves towards the
1456    nearest Mapping with which it can merge. */
1457             if( do1 || (
1458                 nstep1 != -1 && ( nstep2 == -1 || nstep2 > nstep1 ) ) ){
1459                nclass = class1;
1460                i1 = where - 1;
1461                i2 = where;
1462                neighbour = i1;
1463             } else if( do2 || nstep2 != -1 ){
1464                nclass = class2;
1465                i1 = where;
1466                i2 = where + 1;
1467                neighbour = i2;
1468             } else {
1469                nclass = NULL;
1470             }
1471 
1472 /* If there is a target Mapping in the list with which the WinMap could
1473    merge, replace the supplied Mappings with swapped Mappings to bring a
1474    WinMap closer to the target Mapping. */
1475             if( nclass ){
1476 
1477 /* It is possible that the neighbouring Mapping with which we are about to
1478    swap could also merge with the target Mapping. When the neighbouring
1479    Mapping is reconsidered it may well swap the pair back to put itself nearer
1480    the target Mapping. We need to be careful not to end up in an infinite loop
1481    in which the pair of neighbouring Mappings are constantly swapped backwards
1482    and forwards as each attempts to put itself closer to the target Mapping.
1483    To prevent this, we only swap the pair of Mappings if the neighbouring
1484    Mapping could not itself merge with the target Mapping. Check to see
1485    if this is the case by attempting to merge the neighbouring Mapping with
1486    the target Mapping. */
1487                map2 = astClone( (*map_list)[ neighbour ] );
1488                nmapt = *nmap - neighbour;
1489                maplt = *map_list + neighbour;
1490                invlt = *invert_list + neighbour;
1491                result = astMapMerge( map2, 0, series, &nmapt, &maplt, &invlt );
1492                map2 = astAnnul( map2 );
1493 
1494 /* If the above call produced a change in the  Mapping list, return the
1495    remaining number of mappings.. */
1496                if( result != -1 ){
1497                   *nmap = nmapt + neighbour;
1498 
1499 /* Otherwise, if there was no change in the mapping list... */
1500                } else {
1501 
1502                   if( !strcmp( nclass, "MatrixMap" ) ){
1503                      WinMat( (*map_list) + i1, (*invert_list) + i1, where - i1, status );
1504 
1505                   } else if( !strcmp( nclass, "PermMap" ) ){
1506                      WinPerm( (*map_list) + i1, (*invert_list) + i1, where - i1, status );
1507 
1508                   } else if( !strcmp( nclass, "WcsMap" ) ){
1509                      WinWcs( (*map_list) + i1, (*invert_list) + i1, where - i1, status );
1510                   }
1511 
1512 /* Store the index of the first modified Mapping. */
1513                   result = i1;
1514                }
1515 
1516 /* If there is no Mapping available for merging, it may still be
1517    advantageous to swap with a neighbour because the swapped Mapping may
1518    be simpler than the original Mappings. For instance, a PermMap may
1519    strip axes of the WinMap leaving only a UnitMap. Also, the two neighbours
1520    may be able to merge. */
1521             } else if( swaphi || swaplo ) {
1522 
1523 /* Try swapping with each possible neighbour in turn. */
1524                for( i = 0; i < 2; i++ ) {
1525 
1526 /*  Set up the class and pointers for the mappings to be swapped, first
1527     the lower neighbour, then the upper neighbour. */
1528                   if( i == 0 && swaplo ){
1529                      nclass = class1;
1530                      i1 = where - 1;
1531                      i2 = where;
1532 
1533                   } else if( i == 1 && swaphi ){
1534                      nclass = class2;
1535                      i1 = where;
1536                      i2 = where + 1;
1537 
1538                   } else {
1539                      nclass = NULL;
1540                   }
1541 
1542 /* If we have a Mapping to swap with... */
1543                   if( nclass ) {
1544 
1545 /* Take copies of the Mapping and Invert flag arrays so we do not change
1546    the supplied values. */
1547                      mc[ 0 ] = (AstMapping *) astCopy( ( (*map_list) + i1 )[0] );
1548                      mc[ 1 ] = (AstMapping *) astCopy( ( (*map_list) + i1 )[1] );
1549                      ic[ 0 ] = ( (*invert_list) + i1 )[0];
1550                      ic[ 1 ] = ( (*invert_list) + i1 )[1];
1551 
1552 /* Swap these Mappings. */
1553                      if( !strcmp( nclass, "MatrixMap" ) ){
1554                         WinMat( mc, ic, where - i1, status );
1555                      } else if( !strcmp( nclass, "PermMap" ) ){
1556                         WinPerm( mc, ic, where - i1, status );
1557                      } else if( !strcmp( nclass, "WcsMap" ) ){
1558                         WinWcs( mc, ic, where - i1, status );
1559                      }
1560 
1561 /* See if the two neighbouring Mappings can merge now that the nominated
1562    Mapping is no longer in between them. First get a list of Mapping
1563    pointers containing the two Mappings to be merged, and associated
1564    invert flags. */
1565                      if( i == 0 && where != *nmap - 1 ) {
1566                         nc[ 0 ] = astClone( mc[ 1 ] );
1567                         nc[ 1 ] = astClone( (*map_list)[ where + 1 ] );
1568                         inc[ 0 ] = ic[ 1 ];
1569                         inc[ 1 ] = (*invert_list)[ where + 1 ];
1570 
1571                      } else if( i == 1 && where > 0 ) {
1572                         nc[ 0 ] = astClone( (*map_list)[ where - 1 ] );
1573                         nc[ 1 ] = astClone( mc[ 0 ] );
1574                         inc[ 0 ] = (*invert_list)[ where - 1 ];
1575                         inc[ 1 ] = ic[ 0 ];
1576 
1577                      } else {
1578                         nc[ 0 ] = NULL;
1579                         nc[ 1 ] = NULL;
1580                      }
1581 
1582 /* If both neighbours are available, use astMapMerge to see if it is
1583    possible to merge the two Mappings. */
1584                      swap = 0;
1585                      if( nc[ 0 ] && nc[ 1 ] ) {
1586                         nmapt = 2;
1587                         maplt = nc;
1588                         invlt = inc;
1589                         map2 = astClone( nc[ 0 ] );
1590                         swap = astMapMerge( map2, 0, series, &nmapt, &maplt, &invlt );
1591                         map2 = astAnnul( map2 );
1592                         if( swap == -1 ) {
1593                            map2 = astClone( nc[ 1 ] );
1594                            swap = astMapMerge( map2, 1, series, &nmapt, &maplt, &invlt );
1595                            map2 = astAnnul( map2 );
1596                         }
1597                         swap = ( nmapt < 2 ) ? 1 : 0;
1598                      }
1599 
1600 /* Free resources. */
1601                      if( nc[ 0 ] ) nc[ 0 ] = astAnnul( nc[ 0 ] );
1602                      if( nc[ 1 ] ) nc[ 1 ] = astAnnul( nc[ 1 ] );
1603 
1604 /* If the neighbours could not merge, see if either swapped Mapping can
1605    be simplified. */
1606                      if( !swap ) {
1607                         smc0 = astSimplify( mc[0] );
1608                         if(  smc0 != mc[0] ) {
1609                            swap = 1;
1610                         } else {
1611                            smc1 = astSimplify( mc[1] );
1612                            swap = ( smc1 != mc[1] );
1613                            smc1 = astAnnul( smc1 );
1614                         }
1615                         smc0 = astAnnul( smc0 );
1616                      }
1617 
1618 /* If there is some point in swapping the Mappings, swap them in the
1619    supplied lists. Otherwise annul the swapped Mappings. */
1620                      if( swap ) {
1621                         (*map_list)[ i1 ] = astAnnul( (*map_list)[ i1 ] );
1622                         (*map_list)[ i2 ] = astAnnul( (*map_list)[ i2 ] );
1623                         (*map_list)[ i1 ] = mc[ 0 ];
1624                         (*map_list)[ i2 ] = mc[ 1 ];
1625                         (*invert_list)[ i1 ] = ic[ 0 ];
1626                         (*invert_list)[ i2 ] = ic[ 1 ];
1627                         result = i1;
1628                         break;
1629 
1630                      } else {
1631                         mc[ 0 ] = astAnnul( mc[ 0 ] );
1632                         mc[ 1 ] = astAnnul( mc[ 1 ] );
1633                      }
1634                   }
1635                }
1636             }
1637          }
1638 
1639 /* In parallel. */
1640 /* ============ */
1641 /* WinMaps are combined in parallel with neighbouring WinMaps, ZoomMaps and
1642    UnitMaps. */
1643       } else {
1644 
1645 /* We first look to see if the WinMap can be merged with one of its
1646    neighbours, resulting in a reduction of one in the number of Mappings
1647    in the list. WinMaps can only merge directly with another WinMap, a
1648    ZoomMap, or a UnitMap. */
1649          if( class1 && ( !strcmp( class1, "WinMap" ) ||
1650                          !strcmp( class1, "ZoomMap" ) ||
1651                          !strcmp( class1, "UnitMap" ) ) ){
1652             nclass = class1;
1653             i1 = where - 1;
1654             i2 = where;
1655 
1656          } else if( class2 && ( !strcmp( class2, "WinMap" ) ||
1657                                 !strcmp( class2, "ZoomMap" ) ||
1658                                 !strcmp( class2, "UnitMap" ) ) ){
1659             nclass = class2;
1660             i1 = where;
1661             i2 = where + 1;
1662 
1663          } else {
1664             nclass = NULL;
1665          }
1666 
1667 /* If the WinMap can merge with one of its neighbours, create the merged
1668    Mapping. */
1669          if( nclass ){
1670 
1671             if( !strcmp( nclass, "WinMap" ) ){
1672                newwm = WinWin( ( *map_list )[ i1 ], ( *map_list )[ i2 ],
1673                                ( *invert_list )[ i1 ], ( *invert_list )[ i2 ],
1674                                0, status );
1675                invert = 0;
1676 
1677             } else if( !strcmp( nclass, "ZoomMap" ) ){
1678                if( i1 == where ){
1679                   newwm = WinZoom( (AstWinMap *)( *map_list )[ i1 ],
1680                                    (AstZoomMap *)( *map_list )[ i2 ],
1681                            ( *invert_list )[ i1 ], ( *invert_list )[ i2 ], 1, 0, status );
1682                } else {
1683                   newwm = WinZoom( (AstWinMap *)( *map_list )[ i2 ],
1684                                    (AstZoomMap *)( *map_list )[ i1 ],
1685                            ( *invert_list )[ i2 ], ( *invert_list )[ i1 ], 0, 0, status );
1686                }
1687                invert = 0;
1688 
1689             } else {
1690                if( i1 == where ){
1691                   newwm = WinUnit( (AstWinMap *)( *map_list )[ i1 ],
1692                                    (AstUnitMap *)( *map_list )[ i2 ],
1693                                    ( *invert_list )[ i1 ], 1, status );
1694                } else {
1695                   newwm = WinUnit( (AstWinMap *)( *map_list )[ i2 ],
1696                                    (AstUnitMap *)( *map_list )[ i1 ],
1697                                    ( *invert_list )[ i2 ], 0, status );
1698                }
1699                invert = 0;
1700 
1701             }
1702 
1703 /* If succesfull... */
1704             if( astOK ){
1705 
1706 /* Annul the first of the two Mappings, and replace it with the merged
1707    WinMap. Also set the invert flag. */
1708                (void) astAnnul( ( *map_list )[ i1 ] );
1709                ( *map_list )[ i1 ] = (AstMapping *) newwm;
1710                ( *invert_list )[ i1 ] = invert;
1711 
1712 /* Annul the second of the two Mappings, and shuffle down the rest of the
1713    list to fill the gap. */
1714                (void) astAnnul( ( *map_list )[ i2 ] );
1715                for ( i = i2 + 1; i < *nmap; i++ ) {
1716                   ( *map_list )[ i - 1 ] = ( *map_list )[ i ];
1717                   ( *invert_list )[ i - 1 ] = ( *invert_list )[ i ];
1718                }
1719 
1720 /* Clear the vacated element at the end. */
1721                ( *map_list )[ *nmap - 1 ] = NULL;
1722                ( *invert_list )[ *nmap - 1 ] = 0;
1723 
1724 /* Decrement the Mapping count and return the index of the first
1725    modified element. */
1726                ( *nmap )--;
1727                result = i1;
1728 
1729             }
1730          }
1731       }
1732    }
1733 
1734 /* Return the result. */
1735    return result;
1736 }
1737 
MapSplit(AstMapping * this_map,int nin,const int * in,AstMapping ** map,int * status)1738 static int *MapSplit( AstMapping *this_map, int nin, const int *in, AstMapping **map, int *status ){
1739 /*
1740 *  Name:
1741 *     MapSplit
1742 
1743 *  Purpose:
1744 *     Create a Mapping representing a subset of the inputs of an existing
1745 *     WinMap.
1746 
1747 *  Type:
1748 *     Private function.
1749 
1750 *  Synopsis:
1751 *     #include "winmap.h"
1752 *     int *MapSplit( AstMapping *this, int nin, const int *in, AstMapping **map, int *status )
1753 
1754 *  Class Membership:
1755 *     WinMap method (over-rides the protected astMapSplit method
1756 *     inherited from the Mapping class).
1757 
1758 *  Description:
1759 *     This function creates a new Mapping by picking specified inputs from
1760 *     an existing WinMap. This is only possible if the specified inputs
1761 *     correspond to some subset of the WinMap outputs. That is, there
1762 *     must exist a subset of the WinMap outputs for which each output
1763 *     depends only on the selected WinMap inputs, and not on any of the
1764 *     inputs which have not been selected. If this condition is not met
1765 *     by the supplied WinMap, then a NULL Mapping is returned.
1766 
1767 *  Parameters:
1768 *     this
1769 *        Pointer to the WinMap to be split (the WinMap is not actually
1770 *        modified by this function).
1771 *     nin
1772 *        The number of inputs to pick from "this".
1773 *     in
1774 *        Pointer to an array of indices (zero based) for the inputs which
1775 *        are to be picked. This array should have "nin" elements. If "Nin"
1776 *        is the number of inputs of the supplied WinMap, then each element
1777 *        should have a value in the range zero to Nin-1.
1778 *     map
1779 *        Address of a location at which to return a pointer to the new
1780 *        Mapping. This Mapping will have "nin" inputs (the number of
1781 *        outputs may be different to "nin"). A NULL pointer will be
1782 *        returned if the supplied WinMap has no subset of outputs which
1783 *        depend only on the selected inputs.
1784 *     status
1785 *        Pointer to the inherited status variable.
1786 
1787 *  Returned Value:
1788 *     A pointer to a dynamically allocated array of ints. The number of
1789 *     elements in this array will equal the number of outputs for the
1790 *     returned Mapping. Each element will hold the index of the
1791 *     corresponding output in the supplied WinMap. The array should be
1792 *     freed using astFree when no longer needed. A NULL pointer will
1793 *     be returned if no output Mapping can be created.
1794 
1795 *  Notes:
1796 *     - If this function is invoked with the global error status set,
1797 *     or if it should fail for any reason, then NULL values will be
1798 *     returned as the function value and for the "map" pointer.
1799 */
1800 
1801 /* Local Variables: */
1802    AstWinMap *newwm;          /* Pointer to returned WinMap */
1803    AstWinMap *this;           /* Pointer to WinMap structure */
1804    double *a;                 /* Pointer to zero terms */
1805    double *b;                 /* Pointer to scale terms */
1806    int *result;               /* Pointer to returned array */
1807    int i;                     /* Loop count */
1808    int iin;                   /* Mapping input index */
1809    int mnin;                  /* No. of Mapping inputs */
1810    int ok;                    /* Are input indices OK? */
1811 
1812 /* Initialise */
1813    result = NULL;
1814    *map = NULL;
1815 
1816 /* Check the global error status. */
1817    if ( !astOK ) return result;
1818 
1819 /* Get a pointer to the WinMap structure. */
1820    this = (AstWinMap *) this_map;
1821 
1822 /* Allocate memory for the returned array and create a WinMap with the
1823    required number of axes and undefined corners. */
1824    result = astMalloc( sizeof( int )*(size_t) nin );
1825    newwm = astWinMap( nin, NULL, NULL, NULL, NULL, "", status );
1826    *map = (AstMapping *) newwm;
1827 
1828 /* Now get pointers to the scale and zero terms of the supplied WinMap
1829    (these describe the forward transformation, taking into account the
1830    setting of the Invert flag). */
1831    (void) astWinTerms( this , &a, &b );
1832 
1833 /* Check pointers can be used safely. */
1834    if( astOK ) {
1835 
1836 /* Store the required scale and zero terms from the supplied WinMap
1837    in the new WinMap. At the same time check that each axis is valid. */
1838       mnin = astGetNin( this );
1839       ok = 1;
1840       for( i = 0; i < nin; i++ ) {
1841          iin = in[ i ];
1842          if( iin >= 0 && iin < mnin ) {
1843             (newwm->a)[ i ] = a[ iin ];
1844             (newwm->b)[ i ] = b[ iin ];
1845             result[ i ] = iin;
1846          } else {
1847             ok = 0;
1848             break;
1849          }
1850       }
1851 
1852 /* If the "in" array contained any invalid values, free the returned
1853    resources. */
1854       if( !ok ) {
1855          result = astFree( result );
1856          *map = astAnnul( *map );
1857       }
1858    }
1859 
1860 /* Free resources. */
1861    a = astFree( a );
1862    b = astFree( b );
1863 
1864 /* Free returned resources if an error has occurred. */
1865    if( !astOK ) {
1866       result = astFree( result );
1867       *map = astAnnul( *map );
1868    }
1869 
1870 /* Return the list of output indices. */
1871    return result;
1872 }
1873 
PermGet(AstPermMap * map,int ** outperm,int ** inperm,double ** consts,int * status)1874 static void PermGet( AstPermMap *map, int **outperm, int **inperm,
1875                      double **consts, int *status ){
1876 /*
1877 *  Name:
1878 *     PermGet
1879 
1880 *  Purpose:
1881 *     Get the axis permutation and constants array for a PermMap.
1882 
1883 *  Type:
1884 *     Private function.
1885 
1886 *  Synopsis:
1887 *     #include "winmap.h"
1888 *     void PermGet( AstPermMap *map, int **outperm, int **inperm,
1889 *                   double **const, int *status )
1890 
1891 *  Class Membership:
1892 *     WinMap member function
1893 
1894 *  Description:
1895 *     This function returns axis permutation and constants arrays which can
1896 *     be used to create a PermMap which is equivalent to the supplied PermMap.
1897 
1898 *  Parameters:
1899 *     map
1900 *        The PermMap.
1901 *     outperm
1902 *        An address at which to return a popinter to an array of ints
1903 *        holding the output axis permutation array. The array should be
1904 *        released using astFree when no longer needed.
1905 *     inperm
1906 *        An address at which to return a popinter to an array of ints
1907 *        holding the input axis permutation array. The array should be
1908 *        released using astFree when no longer needed.
1909 *     consts
1910 *        An address at which to return a popinter to an array of doubles
1911 *        holding the constants array. The array should be released using
1912 *        astFree when no longer needed.
1913 *     status
1914 *        Pointer to the inherited status variable.
1915 
1916 *  Notes:
1917 *     -  NULL pointers are returned if an error has already occurred, or if
1918 *     this function should fail for any reason.
1919 */
1920 
1921 /* Local Variables: */
1922    AstPointSet *pset1;       /* PointSet holding input positions for PermMap */
1923    AstPointSet *pset2;       /* PointSet holding output positions for PermMap */
1924    double **ptr1;            /* Pointer to pset1 data */
1925    double **ptr2;            /* Pointer to pset2 data */
1926    double *cnst;             /* Pointer to constants array */
1927    double cn;                /* Potential new constant value */
1928    double ip;                /* Potential output axis index */
1929    double op;                /* Potential input axis index */
1930    int *inprm;               /* Pointer to input axis permutation array */
1931    int *outprm;              /* Pointer to output axis permutation array */
1932    int i;                    /* Axis count */
1933    int nc;                   /* Number of constants stored so far */
1934    int nin;                  /* No. of input coordinates for the PermMap */
1935    int nout;                 /* No. of output coordinates for the PermMap */
1936 
1937 /* Initialise. */
1938    if( outperm ) *outperm = NULL;
1939    if( inperm ) *inperm = NULL;
1940    if( consts ) *consts = NULL;
1941 
1942 /* Check the global error status and the supplied pointers. */
1943    if ( !astOK || !outperm || !inperm || !consts ) return;
1944 
1945 /* Initialise variables to avoid "used of uninitialised variable"
1946    messages from dumb compilers. */
1947    nc = 0;
1948 
1949 /* Get the number of input and output axes for the supplied PermMap. */
1950    nin = astGetNin( map );
1951    nout = astGetNout( map );
1952 
1953 /* Allocate the memory for the returned arrays. */
1954    outprm = (int *) astMalloc( sizeof( int )* (size_t) nout );
1955    inprm = (int *) astMalloc( sizeof( int )* (size_t) nin );
1956    cnst = (double *) astMalloc( sizeof( double )* (size_t) ( nout + nin ) );
1957 
1958 /* Returned the pointers to these arrays.*/
1959    *outperm = outprm;
1960    *inperm = inprm;
1961    *consts = cnst;
1962 
1963 /* Create two PointSets, each holding two points, which can be used for
1964    input and output positions with the PermMap. */
1965    pset1 = astPointSet( 2, nin, "", status );
1966    pset2 = astPointSet( 2, nout, "", status );
1967 
1968 /* Set up the two input positions to be [0,1,2...] and [-1,-1,-1,...]. The
1969    first position is used to enumerate the axes, and the second is used to
1970    check for constant axis values. */
1971    ptr1 = astGetPoints( pset1 );
1972    if( astOK ){
1973       for( i = 0; i < nin; i++ ){
1974          ptr1[ i ][ 0 ] = ( double ) i;
1975          ptr1[ i ][ 1 ] = -1.0;
1976       }
1977    }
1978 
1979 /* Use the PermMap to transform these positions in the forward direction. */
1980    (void) astTransform( map, pset1, 1, pset2 );
1981 
1982 /* Look at the mapped positions to determine the output axis permutation
1983    array. */
1984    ptr2 = astGetPoints( pset2 );
1985    if( astOK ){
1986 
1987 /* No constant axis valeus found yet. */
1988       nc = 0;
1989 
1990 /* Do each output axis. */
1991       for( i = 0; i < nout; i++ ){
1992 
1993 /* If the output axis value is copied from an input axis value, the index
1994    of the appropriate input axis will be in the mapped first position. */
1995          op = ptr2[ i ][ 0 ];
1996 
1997 /* If the output axis value is assigned a constant value, the result of
1998    mapping the two different input axis values will be the same. */
1999          cn = ptr2[ i ][ 1 ];
2000          if( op == cn ) {
2001 
2002 /* We have found another constant. Store it in the constants array, and
2003    store the index of the constant in the output axis permutation array. */
2004             cnst[ nc ] = cn;
2005             outprm[ i ] = -( nc + 1 );
2006             nc++;
2007 
2008 /* If the output axis values are different, then the output axis value
2009    must be copied from the input axis value. */
2010          } else {
2011             outprm[ i ] = (int) ( op + 0.5 );
2012          }
2013       }
2014    }
2015 
2016 /* Now do the same thing to determine the input permutation array. */
2017    if( astOK ){
2018       for( i = 0; i < nout; i++ ){
2019          ptr2[ i ][ 0 ] = ( double ) i;
2020          ptr2[ i ][ 1 ] = -1.0;
2021       }
2022    }
2023 
2024    (void) astTransform( map, pset2, 0, pset1 );
2025 
2026    if( astOK ){
2027 
2028       for( i = 0; i < nin; i++ ){
2029 
2030          ip = ptr1[ i ][ 0 ];
2031          cn = ptr1[ i ][ 1 ];
2032          if( ip == cn ) {
2033 
2034             cnst[ nc ] = cn;
2035             inprm[ i ] = -( nc + 1 );
2036             nc++;
2037 
2038          } else {
2039             inprm[ i ] = (int) ( ip + 0.5 );
2040          }
2041       }
2042    }
2043 
2044 /* Annul the PointSets. */
2045    pset1 = astAnnul( pset1 );
2046    pset2 = astAnnul( pset2 );
2047 
2048 /* If an error has occurred, attempt to free the returned arrays. */
2049    if( !astOK ) {
2050       *outperm = (int *) astFree( (void *) *outperm );
2051       *inperm = (int *) astFree( (void *) *inperm );
2052       *consts = (double *) astFree( (void *) *consts );
2053    }
2054 
2055 /* Return. */
2056    return;
2057 }
2058 
Rate(AstMapping * this,double * at,int ax1,int ax2,int * status)2059 static double Rate( AstMapping *this, double *at, int ax1, int ax2, int *status ){
2060 /*
2061 *  Name:
2062 *     Rate
2063 
2064 *  Purpose:
2065 *     Calculate the rate of change of a Mapping output.
2066 
2067 *  Type:
2068 *     Private function.
2069 
2070 *  Synopsis:
2071 *     #include "winmap.h"
2072 *     result = Rate( AstMapping *this, double *at, int ax1, int ax2, int *status )
2073 
2074 *  Class Membership:
2075 *     WinMap member function (overrides the astRate method inherited
2076 *     from the Mapping class ).
2077 
2078 *  Description:
2079 *     This function returns the rate of change of a specified output of
2080 *     the supplied Mapping with respect to a specified input, at a
2081 *     specified input position.
2082 
2083 *  Parameters:
2084 *     this
2085 *        Pointer to the Mapping to be applied.
2086 *     at
2087 *        The address of an array holding the axis values at the position
2088 *        at which the rate of change is to be evaluated. The number of
2089 *        elements in this array should equal the number of inputs to the
2090 *        Mapping.
2091 *     ax1
2092 *        The index of the Mapping output for which the rate of change is to
2093 *        be found (output numbering starts at 0 for the first output).
2094 *     ax2
2095 *        The index of the Mapping input which is to be varied in order to
2096 *        find the rate of change (input numbering starts at 0 for the first
2097 *        input).
2098 *     status
2099 *        Pointer to the inherited status variable.
2100 
2101 *  Returned Value:
2102 *     The rate of change of Mapping output "ax1" with respect to input
2103 *     "ax2", evaluated at "at", or AST__BAD if the value cannot be
2104 *     calculated.
2105 
2106 */
2107 
2108 /* Local Variables: */
2109    AstWinMap *map;
2110    double result;
2111 
2112 /* Check inherited status */
2113    if( !astOK ) return AST__BAD;
2114 
2115 /* Get a pointer to the WinMap structure. */
2116    map = (AstWinMap *) this;
2117 
2118 /* If the input and output axes are not equal the result is zero. */
2119    if( ax1 != ax2 ) {
2120       result = 0.0;
2121 
2122 /* Otherwise, return the scale factor for the axis, taking the reciprocal
2123    if the WinMap has been inverted. */
2124    } else {
2125       result = ( map->b )[ ax1 ];
2126       if( astGetInvert( map ) ) {
2127          if( result != 0.0 && result != AST__BAD ) {
2128             result = 1.0/result;
2129          } else {
2130             result = AST__BAD;
2131          }
2132       }
2133    }
2134 
2135 /* Return the result. */
2136    return result;
2137 }
2138 
SetAttrib(AstObject * this_object,const char * setting,int * status)2139 static void SetAttrib( AstObject *this_object, const char *setting, int *status ) {
2140 /*
2141 *  Name:
2142 *     SetAttrib
2143 
2144 *  Purpose:
2145 *     Set an attribute value for a WinMap.
2146 
2147 *  Type:
2148 *     Private function.
2149 
2150 *  Synopsis:
2151 *     #include "winmap.h"
2152 *     void SetAttrib( AstObject *this, const char *setting )
2153 
2154 *  Class Membership:
2155 *     WinMap member function (over-rides the astSetAttrib protected
2156 *     method inherited from the Mapping class).
2157 
2158 *  Description:
2159 *     This function assigns an attribute value for a WinMap, the
2160 *     attribute and its value being specified by means of a string of
2161 *     the form:
2162 *
2163 *        "attribute= value "
2164 *
2165 *     Here, "attribute" specifies the attribute name and should be in
2166 *     lower case with no white space present. The value to the right
2167 *     of the "=" should be a suitable textual representation of the
2168 *     value to be assigned and this will be interpreted according to
2169 *     the attribute's data type.  White space surrounding the value is
2170 *     only significant for string attributes.
2171 
2172 *  Parameters:
2173 *     this
2174 *        Pointer to the WinMap.
2175 *     setting
2176 *        Pointer to a null-terminated string specifying the new attribute
2177 *        value.
2178 */
2179 
2180 /* Local Variables: */
2181    AstWinMap *this;             /* Pointer to the WinMap structure */
2182    int len;                      /* Length of setting string */
2183 
2184 /* Check the global error status. */
2185    if ( !astOK ) return;
2186 
2187 /* Obtain a pointer to the WinMap structure. */
2188    this = (AstWinMap *) this_object;
2189 
2190 /* Obtain the length of the setting string. */
2191    len = (int) strlen( setting );
2192 
2193 /* The WinMap class currently has no attributes, so pass it on to the parent
2194    method for further interpretation. */
2195    (*parent_setattrib)( this_object, setting, status );
2196 
2197 }
2198 
TestAttrib(AstObject * this_object,const char * attrib,int * status)2199 static int TestAttrib( AstObject *this_object, const char *attrib, int *status ) {
2200 /*
2201 *  Name:
2202 *     TestAttrib
2203 
2204 *  Purpose:
2205 *     Test if a specified attribute value is set for a WinMap.
2206 
2207 *  Type:
2208 *     Private function.
2209 
2210 *  Synopsis:
2211 *     #include "winmap.h"
2212 *     int TestAttrib( AstObject *this, const char *attrib, int *status )
2213 
2214 *  Class Membership:
2215 *     WinMap member function (over-rides the astTestAttrib protected
2216 *     method inherited from the Mapping class).
2217 
2218 *  Description:
2219 *     This function returns a boolean result (0 or 1) to indicate whether
2220 *     a value has been set for one of a WinMap's attributes.
2221 
2222 *  Parameters:
2223 *     this
2224 *        Pointer to the WinMap.
2225 *     attrib
2226 *        Pointer to a null-terminated string specifying the attribute
2227 *        name.  This should be in lower case with no surrounding white
2228 *        space.
2229 *     status
2230 *        Pointer to the inherited status variable.
2231 
2232 *  Returned Value:
2233 *     One if a value has been set, otherwise zero.
2234 
2235 *  Notes:
2236 *     - A value of zero will be returned if this function is invoked
2237 *     with the global status set, or if it should fail for any reason.
2238 */
2239 
2240 /* Local Variables: */
2241    AstWinMap *this;             /* Pointer to the WinMap structure */
2242    int result;                   /* Result value to return */
2243 
2244 /* Initialise. */
2245    result = 0;
2246 
2247 /* Check the global error status. */
2248    if ( !astOK ) return result;
2249 
2250 /* Obtain a pointer to the WinMap structure. */
2251    this = (AstWinMap *) this_object;
2252 
2253 /* The WinMap class currently has no attributes, so pass it on to the parent
2254    method for further interpretation. */
2255    result = (*parent_testattrib)( this_object, attrib, status );
2256 
2257 /* Return the result, */
2258    return result;
2259 }
2260 
Transform(AstMapping * this,AstPointSet * in,int forward,AstPointSet * out,int * status)2261 static AstPointSet *Transform( AstMapping *this, AstPointSet *in,
2262                                int forward, AstPointSet *out, int *status ) {
2263 /*
2264 *  Name:
2265 *     Transform
2266 
2267 *  Purpose:
2268 *     Apply a WinMap to transform a set of points.
2269 
2270 *  Type:
2271 *     Private function.
2272 
2273 *  Synopsis:
2274 *     #include "winmap.h"
2275 *     AstPointSet *Transform( AstMapping *this, AstPointSet *in,
2276 *                             int forward, AstPointSet *out, int *status )
2277 
2278 *  Class Membership:
2279 *     WinMap member function (over-rides the astTransform protected
2280 *     method inherited from the Mapping class).
2281 
2282 *  Description:
2283 *     This function takes a WinMap and a set of points encapsulated in a
2284 *     PointSet and transforms the points so as to map them into the
2285 *     required window.
2286 
2287 *  Parameters:
2288 *     this
2289 *        Pointer to the WinMap.
2290 *     in
2291 *        Pointer to the PointSet holding the input coordinate data.
2292 *     forward
2293 *        A non-zero value indicates that the forward coordinate transformation
2294 *        should be applied, while a zero value requests the inverse
2295 *        transformation.
2296 *     out
2297 *        Pointer to a PointSet which will hold the transformed (output)
2298 *        coordinate values. A NULL value may also be given, in which case a
2299 *        new PointSet will be created by this function.
2300 *     status
2301 *        Pointer to the inherited status variable.
2302 
2303 *  Returned Value:
2304 *     Pointer to the output (possibly new) PointSet.
2305 
2306 *  Notes:
2307 *     -  A null pointer will be returned if this function is invoked with the
2308 *     global error status set, or if it should fail for any reason.
2309 *     -  The number of coordinate values per point in the input PointSet must
2310 *     match the number of coordinates for the WinMap being applied.
2311 *     -  If an output PointSet is supplied, it must have space for sufficient
2312 *     number of points and coordinate values per point to accommodate the
2313 *     result. Any excess space will be ignored.
2314 */
2315 
2316 /* Local Variables: */
2317    AstPointSet *result;          /* Pointer to output PointSet */
2318    AstWinMap *map;               /* Pointer to WinMap to be applied */
2319    const char *class;            /* Object class */
2320    double **ptr_in;              /* Pointer to input coordinate data */
2321    double **ptr_out;             /* Pointer to output coordinate data */
2322    double *axin;                 /* Pointer to next input axis value */
2323    double *axout;                /* Pointer to next output axis value */
2324    double *a;                    /* Pointer to next constant term */
2325    double *b;                    /* Pointer to next multiplicative term */
2326    double aa;                    /* Constant term */
2327    double bb;                    /* Multiplicative term */
2328    int coord;                    /* Loop counter for coordinates */
2329    int def;                      /* Is mapping defined? */
2330    int ncoord;                   /* Number of coordinates per point */
2331    int npoint;                   /* Number of points */
2332    int point;                    /* Loop counter for points */
2333 
2334 /* Check the global error status. */
2335    if ( !astOK ) return NULL;
2336 
2337 /* Initialise variables to avoid "used of uninitialised variable"
2338    messages from dumb compilers. */
2339    aa = 0.0;
2340    bb = 0.0;
2341 
2342 /* Obtain a pointer to the WinMap. */
2343    map = (AstWinMap *) this;
2344 
2345 /* Apply the parent mapping using the stored pointer to the Transform member
2346    function inherited from the parent Mapping class. This function validates
2347    all arguments and generates an output PointSet if necessary, but does not
2348    actually transform any coordinate values. */
2349    result = (*parent_transform)( this, in, forward, out, status );
2350 
2351 /* We will now extend the parent astTransform method by performing the
2352    calculations needed to generate the output coordinate values. */
2353 
2354 /* Determine the numbers of points and coordinates per point from the input
2355    PointSet and obtain pointers for accessing the input and output coordinate
2356    values. */
2357    ncoord = astGetNcoord( in );
2358    npoint = astGetNpoint( in );
2359    ptr_in = astGetPoints( in );
2360    ptr_out = astGetPoints( result );
2361 
2362 /* Determine whether to apply the forward or inverse mapping, according to the
2363    direction specified and whether the mapping has been inverted. */
2364    if ( astGetInvert( map ) ) forward = !forward;
2365 
2366 /* Report an error if the WinMap does not contain any scales or shifts. */
2367    if( !(map->a && map->b) && astOK ){
2368       class = astGetClass( this );
2369       astError( AST__BADWM, "astTransform(%s): The supplied %s does not "
2370                 "contain any window information.", status, class, class );
2371    }
2372 
2373 /* Perform coordinate arithmetic. */
2374 /* ------------------------------ */
2375    if( astOK ){
2376 
2377 /* Store pointers to the shift and scale for the next axis. */
2378       a = map->a;
2379       b = map->b;
2380 
2381 /* Apply the mapping to each axis. */
2382       for( coord = 0; coord < ncoord; coord++ ){
2383 
2384 /* If either the scale or shift is bad indicate that the mapping is
2385    not defined on this axis. */
2386          if( *a == AST__BAD || *b == AST__BAD ){
2387             def = 0;
2388 
2389 /* Otherwise, get the scale and offset factors for this axis, taking account of
2390    whether the mapping is inverted or not. If the mapping is undefined, set
2391    the "def" flag to indicate this. */
2392          } else {
2393             aa = *a;
2394             bb = *b;
2395 
2396             if( forward ){
2397                def = 1;
2398 
2399             } else if( bb != 0.0 ){
2400                bb = 1.0/bb;
2401                aa = -aa*bb;
2402                def = 1;
2403 
2404             } else {
2405                def = 0;
2406             }
2407 
2408          }
2409 
2410 /* Store pointers to the first inpout and output values on this axis. */
2411          axin = ptr_in[ coord ];
2412          axout = ptr_out[ coord ];
2413 
2414 /* If the mapping is defined, apply it to the supplied points. */
2415          if( def ){
2416 
2417             for( point = 0; point < npoint; point++ ){
2418                if( *axin != AST__BAD ){
2419                   *(axout++) = aa + bb*(*axin);
2420                } else {
2421                   *(axout++) = AST__BAD;
2422                }
2423                axin++;
2424             }
2425 
2426 /* If the mapping is not defined, store bad values on this axis in the
2427    returned points. */
2428          } else {
2429             for( point = 0; point < npoint; point++ ) *(axout++) = AST__BAD;
2430          }
2431 
2432 /* Point to the scale and shift for the next axis. */
2433          a++;
2434          b++;
2435       }
2436 
2437    }
2438 
2439 /* Return a pointer to the output PointSet. */
2440    return result;
2441 }
2442 
WinMat(AstMapping ** maps,int * inverts,int iwm,int * status)2443 static void WinMat( AstMapping **maps, int *inverts, int iwm, int *status ){
2444 /*
2445 *  Name:
2446 *     WinMat
2447 
2448 *  Purpose:
2449 *     Swap a WinMap and a MatrixMap.
2450 
2451 *  Type:
2452 *     Private function.
2453 
2454 *  Synopsis:
2455 *     #include "winmap.h"
2456 *     void WinMat( AstMapping **maps, int *inverts, int iwm, int *status )
2457 
2458 *  Class Membership:
2459 *     WinMap member function
2460 
2461 *  Description:
2462 *     A list of two Mappings is supplied containing a WinMap and a
2463 *     MatrixMap. These Mappings are annulled, and replaced with
2464 *     another pair of Mappings consisting of a WinMap and a MatrixMap
2465 *     in the opposite order. These Mappings are chosen so that their
2466 *     combined effect is the same as the original pair of Mappings.
2467 *     The scale factors in the returned WinMap are always unity (i.e.
2468 *     the differences in scaling get absorbed into the returned
2469 *     MatrixMap).
2470 
2471 *  Parameters:
2472 *     maps
2473 *        A pointer to an array of two Mapping pointers.
2474 *     inverts
2475 *        A pointer to an array of two invert flags.
2476 *     iwm
2477 *        The index within "maps" of the WinMap.
2478 *     status
2479 *        Pointer to the inherited status variable.
2480 
2481 */
2482 
2483 /* Local Variables: */
2484    AstMatrixMap *m1;             /* Pointer to Diagonal scale factor MatrixMap */
2485    AstMatrixMap *m2;             /* Pointer to returned MatrixMap */
2486    AstMatrixMap *sm2;            /* Pointer to simplified returned MatrixMap */
2487    AstMatrixMap *mm;             /* Pointer to the supplied MatrixMap */
2488    AstPointSet *pset1;           /* Shift terms from supplied WinMap */
2489    AstPointSet *pset2;           /* Shift terms for returned WinMap */
2490    AstWinMap *w1;                /* Pointer to the returned WinMap */
2491    AstWinMap *sw1;               /* Pointer to the simplified returned WinMap */
2492    AstWinMap *wm;                /* Pointer to the supplied WinMap */
2493    double **ptr1;                /* Pointer to pset1 data */
2494    double **ptr2;                /* Pointer to pset2 data */
2495    double *a;                    /* Array of shift terms from supplied WinMap */
2496    double *aa;                   /* Pointer to next shift term */
2497    double *b;                    /* Array of scale terms from supplied WinMap */
2498    double *bb;                   /* Pointer to next scale term */
2499    int i;                        /* Axis count */
2500    int nin;                      /* No. of axes in supplied WinMap */
2501    int nout;                     /* No. of axes in returned WinMap */
2502    int old_minv;                 /* Invert value for the supplied MatrixMap */
2503    int old_winv;                 /* Invert value for the supplied WinMap */
2504 
2505 /* Check the global error status. */
2506    if ( !astOK ) return;
2507 
2508 /* Store pointers to the supplied WinMap and the MatrixMap. */
2509    wm = (AstWinMap *) maps[ iwm ];
2510    mm = (AstMatrixMap *) maps[ 1 - iwm ];
2511 
2512 /* Temporarily set the Invert attribute of the supplied Mappings to the
2513    supplied values. */
2514    old_winv = astGetInvert( wm );
2515    astSetInvert( wm, inverts[ iwm ] );
2516 
2517    old_minv = astGetInvert( mm );
2518    astSetInvert( mm, inverts[ 1 - iwm ] );
2519 
2520 /* Get copies of the shift and scale terms used by the WinMap. This
2521    also returns the number of axes in the WinMap. */
2522    nin = astWinTerms( wm, &a, &b );
2523 
2524 /* Create a diagonal MatrixMap holding the scale factors from the
2525    supplied WinMap. */
2526    m1 = astMatrixMap( nin, nin, 1, b, "", status );
2527 
2528 /* Create a PointSet holding a single position given by the shift terms
2529    in the supplied WinMap. */
2530    pset1 = astPointSet( 1, nin, "", status );
2531    ptr1 = astGetPoints( pset1 );
2532    if( astOK ){
2533       aa = a;
2534       for( i = 0; i < nin; i++ ) ptr1[ i ][ 0 ] = *(aa++);
2535    }
2536 
2537 /* First deal with cases when the WinMap is applied first, followed by
2538    the MatrixMap. */
2539    if( iwm == 0 ){
2540 
2541 /* Multiply the diagonal matrix holding the WinMap scale factors by the
2542    supplied matrix. The resulting MatrixMap is the one to return in the
2543    map list. */
2544       m2 = astMtrMult( m1, mm );
2545 
2546 /* Transform the position given by the shift terms from the supplied
2547    WinMap using the supplied MatrixMap to get the shift terms for
2548    the returned WinMap. */
2549       pset2 = astTransform( mm, pset1, 1, NULL );
2550 
2551 /* Now deal with cases when the MatrixMap is applied first, followed by
2552    the WinMap. */
2553    } else {
2554 
2555 /* Multiply the supplied MatrixMap by the diagonal matrix holding scale
2556    factors from the supplied WinMap. The resulting MatrixMap is the one to
2557    return in the map list. */
2558       m2 = astMtrMult( mm, m1 );
2559 
2560 /* Transform the position given by the shift terms from the supplied
2561    WinMap using the inverse of the returned MatrixMap to get the shift
2562    terms for the returned WinMap. */
2563       pset2 = astTransform( m2, pset1, 0, NULL );
2564 
2565    }
2566 
2567 /* Re-instate the original value of the Invert attributes of the supplied
2568    Mappings. */
2569    astSetInvert( wm, old_winv );
2570    astSetInvert( mm, old_minv );
2571 
2572 /* Get pointers to the shift terms for the returned WinMap. */
2573    ptr2 = astGetPoints( pset2 );
2574 
2575 /* Create the returned WinMap, initially with undefined corners. The number of
2576    axes in the WinMap must equal the number of shift terms. */
2577    nout = astGetNcoord( pset2 );
2578    w1 = astWinMap( nout, NULL, NULL, NULL, NULL, "", status );
2579 
2580 /* If succesful, store the scale and shift terms in the WinMap. The scale
2581    terms are always unity. */
2582    if( astOK ){
2583       bb = w1->b;
2584       aa = w1->a;
2585       for( i = 0; i < nout; i++ ) {
2586          *(bb++) = 1.0;
2587          *(aa++) = ptr2[ i ][ 0 ];
2588       }
2589 
2590 /* Replace the supplied Mappings and invert flags with the ones found
2591    above. Remember that the order of the Mappings is now swapped */
2592       (void) astAnnul( maps[ 0 ] );
2593       (void) astAnnul( maps[ 1 ] );
2594 
2595       sw1 = astSimplify( w1 );
2596       w1 = astAnnul( w1 );
2597 
2598       maps[ 1 - iwm ] = (AstMapping *) sw1;
2599       inverts[ 1 - iwm  ] = astGetInvert( sw1 );
2600 
2601       sm2 = astSimplify( m2 );
2602       m2 = astAnnul( m2 );
2603 
2604       maps[ iwm ] = (AstMapping *) sm2;
2605       inverts[ iwm  ] = astGetInvert( sm2 );
2606 
2607    }
2608 
2609 /* Annul the MatrixMap and PointSet holding the scale and shift terms from the
2610    supplied WinMap. */
2611    m1 = astAnnul( m1 );
2612    pset1 = astAnnul( pset1 );
2613    pset2 = astAnnul( pset2 );
2614 
2615 /* Free the copies of the scale and shift terms from the supplied WinMap. */
2616    b = (double *) astFree( (void *) b );
2617    a = (double *) astFree( (void *) a );
2618 
2619 /* Return. */
2620    return;
2621 }
2622 
WinWcs(AstMapping ** maps,int * inverts,int iwm,int * status)2623 static void WinWcs( AstMapping **maps, int *inverts, int iwm, int *status ){
2624 /*
2625 *  Name:
2626 *     WinWcs
2627 
2628 *  Purpose:
2629 *     Swap a WinMap and a WcsMap.
2630 
2631 *  Type:
2632 *     Private function.
2633 
2634 *  Synopsis:
2635 *     #include "winmap.h"
2636 *     void WinWcs( AstMapping **maps, int *inverts, int iwm, int *status )
2637 
2638 *  Class Membership:
2639 *     WinMap member function
2640 
2641 *  Description:
2642 *     A list of two Mappings is supplied containing a WinMap and a
2643 *     WcsMap. These Mappings are swapped.
2644 
2645 *  Parameters:
2646 *     maps
2647 *        A pointer to an array of two Mapping pointers.
2648 *     inverts
2649 *        A pointer to an array of two invert flags.
2650 *     iwm
2651 *        The index within "maps" of the WinMap.
2652 *     status
2653 *        Pointer to the inherited status variable.
2654 
2655 */
2656 
2657 /* Local Variables: */
2658    AstMapping *m1;               /* Pointer to a Mapping */
2659    int inv;                      /* Invert value */
2660 
2661 /* Check the global error status. */
2662    if ( !astOK ) return;
2663 
2664 /* Simply swap the values (the CanSwap function will have checked that
2665    the WcsMap and WinMap can simply be swapped). */
2666    m1 = maps[ 0 ];
2667    maps[ 0 ] = maps[ 1 ];
2668    maps[ 1 ] = m1;
2669 
2670    inv = inverts[ 0 ];
2671    inverts[ 0 ] = inverts[ 1 ];
2672    inverts[ 1 ] = inv;
2673 
2674 /* Return. */
2675    return;
2676 }
2677 
WinPerm(AstMapping ** maps,int * inverts,int iwm,int * status)2678 static void WinPerm( AstMapping **maps, int *inverts, int iwm, int *status ){
2679 /*
2680 *  Name:
2681 *     WinPerm
2682 
2683 *  Purpose:
2684 *     Swap a WinMap and a PermMap.
2685 
2686 *  Type:
2687 *     Private function.
2688 
2689 *  Synopsis:
2690 *     #include "winmap.h"
2691 *     void WinPerm( AstMapping **maps, int *inverts, int iwm, int *status )
2692 
2693 *  Class Membership:
2694 *     WinMap member function
2695 
2696 *  Description:
2697 *     A list of two Mappings is supplied containing a WinMap and a
2698 *     PermMap. These Mappings are annulled, and replaced with
2699 *     another pair of Mappings consisting of a WinMap and a PermMap
2700 *     in the opposite order. These Mappings are chosen so that their
2701 *     combined effect is the same as the original pair of Mappings.
2702 
2703 *  Parameters:
2704 *     maps
2705 *        A pointer to an array of two Mapping pointers.
2706 *     inverts
2707 *        A pointer to an array of two invert flags.
2708 *     iwm
2709 *        The index within "maps" of the WinMap.
2710 *     status
2711 *        Pointer to the inherited status variable.
2712 
2713 *  Notes:
2714 *     -  All links between input and output axes in the PermMap must
2715 *     be bi-directional, but there can be unconnected axes, and there
2716 *     need not be the same number of input and output axes.
2717 
2718 */
2719 
2720 /* Local Variables: */
2721    AstPermMap *pm;               /* Pointer to the supplied PermMap */
2722    AstPermMap *p1;               /* Pointer to the returned PermMap */
2723    AstPermMap *sp1;              /* Pointer to the simplified returned PermMap */
2724    AstWinMap *w1;                /* Pointer to the returned WinMap */
2725    AstWinMap *sw1;               /* Pointer to the simplified returned PermMap */
2726    AstWinMap *wm;                /* Pointer to the supplied WinMap */
2727    double *a;                    /* Array of shift terms from supplied WinMap */
2728    double *aa;                   /* Pointer to next shift term */
2729    double *b;                    /* Array of scale terms from supplied WinMap */
2730    double *bb;                   /* Pointer to next scale term */
2731    double *consts;               /* Pointer to constants array */
2732    double c;                     /* A constant value */
2733    int *inperm;                  /* Pointer to input axis permutation array */
2734    int *outperm;                 /* Pointer to output axis permutation array */
2735    int i;                        /* Axis count */
2736    int j;                        /* Axis index */
2737    int nin;                      /* No. of axes in supplied WinMap */
2738    int npin;                     /* No. of input axes in supplied PermMap */
2739    int npout;                    /* No. of output axes in supplied PermMap */
2740    int old_pinv;                 /* Invert value for the supplied PermMap */
2741    int old_winv;                 /* Invert value for the supplied WinMap */
2742 
2743 
2744 /* Check the global error status. */
2745    if ( !astOK ) return;
2746 
2747 /* Initialise variables to avoid "used of uninitialised variable"
2748    messages from dumb compilers. */
2749    p1 = NULL;
2750    w1 = NULL;
2751 
2752 /* Store pointers to the supplied WinMap and the PermMap. */
2753    wm = (AstWinMap *) maps[ iwm ];
2754    pm = (AstPermMap *) maps[ 1 - iwm ];
2755 
2756 /* Temporarily set the Invert attribute of the supplied Mappings to the
2757    supplied values. */
2758    old_winv = astGetInvert( wm );
2759    astSetInvert( wm, inverts[ iwm ] );
2760 
2761    old_pinv = astGetInvert( pm );
2762    astSetInvert( pm, inverts[ 1 - iwm ] );
2763 
2764 /* Get copies of the shift and scale terms used by the WinMap. This
2765    also returns the number of axes in the WinMap. */
2766    nin = astWinTerms( wm, &a, &b );
2767 
2768 /* Get the axis permutation and constants arrays representing the
2769    PermMap. Note, no constants are used more than once in the returned
2770    arrays (i.e. duplicate constants are returned in "consts" if more than
2771    one axis uses a given constant). */
2772    PermGet( pm, &outperm, &inperm, &consts, status );
2773 
2774    if( astOK ) {
2775 
2776 /* Get the number of input and output axes in the PermMap. */
2777       npin = astGetNin( pm );
2778       npout = astGetNout( pm );
2779 
2780 /* First consider cases where the WinMap is applied first, followed by the
2781    PermMap. */
2782       if( iwm == 0 ) {
2783 
2784 /* Create the new WinMap, initially with undefined corners. Its number
2785    of axes will equal the number of output axes of the PermMap. */
2786          w1 = astWinMap( npout, NULL, NULL, NULL, NULL, "", status );
2787 
2788 /* Get pointers to the scale and shift terms for the new WinMap. */
2789          bb = w1->b;
2790          aa = w1->a;
2791 
2792 /* Thinking of the forward CmpMap first, consider each of the output axes of
2793    the PermMap. */
2794          for( i = 0; i < npout; i++ ){
2795 
2796 /* If the value for this output axis is derived from an input axis, copy the
2797    scale and shift terms from the corresponding input axis to the new
2798    WinMap. */
2799             j = outperm[ i ];
2800             if( j >= 0 && j < nin ) {
2801                aa[ i ] = a[ j ];
2802                bb[ i ] = b[ j ];
2803 
2804 /* If this output axis is assigned a constant value, use zero and one for
2805    the shift and scale in order to preserve the constant value produced
2806    by the PermMap. */
2807             } else {
2808                aa[ i ] = 0.0;
2809                bb[ i ] = 1.0;
2810             }
2811 
2812          }
2813 
2814 /* Now consider the inverse CmpMap. Any constants produced by the inverse
2815    PermMap would previously have been scaled by the inverse WinMap. Since
2816    there will be no inverse WinMap to perform this scaling in the returned
2817    Mappings, we need to change the constant values to be the values after
2818    the scaling which would have been applied by the WinMap. Consider each
2819    of the input axes of the PermMap.*/
2820          for( i = 0; i < npin; i++ ){
2821 
2822 /* Skip axes which are not assigned a constant value. */
2823             if( inperm[ i ] < 0 ) {
2824 
2825 /* Scale the constant term associated with this input axis using the
2826    inverse WinMap unless it is AST__BAD. */
2827                c = consts[ -inperm[ i ] - 1 ];
2828                if( c != AST__BAD ) {
2829 
2830                   if( a[ i ] != AST__BAD && b[ i ] != AST__BAD &&
2831                       b[ i ] != 0.0 ) {
2832                      consts[ -inperm[ i ] - 1 ] = ( c - a[ i ] )/b[ i ];
2833                   } else {
2834                      consts[ -inperm[ i ] - 1 ] = AST__BAD;
2835                   }
2836 
2837                }
2838 
2839             }
2840 
2841          }
2842 
2843 /* Now consider cases where the PermMap is applied first, followed by the
2844    WinMap. */
2845       } else {
2846 
2847 /* Create the new WinMap, initially with undefined corners. Its number
2848    of axes will equal the number of input axes of the PermMap. */
2849          w1 = astWinMap( npin, NULL, NULL, NULL, NULL, "", status );
2850 
2851 /* Get pointers to the scale and shift terms for the new WinMap. */
2852          bb = w1->b;
2853          aa = w1->a;
2854 
2855 /* Thinking first about the inverse WinMap, consider each of the input axes
2856    of the PermMap. */
2857          for( i = 0; i < npin; i++ ){
2858 
2859 /* If the value for this input axis is derived from an output axis, copy the
2860    scale and shift terms from the corresponding output axis to the new
2861    WinMap. */
2862             j = inperm[ i ];
2863             if( j >= 0 && j < nin ) {
2864                aa[ i ] = a[ j ];
2865                bb[ i ] = b[ j ];
2866 
2867 /* If this input axis is assigned a constant value, use zero and one for
2868    the shift and scale in order to preserve the constant value produced
2869    by the PermMap. */
2870             } else {
2871                aa[ i ] = 0.0;
2872                bb[ i ] = 1.0;
2873             }
2874 
2875          }
2876 
2877 /* Now consider the forward WinMap. Any constants produced by the forward
2878    PermMap would previously have been scaled by the forward WinMap. Since
2879    there will be no forward WinMap to perform this scaling in the returned
2880    Mappings, we need to change the constant values to be the values after
2881    the scaling which would have been applied by the WinMap. Consider each
2882    of the output axes of the PermMap.*/
2883          for( i = 0; i < npout; i++ ){
2884 
2885 /* Skip axes which are not assigned a constant value. */
2886             if( outperm[ i ] < 0 ) {
2887 
2888 /* Scale the constant term associated with this input axis using the
2889    forward WinMap unless it is AST__BAD. */
2890                c = consts[ -outperm[ i ] - 1 ];
2891                if( c != AST__BAD ) {
2892 
2893                   if( a[ i ] != AST__BAD && b[ i ] != AST__BAD ) {
2894                      consts[ -outperm[ i ] - 1 ] = a[ i ] + c*b[ i ];
2895                   } else {
2896                      consts[ -outperm[ i ] - 1 ] = AST__BAD;
2897                   }
2898 
2899                }
2900 
2901             }
2902 
2903          }
2904 
2905       }
2906 
2907 /* Create a new PermMap (since the constants may have changed). */
2908       p1 = astPermMap( npin, inperm, npout, outperm, consts, "", status );
2909 
2910 /* Free the axis permutation and constants arrays. */
2911       outperm = (int *) astFree( (void *) outperm );
2912       inperm = (int *) astFree( (void *) inperm );
2913       consts = (double *) astFree( (void *) consts );
2914    }
2915 
2916 /* Re-instate the original value of the Invert attributes of the supplied
2917    Mappings. */
2918    astSetInvert( wm, old_winv );
2919    astSetInvert( pm, old_pinv );
2920 
2921 /* Replace the supplied Mappings with the ones created above, swapping the
2922    order. */
2923    if( astOK ){
2924       (void) astAnnul( wm );
2925       (void) astAnnul( pm );
2926 
2927       sp1 = astSimplify( p1 );
2928       p1 = astAnnul( p1 );
2929 
2930       sw1 = astSimplify( w1 );
2931       w1 = astAnnul( w1 );
2932 
2933       maps[ iwm ] = (AstMapping *) sp1;
2934       inverts[ iwm ] = 0;
2935 
2936       maps[ 1 - iwm ] = (AstMapping *) sw1;
2937       inverts[ 1 - iwm  ] = astGetInvert( sw1 );
2938    }
2939 
2940 /* Free the copies of the scale and shift terms from the supplied WinMap. */
2941    b = (double *) astFree( (void *) b );
2942    a = (double *) astFree( (void *) a );
2943 
2944 /* Return. */
2945    return;
2946 }
2947 
WinTerms(AstWinMap * this,double ** shift,double ** scale,int * status)2948 static int WinTerms( AstWinMap *this, double **shift, double **scale, int *status ){
2949 /*
2950 *+
2951 *  Name:
2952 *     astWinTerms
2953 
2954 *  Purpose:
2955 *     Obtain the scale and shift terms used by a WinMap.
2956 
2957 *  Type:
2958 *     Protected virtual function.
2959 
2960 *  Synopsis:
2961 *     #include "winmap.h"
2962 *     int astWinTerms( AstWinMap *this, double **shift, double **scale )
2963 
2964 *  Class Membership:
2965 *     WinMap mewthod.
2966 
2967 *  Description:
2968 *     This function returns copies of the scale and shift terms used by a
2969 *     WinMap when transforming points. Each axis of the WinMap has a scale
2970 *     term B, and a shift term A, and the transformation of a point is done
2971 *     by applying these to each input axis value X in turn, to get the
2972 *     output axis value B.X + A. The returned terms take into account the
2973 *     current setting of the Invert attribute of the WinMap.
2974 
2975 *  Parameters:
2976 *     this
2977 *        Pointer to the WinMap.
2978 *     shift
2979 *        The address of a location at which to return a pointer to the
2980 *        start of a dynamically allocated array holding the shift terms
2981 *        for each axis.
2982 *     scale
2983 *        The address of a location at which to return a pointer to the
2984 *        start of a dynamically allocated array holding the scale terms
2985 *        for each axis.
2986 
2987 *  Returned Value:
2988 *     The number of axes in the WinMap. This is the same as the number of
2989 *     elements in the returned arrays.
2990 
2991 *  Notes:
2992 *     -  The returned arrays should be released using astFree when no
2993 *     longer needed.
2994 *     -  NULL pointers can be supplied for "scale" or "shift" if the
2995 *     corresponding arrays are not required.
2996 *     -  A value of zero will be returned, together with NULL pointers
2997 *     for "scale" and "shift" if this function is invoked with the
2998 *     global error status set, or if it should fail for any reason.
2999 *-
3000 */
3001 
3002 /* Local Variables: */
3003    double *a;             /* Pointer to a copy of the shift term array */
3004    double *aa;            /* Pointer to the next shift term */
3005    double *b;             /* Pointer to a copy of the scale term array */
3006    double *bb;            /* Pointer to the next scale term */
3007    int i;                 /* Axis count */
3008    int result;            /* The returned number of axes */
3009    size_t absize;         /* Size of shift and scale arrays */
3010 
3011 /* Initialise. */
3012    result = 0;
3013    if( scale ) *scale = NULL;
3014    if( shift ) *shift = NULL;
3015 
3016 /* Check the global status. */
3017    if ( !astOK ) return result;
3018 
3019 /* Get the number of axes in the WinMap. */
3020    result = astGetNin( this );
3021 
3022 /* Create copies of the scale and shift terms from the WinMap. */
3023    absize = sizeof( double )*(size_t) result;
3024    b = (double *) astStore( NULL, (void *) this->b, absize );
3025    a = (double *) astStore( NULL, (void *) this->a, absize );
3026 
3027 /* Check the pointers can be used. */
3028    if( astOK ){
3029 
3030 /* If the WinMap is inverted, replace the scale and shift terms
3031    by the corresponding values for the inverted mapping. */
3032       if( astGetInvert( this ) ){
3033          bb = b;
3034          aa = a;
3035 
3036          for( i = 0; i < result; i++ ){
3037             if( *aa != AST__BAD && *bb != 0.0 && *bb != AST__BAD ){
3038                *bb = 1.0/(*bb);
3039                *aa *= -(*bb);
3040             } else {
3041                *bb = AST__BAD;
3042                *aa = AST__BAD;
3043             }
3044 
3045             aa++;
3046             bb++;
3047 
3048          }
3049       }
3050 
3051 /* Store the required pointers, and free arrays which are not required. */
3052       if( scale ){
3053          *scale = b;
3054       } else {
3055          b = (double *) astFree( (void *) b );
3056       }
3057 
3058       if( shift ){
3059          *shift = a;
3060       } else {
3061          a = (double *) astFree( (void *) a );
3062       }
3063 
3064    }
3065 
3066 /* If an error has occurred, free the arrays and return zero. */
3067    if( !astOK ){
3068       if( scale ) *scale = (double *) astFree( (void *) *scale );
3069       if( shift ) *shift = (double *) astFree( (void *) *shift );
3070       result = 0;
3071    }
3072 
3073 /* Return the answer. */
3074    return result;
3075 
3076 }
3077 
WinUnit(AstWinMap * wm,AstUnitMap * um,int winv,int win1,int * status)3078 static AstWinMap *WinUnit( AstWinMap *wm, AstUnitMap *um, int winv,
3079                            int win1, int *status ){
3080 /*
3081 *  Name:
3082 *     WinUnit
3083 
3084 *  Purpose:
3085 *     Create a WinMap by merging a WinMap and a UnitMap in parallel.
3086 
3087 *  Type:
3088 *     Private function.
3089 
3090 *  Synopsis:
3091 *     #include "winmap.h"
3092 *     AstWinMap *WinUnit( AstWinMap *wm, AstUnitMap *um, int winv, int win1, int *status )
3093 
3094 *  Class Membership:
3095 *     WinMap member function
3096 
3097 *  Description:
3098 *     This function creates a new WinMap which performs a mapping
3099 *     equivalent to applying the two supplied Mappings in parallel in
3100 *     the directions specified by the "invert" flag (the Invert
3101 *     attribute of the supplied WinMap is ignored).
3102 
3103 *  Parameters:
3104 *     wm
3105 *        A pointer to the WinMap.
3106 *     um
3107 *        A pointer to the UnitMap.
3108 *     winv
3109 *        The invert flag to use with wm. A value of zero causes the forward
3110 *        mapping to be used, and a non-zero value causes the inverse
3111 *        mapping to be used.
3112 *     win1
3113 *        Indicates the order in which the Mappings should be applied.
3114 *
3115 *        If win1 is non-zero:
3116 *           "wm" applies to the lower axis indices and "um" to the upper
3117 *           axis indices.
3118 *
3119 *        If win1 is zero:
3120 *           "um" applies to the lower axis indices and "wm" to the upper
3121 *           axis indices.
3122 *     status
3123 *        Pointer to the inherited status variable.
3124 
3125 *  Returned Value:
3126 *     Pointer to the new WinMap.
3127 
3128 *  Notes:
3129 *     -  The forward direction of the returned WinMap is equivalent to the
3130 *     combined effect of the two supplied Mappings, operating in the
3131 *     directions specified by "winv".
3132 *     -  A null pointer will be returned if this function is invoked with the
3133 *     global error status set, or if it should fail for any reason.
3134 */
3135 
3136 /* Local Variables: */
3137    AstWinMap *result;            /* Pointer to output WinMap */
3138    double *a;                    /* Pointer to shift term array */
3139    double *aa;                   /* Pointer to next shift term */
3140    double *ar;                   /* Pointer to next shift term in result */
3141    double *b;                    /* Pointer to scale term array */
3142    double *bb;                   /* Pointer to next scale term */
3143    double *br;                   /* Pointer to next scale term in result */
3144    int i;                        /* Axis index */
3145    int ninw;                     /* No. of axes in the WinMap */
3146    int ninu;                     /* No. of axes in the UnitMap */
3147    int old_winv;                 /* Original setting of WinMap Invert attribute */
3148 
3149 /* Check the global error status. */
3150    if ( !astOK ) return NULL;
3151 
3152 /* Initialise the returned pointer. */
3153    result = NULL;
3154 
3155 /* Temporarily set the Invert attribute of the WinMap to the supplied
3156    value. */
3157    old_winv = astGetInvert( wm );
3158    astSetInvert( wm, winv );
3159 
3160 /* Create copies of the scale and shift terms from the WinMap, and store the
3161    number of axes in it. */
3162    ninw = astWinTerms( wm, &a, &b );
3163 
3164 /* Get the number of axes in the UnitMap. */
3165    ninu = astGetNin( um );
3166 
3167 /* Create the merged WinMap with unspecified corners. */
3168    result = astWinMap( ninw + ninu, NULL, NULL, NULL, NULL, "", status );
3169 
3170 /* Check the pointers can be used. */
3171    if( astOK ){
3172 
3173 /* If the WinMap applies to the lower axis indices... */
3174       if( win1 ){
3175 
3176 /* Use the scale and shift terms from the WinMap for the lower axes of
3177    the new WinMap. */
3178          aa = a;
3179          bb = b;
3180          ar = result->a;
3181          br = result->b;
3182 
3183          for( i = 0; i < ninw; i++ ){
3184             *(ar++) = *(aa++);
3185             *(br++) = *(bb++);
3186          }
3187 
3188 /* Use the scale factor to 1.0 and the shift term to zero for the upper axes
3189    of the new WinMap. */
3190          for( i = 0; i < ninu; i++ ){
3191             *(ar++) = 0.0;
3192             *(br++) = 1.0;
3193          }
3194 
3195 /* If the WinMap applies to the upper axis indices... */
3196       } else {
3197 
3198 /* Use the scale factor to 1.0 and the shift term to zero for the lower axes
3199    of the new WinMap. */
3200          ar = result->a;
3201          br = result->b;
3202 
3203          for( i = 0; i < ninu; i++ ){
3204             *(ar++) = 0.0;
3205             *(br++) = 1.0;
3206          }
3207 
3208 /* Use the scale and shift terms from the WinMap for the upper axes of
3209    the new WinMap. */
3210          aa = a;
3211          bb = b;
3212 
3213          for( i = 0; i < ninw; i++ ){
3214             *(ar++) = *(aa++);
3215             *(br++) = *(bb++);
3216          }
3217       }
3218    }
3219 
3220 /* Free the copies of the scale and shift terms from the supplied WinMap. */
3221    b = (double *) astFree( (void *) b );
3222    a = (double *) astFree( (void *) a );
3223 
3224 /* Re-instate the original setting of the Invert attribute for the
3225    supplied WinMap. */
3226    astSetInvert( wm, old_winv );
3227 
3228 /* If an error has occurred, annull the returned WinMap. */
3229    if( !astOK ) result = astAnnul( result );
3230 
3231 /* Return a pointer to the output WinMap. */
3232    return result;
3233 }
3234 
WinWin(AstMapping * map1,AstMapping * map2,int inv1,int inv2,int series,int * status)3235 static AstWinMap *WinWin( AstMapping *map1, AstMapping *map2, int inv1,
3236                           int inv2, int series, int *status ){
3237 /*
3238 *  Name:
3239 *     WinWin
3240 
3241 *  Purpose:
3242 *     Create a merged WinMap from two supplied WinMaps.
3243 
3244 *  Type:
3245 *     Private function.
3246 
3247 *  Synopsis:
3248 *     #include "winmap.h"
3249 *     AstWinMap *WinWin( AstMapping *map1, AstMapping *map2, int inv1,
3250 *                        int inv2, int series, int *status )
3251 
3252 *  Class Membership:
3253 *     WinMap member function
3254 
3255 *  Description:
3256 *     This function creates a new WinMap which performs a mapping
3257 *     equivalent to applying the two supplied WinMaps either in series
3258 *     or parallel in the directions specified by the "invert" flags
3259 *     (the Invert attributes of the supplied WinMaps are ignored).
3260 
3261 *  Parameters:
3262 *     map1
3263 *        A pointer to the WinMap to apply first (if in series), or to the
3264 *        lower axis indices (if in parallel)
3265 *     map2
3266 *        A pointer to the WinMap to apply second (if in series), or to the
3267 *        upper axis indices (if in parallel)
3268 *     inv1
3269 *        The invert flag to use with map1. A value of zero causes the forward
3270 *        mapping to be used, and a non-zero value causes the inverse
3271 *        mapping to be used.
3272 *     inv2
3273 *        The invert flag to use with map2.
3274 *     series
3275 *        If non-zero, then the supplied WinMaps are combined in series.
3276 *        Otherwise, they are combined in parallel.
3277 *     status
3278 *        Pointer to the inherited status variable.
3279 
3280 *  Returned Value:
3281 *     Pointer to the new WinMap.
3282 
3283 *  Notes:
3284 *     -  The forward direction of the returned WinMap is equivalent to the
3285 *     combined effect of the two supplied WinMap, operating in the
3286 *     directions specified by "inv1" and "inv2".
3287 *     -  A null pointer will be returned if this function is invoked with the
3288 *     global error status set, or if it should fail for any reason.
3289 */
3290 
3291 /* Local Variables: */
3292    AstWinMap *result;            /* Pointer to output WinMap */
3293    AstWinMap *wm1;               /* Pointer to the first supplied WinMap */
3294    AstWinMap *wm2;               /* Pointer to the second supplied WinMap */
3295    double *a[ 2 ];               /* Pointers to shift term arrays */
3296    double *a0;                   /* Pointer to next shift term from WinMap 1 */
3297    double *a1;                   /* Pointer to next shift term from WinMap 2 */
3298    double *ar;                   /* Pointer to next shift term in result */
3299    double *b[ 2 ];               /* Pointers to scale term arrays */
3300    double *b0;                   /* Pointer to next scale term from WinMap 1 */
3301    double *b1;                   /* Pointer to next scale term from WinMap 2 */
3302    double *br;                   /* Pointer to next scale term in result */
3303    int cancel;                   /* Do the two WinMaps cancel out? */
3304    int i;                        /* Axis index */
3305    int invert[ 2 ];              /* Array of invert flags */
3306    int nin[ 2 ];                 /* No. of axes in the two WinMaps */
3307 
3308 /* Check the global error status. */
3309    if ( !astOK ) return NULL;
3310 
3311 /* Initialise the returned pointer. */
3312    result = NULL;
3313 
3314 /* Store pointers to the WinMaps. */
3315    wm1 = (AstWinMap *) map1;
3316    wm2 = (AstWinMap *) map2;
3317 
3318 /* Temporarily set their Invert attributes to the supplied values. */
3319    invert[ 0 ] = astGetInvert( wm1 );
3320    astSetInvert( wm1, inv1 );
3321 
3322    invert[ 1 ] = astGetInvert( wm2 );
3323    astSetInvert( wm2, inv2 );
3324 
3325 /* Create copies of the scale and shift terms from the two WinMaps,
3326    and store the number of axes in each WinMap. The scale and shift terms
3327    returned take into account the setting of the Invert attribute. */
3328    nin[ 0 ] = astWinTerms( wm1, a, b );
3329    nin[ 1 ] = astWinTerms( wm2, a + 1, b + 1 );
3330 
3331 /* Check the pointers can be used. */
3332    if( astOK ){
3333 
3334 /* Series */
3335 /* ====== */
3336       if( series ){
3337 
3338 /* Check for equal and opposite WinMaps. Do this explicitly using the
3339    supplied Mappings rather than the values returned by astWinTerms to
3340    avoid the affects of rounding error sin the inversions performed by
3341    astWinTerms. */
3342          if( ( inv1 == 0 ) != ( inv2 == 0 ) ) {
3343             cancel = 1;
3344             for( i = 0; i < nin[ 0 ]; i++ ){
3345                if( !EQUAL( (wm1->a)[ i ], (wm2->a)[ i ] ) ||
3346                    !EQUAL( (wm1->b)[ i ], (wm2->b)[ i ] ) ) {
3347                   cancel = 0;
3348                   break;
3349                }
3350             }
3351          } else {
3352             cancel = 0;
3353          }
3354 
3355 /* If they cancel, just put unit values into the WinMap. */
3356          if( cancel ) {
3357             a0 = a[ 0 ];
3358             b0 = b[ 0 ];
3359             for( i = 0; i < nin[ 0 ]; i++ ){
3360                *(a0++) = 0.0;
3361                *(b0++) = 1.0;
3362             }
3363 
3364 /* Otherwise, merge the scale and shift terms for the two WinMaps, overwriting
3365    the terms for the first WinMap. To be merged in series, both WinMaps must
3366    have the same number of axes, so it matters not whether we use nin[ 0 ]
3367    or nin[ 1 ] to specify the number of axes. */
3368          } else {
3369             a0 = a[ 0 ];
3370             b0 = b[ 0 ];
3371             a1 = a[ 1 ];
3372             b1 = b[ 1 ];
3373             for( i = 0; i < nin[ 0 ]; i++ ){
3374 
3375                if( *a0 != AST__BAD && *b0 != AST__BAD &&
3376                    *a1 != AST__BAD && *b1 != AST__BAD ){
3377 
3378                   *a0 *= (*b1);
3379                   *a0 += (*a1);
3380                   *b0 *= (*b1);
3381 
3382                } else {
3383                   *a0 = AST__BAD;
3384                   *b0 = AST__BAD;
3385                   *a1 = AST__BAD;
3386                   *b1 = AST__BAD;
3387                }
3388 
3389 /* Move on to the next axis. */
3390                a0++;
3391                b0++;
3392                a1++;
3393                b1++;
3394             }
3395          }
3396 
3397 /* Create the merged WinMap with unspecified corners. */
3398          result = astWinMap( nin[ 0 ], NULL, NULL, NULL, NULL, "", status );
3399 
3400 /* Store the merged scale and shift terms in the new WinMap. The forward
3401    transformation of this WinMap then corresponds to the combination of the
3402    two supplied WinMaps, taking into account their invert flags. */
3403          a0 = a[ 0 ];
3404          b0 = b[ 0 ];
3405          ar = result->a;
3406          br = result->b;
3407          for( i = 0; i < nin[ 0 ]; i++ ){
3408             *(ar++) = *(a0++);
3409             *(br++) = *(b0++);
3410          }
3411 
3412 /* Parallel */
3413 /* ======== */
3414       } else {
3415 
3416 /* Create the merged WinMap with unspecified corners. */
3417          result = astWinMap( nin[ 0 ] + nin[ 1 ], NULL, NULL, NULL, NULL, "", status );
3418 
3419 /* Copy the scale and shift terms into the new WinMap. */
3420          a0 = a[ 0 ];
3421          b0 = b[ 0 ];
3422          a1 = a[ 1 ];
3423          b1 = b[ 1 ];
3424          ar = result->a;
3425          br = result->b;
3426 
3427          for( i = 0; i < nin[ 0 ]; i++ ){
3428             *(ar++) = *(a0++);
3429             *(br++) = *(b0++);
3430          }
3431 
3432          for( i = 0; i < nin[ 1 ]; i++ ){
3433             *(ar++) = *(a1++);
3434             *(br++) = *(b1++);
3435          }
3436       }
3437    }
3438 
3439 /* Re-instate the original settings of the Invert attributes for the
3440    supplied WinMaps. */
3441    astSetInvert( wm1, invert[ 0 ] );
3442    astSetInvert( wm2, invert[ 1 ] );
3443 
3444 /* Free the memory. */
3445    a[ 0 ] = (double *) astFree( (void *) a[ 0 ] );
3446    b[ 0 ] = (double *) astFree( (void *) b[ 0 ] );
3447    a[ 1 ] = (double *) astFree( (void *) a[ 1 ] );
3448    b[ 1 ] = (double *) astFree( (void *) b[ 1 ] );
3449 
3450 /* If an error has occurred, annull the returned WinMap. */
3451    if( !astOK ) result = astAnnul( result );
3452 
3453 /* Return a pointer to the output WinMap. */
3454    return result;
3455 }
3456 
WinZoom(AstWinMap * wm,AstZoomMap * zm,int winv,int zinv,int win1,int series,int * status)3457 static AstWinMap *WinZoom( AstWinMap *wm, AstZoomMap *zm, int winv,
3458                            int zinv, int win1, int series, int *status ){
3459 /*
3460 *  Name:
3461 *     WinZoom
3462 
3463 *  Purpose:
3464 *     Create a WinMap by merging a WinMap and a ZoomMap.
3465 
3466 *  Type:
3467 *     Private function.
3468 
3469 *  Synopsis:
3470 *     #include "winmap.h"
3471 *     AstWinMap *WinZoom( AstWinMap *wm, AstZoomMap *zm, int winv,
3472 *                         int zinv, int win1, int series, int *status )
3473 
3474 *  Class Membership:
3475 *     WinMap member function
3476 
3477 *  Description:
3478 *     This function creates a new WinMap which performs a mapping
3479 *     equivalent to applying the two supplied Mappings in series or
3480 *     parallel in the directions specified by the "invert" flags (the
3481 *     Invert attributes of the supplied WinMaps are ignored).
3482 
3483 *  Parameters:
3484 *     wm
3485 *        A pointer to the WinMap.
3486 *     zm
3487 *        A pointer to the ZoomMap.
3488 *     winv
3489 *        The invert flag to use with wm. A value of zero causes the forward
3490 *        mapping to be used, and a non-zero value causes the inverse
3491 *        mapping to be used.
3492 *     zinv
3493 *        The invert flag to use with zm.
3494 *     win1
3495 *        Indicates the order in which the Mappings should be applied.
3496 *
3497 *        If win1 is non-zero:
3498 *           If in series:
3499 *              "wm" is applied first followed by "zm".
3500 *           If in parallel:
3501 *              "wm" applies to the lower axis indices and "zm" to the upper
3502 *              axis indices.
3503 *
3504 *        If win1 is zero:
3505 *           If in series:
3506 *              "zm" is applied first followed by "wm".
3507 *           If in parallel:
3508 *              "zm" applies to the lower axis indices and "wm" to the upper
3509 *              axis indices.
3510 *     series
3511 *        Should be supplied non-zero if the Mappings are to be combined in
3512 *        series.
3513 *     status
3514 *        Pointer to the inherited status variable.
3515 
3516 *  Returned Value:
3517 *     Pointer to the new WinMap.
3518 
3519 *  Notes:
3520 *     -  The forward direction of the returned WinMap is equivalent to the
3521 *     combined effect of the two supplied Mappings, operating in the
3522 *     directions specified by "zinv" and "winv".
3523 *     -  A null pointer will be returned if this function is invoked with the
3524 *     global error status set, or if it should fail for any reason.
3525 */
3526 
3527 /* Local Variables: */
3528    AstWinMap *result;            /* Pointer to output WinMap */
3529    double *a;                    /* Pointer to shift term array */
3530    double *aa;                   /* Pointer to next shift term */
3531    double *ar;                   /* Pointer to next shift term in result */
3532    double *b;                    /* Pointer to scale term array */
3533    double *bb;                   /* Pointer to next scale term */
3534    double *br;                   /* Pointer to next scale term in result */
3535    double zfac;                  /* Zoom factor */
3536    int i;                        /* Axis index */
3537    int ninw;                     /* No. of axes in the WinMap */
3538    int ninz;                     /* No. of axes in the ZoomMap */
3539    int old_winv;                 /* Original setting of WinMap Invert attribute */
3540    int old_zinv;                 /* Original setting of ZoomMap Invert attribute */
3541 
3542 /* Check the global error status. */
3543    if ( !astOK ) return NULL;
3544 
3545 /* Initialise the returned pointer. */
3546    result = NULL;
3547 
3548 /* Temporarily set the Invert attributes of both Mappings to the supplied
3549    values. */
3550    old_winv = astGetInvert( wm );
3551    astSetInvert( wm, winv );
3552 
3553    old_zinv = astGetInvert( zm );
3554    astSetInvert( zm, zinv );
3555 
3556 /* Get the zoom factor implemented by the ZoomMap. Invert it if necessary
3557    since astGetZoom does not take account of the Invert setting.  */
3558    zfac = astGetZoom( zm );
3559    if( zinv ) zfac = 1.0 / zfac;
3560 
3561 /* Create copies of the scale and shift terms from the WinMap, and store the
3562    number of axes in it. */
3563    ninw = astWinTerms( wm, &a, &b );
3564 
3565 /* Check the pointers can be used. */
3566    if( astOK ){
3567 
3568 /* First do series mode... */
3569       if( series ) {
3570 
3571 /* Modify the WinMap scale and shift terms by the zoom factor. How this is
3572    done depends on which way round the Mappings are applied. */
3573          bb = b;
3574          aa = a;
3575 
3576          for( i = 0; i < ninw; i++ ){
3577 
3578             if( *aa != AST__BAD && *bb != AST__BAD && zfac != AST__BAD ){
3579                *bb *= zfac;
3580                if( win1 ) *aa *= zfac;
3581             } else {
3582                *bb = AST__BAD;
3583                *aa = AST__BAD;
3584             }
3585 
3586             aa++;
3587             bb++;
3588          }
3589 
3590 /* Create the merged WinMap with unspecified corners. */
3591          result = astWinMap( ninw, NULL, NULL, NULL, NULL, "", status );
3592 
3593 /* Store the merged scale and shift terms in the new WinMap. The forward
3594    transformation of this WinMap then corresponds to the combination of the
3595    two supplied Mappings, taking into account their invert flags. */
3596          aa = a;
3597          bb = b;
3598          ar = result->a;
3599          br = result->b;
3600          for( i = 0; i < ninw; i++ ){
3601             *(ar++) = *(aa++);
3602             *(br++) = *(bb++);
3603          }
3604 
3605 /* Now do parallel mode... */
3606       } else {
3607 
3608 /* Get the number of axes in the ZoomMap. */
3609          ninz = astGetNin( zm );
3610 
3611 /* Create the merged WinMap with unspecified corners. */
3612          result = astWinMap( ninw + ninz, NULL, NULL, NULL, NULL, "", status );
3613 
3614 /* If the WinMap applies to the lower axis indices... */
3615          if( win1 ) {
3616 
3617 /* Use the scale and shift terms from the WinMap for the lower axes of
3618    the new WinMap. */
3619             aa = a;
3620             bb = b;
3621             ar = result->a;
3622             br = result->b;
3623 
3624             for( i = 0; i < ninw; i++ ){
3625                *(ar++) = *(aa++);
3626                *(br++) = *(bb++);
3627             }
3628 
3629 /* Use the scale factor (with zero shift) from the ZoomMap for the upper axes
3630    of the new WinMap. */
3631             for( i = 0; i < ninz; i++ ){
3632                *(ar++) = 0.0;
3633                *(br++) = zfac;
3634             }
3635 
3636 /* If the WinMap applies to the upper axis indices... */
3637          } else {
3638 
3639 /* Use the scale factor (with zero shift) from the ZoomMap for the lower axes
3640    of the new WinMap. */
3641             ar = result->a;
3642             br = result->b;
3643 
3644             for( i = 0; i < ninz; i++ ){
3645                *(ar++) = 0.0;
3646                *(br++) = zfac;
3647             }
3648 
3649 /* Use the scale and shift terms from the WinMap for the upper axes of
3650    the new WinMap. */
3651             aa = a;
3652             bb = b;
3653 
3654             for( i = 0; i < ninw; i++ ){
3655                *(ar++) = *(aa++);
3656                *(br++) = *(bb++);
3657             }
3658          }
3659       }
3660    }
3661 
3662 /* Free the copies of the scale and shift terms from the supplied WinMap. */
3663    b = (double *) astFree( (void *) b );
3664    a = (double *) astFree( (void *) a );
3665 
3666 /* Re-instate the original settings of the Invert attribute for the
3667    supplied Mappings. */
3668    astSetInvert( wm, old_winv );
3669    astSetInvert( zm, old_zinv );
3670 
3671 /* If an error has occurred, annull the returned WinMap. */
3672    if( !astOK ) result = astAnnul( result );
3673 
3674 /* Return a pointer to the output WinMap. */
3675    return result;
3676 }
3677 
3678 /* Functions which access class attributes. */
3679 /* ---------------------------------------- */
3680 /* Implement member functions to access the attributes associated with
3681    this class using the macros defined for this purpose in the
3682    "object.h" file. For a description of each attribute, see the class
3683    interface (in the associated .h file). */
3684 
3685 /* Copy constructor. */
3686 /* ----------------- */
Copy(const AstObject * objin,AstObject * objout,int * status)3687 static void Copy( const AstObject *objin, AstObject *objout, int *status ) {
3688 /*
3689 *  Name:
3690 *     Copy
3691 
3692 *  Purpose:
3693 *     Copy constructor for WinMap objects.
3694 
3695 *  Type:
3696 *     Private function.
3697 
3698 *  Synopsis:
3699 *     void Copy( const AstObject *objin, AstObject *objout, int *status )
3700 
3701 *  Description:
3702 *     This function implements the copy constructor for WinMap objects.
3703 
3704 *  Parameters:
3705 *     objin
3706 *        Pointer to the WinMap to be copied.
3707 *     objout
3708 *        Pointer to the WinMap being constructed.
3709 *     status
3710 *        Pointer to the inherited status variable.
3711 
3712 */
3713 
3714 /* Local Variables: */
3715    AstWinMap *out;              /* Pointer to output WinMap */
3716    AstWinMap *in;               /* Pointer to input WinMap */
3717    int ncoord;                  /* No. of axes for the mapping */
3718 
3719 /* Check the global error status. */
3720    if ( !astOK ) return;
3721 
3722 /* Obtain a pointer to the input and output WinMaps. */
3723    in= (AstWinMap *) objin;
3724    out = (AstWinMap *) objout;
3725 
3726 /* Get the number of coordinates mapped by the WinMap. */
3727    ncoord = astGetNin( in );
3728 
3729 /* Allocate memory holding copies of the scales and shifts window defining the
3730    mapping. */
3731    out->a = (double *) astStore( NULL, (void *) in->a,
3732                                   sizeof(double)*(size_t)ncoord );
3733    out->b = (double *) astStore( NULL, (void *) in->b,
3734                                   sizeof(double)*(size_t)ncoord );
3735 
3736 /* If an error occurred, free any allocated memory. */
3737    if ( !astOK ) {
3738       out->a = (double *) astFree( (void *) out->a );
3739       out->b = (double *) astFree( (void *) out->b );
3740    }
3741 
3742 }
3743 
3744 /* Destructor. */
3745 /* ----------- */
Delete(AstObject * obj,int * status)3746 static void Delete( AstObject *obj, int *status ) {
3747 /*
3748 *  Name:
3749 *     Delete
3750 
3751 *  Purpose:
3752 *     Destructor for WinMap objects.
3753 
3754 *  Type:
3755 *     Private function.
3756 
3757 *  Synopsis:
3758 *     void Delete( AstObject *obj, int *status )
3759 
3760 *  Description:
3761 *     This function implements the destructor for WinMap objects.
3762 
3763 *  Parameters:
3764 *     obj
3765 *        Pointer to the WinMap to be deleted.
3766 *     status
3767 *        Pointer to the inherited status variable.
3768 
3769 *  Notes:
3770 *     - This destructor does nothing and exists only to maintain a
3771 *     one-to-one correspondence between destructors and copy
3772 *     constructors.
3773 */
3774 
3775 /* Local Variables: */
3776    AstWinMap *this;              /* Pointer to WinMap */
3777 
3778 /* Obtain a pointer to the WinMap structure. */
3779    this = (AstWinMap *) obj;
3780 
3781 /* Free the memory holding the scales and shifts. */
3782    this->a = (double *) astFree( (void *) this->a );
3783    this->b = (double *) astFree( (void *) this->b );
3784 
3785 }
3786 
3787 /* Dump function. */
3788 /* -------------- */
Dump(AstObject * this_object,AstChannel * channel,int * status)3789 static void Dump( AstObject *this_object, AstChannel *channel, int *status ) {
3790 /*
3791 *  Name:
3792 *     Dump
3793 
3794 *  Purpose:
3795 *     Dump function for WinMap objects.
3796 
3797 *  Type:
3798 *     Private function.
3799 
3800 *  Synopsis:
3801 *     void Dump( AstObject *this, AstChannel *channel, int *status )
3802 
3803 *  Description:
3804 *     This function implements the Dump function which writes out data
3805 *     for the WinMap class to an output Channel.
3806 
3807 *  Parameters:
3808 *     this
3809 *        Pointer to the WinMap whose data are being written.
3810 *     channel
3811 *        Pointer to the Channel to which the data are being written.
3812 *     status
3813 *        Pointer to the inherited status variable.
3814 */
3815 
3816 /* Local Constants: */
3817 #define COMMENT_LEN 50           /* Maximum length of a comment string */
3818 #define KEY_LEN 50               /* Maximum length of a keyword */
3819 
3820 /* Local Variables: */
3821    AstWinMap *this;              /* Pointer to the WinMap structure */
3822    char buff[ KEY_LEN + 1 ];     /* Buffer for keyword string */
3823    char comment[ COMMENT_LEN + 1 ]; /* Buffer for comment string */
3824    int axis;                     /* Axis index */
3825    int ncoord;                   /* No. of axes for mapping */
3826 
3827 /* Check the global error status. */
3828    if ( !astOK ) return;
3829 
3830 /* Obtain a pointer to the WinMap structure. */
3831    this = (AstWinMap *) this_object;
3832 
3833 /* Get the number of coordinates to be mapped. */
3834    ncoord = astGetNin( this );
3835 
3836 /* Write out values representing the instance variables for the
3837    WinMap class.  Accompany these with appropriate comment strings,
3838    possibly depending on the values being written.*/
3839 
3840 /* The scales and shifts. */
3841    for( axis = 0; axis < ncoord; axis++ ){
3842       (void) sprintf( buff, "Sft%d", axis + 1 );
3843       (void) sprintf( comment, "Shift for axis %d", axis + 1 );
3844       astWriteDouble( channel, buff, (this->a)[ axis ] != 0.0, 0,
3845                       (this->a)[ axis ], comment );
3846       (void) sprintf( buff, "Scl%d", axis + 1 );
3847       (void) sprintf( comment, "Scale factor for axis %d", axis + 1 );
3848       astWriteDouble( channel, buff, (this->b)[ axis ] != 1.0, 0,
3849                       (this->b)[ axis ], comment );
3850    }
3851 
3852 /* Undefine macros local to this function. */
3853 #undef COMMENT_LEN
3854 #undef KEY_LEN
3855 }
3856 
3857 /* Standard class functions. */
3858 /* ========================= */
3859 /* Implement the astIsAWinMap and astCheckWinMap functions using the macros
3860    defined for this purpose in the "object.h" header file. */
astMAKE_ISA(WinMap,Mapping)3861 astMAKE_ISA(WinMap,Mapping)
3862 astMAKE_CHECK(WinMap)
3863 
3864 AstWinMap *astWinMap_( int ncoord, const double c1_in[], const double c2_in[],
3865                        const double c1_out[], const double c2_out[],
3866                        const char *options, int *status, ...) {
3867 /*
3868 *++
3869 *  Name:
3870 c     astWinMap
3871 f     AST_WINMAP
3872 
3873 *  Purpose:
3874 *     Create a WinMap.
3875 
3876 *  Type:
3877 *     Public function.
3878 
3879 *  Synopsis:
3880 c     #include "winmap.h"
3881 c     AstWinMap *astWinMap( int ncoord,
3882 c                           const double ina[],  const double inb[],
3883 c                           const double outa[], const double outb[],
3884 c                           const char *options, ... )
3885 f     RESULT = AST_WINMAP( NCOORD, INA, INB, OUTA, OUTB, OPTIONS, STATUS )
3886 
3887 *  Class Membership:
3888 *     WinMap constructor.
3889 
3890 *  Description:
3891 *     This function creates a new WinMap and optionally initialises its
3892 *     attributes.
3893 *
3894 *     A Winmap is a linear Mapping which transforms a rectangular
3895 *     window in one coordinate system into a similar window in another
3896 *     coordinate system by scaling and shifting each axis (the window
3897 *     edges being parallel to the coordinate axes).
3898 *
3899 *     A WinMap is specified by giving the coordinates of two opposite
3900 *     corners (A and B) of the window in both the input and output
3901 *     coordinate systems.
3902 
3903 *  Parameters:
3904 c     ncoord
3905 f     NCOORD = INTEGER (Given)
3906 *        The number of coordinate values for each point to be
3907 *        transformed (i.e. the number of dimensions of the space in
3908 *        which the points will reside). The same number is applicable
3909 *        to both input and output points.
3910 c     ina
3911 f     INA( NCOORD ) = DOUBLE PRECISION (Given)
3912 c        An array containing the "ncoord"
3913 f        An array containing the
3914 *        coordinates of corner A of the window in the input coordinate
3915 *        system.
3916 c     inb
3917 f     INB( NCOORD ) = DOUBLE PRECISION (Given)
3918 c        An array containing the "ncoord"
3919 f        An array containing the
3920 *        coordinates of corner B of the window in the input coordinate
3921 *        system.
3922 c     outa
3923 f     OUTA( NCOORD ) = DOUBLE PRECISION (Given)
3924 c        An array containing the "ncoord"
3925 f        An array containing the
3926 *        coordinates of corner A of the window in the output coordinate
3927 *        system.
3928 c     outb
3929 f     OUTB( NCOORD ) = DOUBLE PRECISION (Given)
3930 c        An array containing the "ncoord"
3931 f        An array containing the
3932 *        coordinates of corner B of the window in the output coordinate
3933 *        system.
3934 c     options
3935 f     OPTIONS = CHARACTER * ( * ) (Given)
3936 c        Pointer to a null-terminated string containing an optional
3937 c        comma-separated list of attribute assignments to be used for
3938 c        initialising the new WinMap. The syntax used is identical to
3939 c        that for the astSet function and may include "printf" format
3940 c        specifiers identified by "%" symbols in the normal way.
3941 f        A character string containing an optional comma-separated
3942 f        list of attribute assignments to be used for initialising the
3943 f        new WinMap. The syntax used is identical to that for the
3944 f        AST_SET routine.
3945 c     ...
3946 c        If the "options" string contains "%" format specifiers, then
3947 c        an optional list of additional arguments may follow it in
3948 c        order to supply values to be substituted for these
3949 c        specifiers. The rules for supplying these are identical to
3950 c        those for the astSet function (and for the C "printf"
3951 c        function).
3952 f     STATUS = INTEGER (Given and Returned)
3953 f        The global status.
3954 
3955 *  Returned Value:
3956 c     astWinMap()
3957 f     AST_WINMAP = INTEGER
3958 *        A pointer to the new WinMap.
3959 
3960 *  Notes:
3961 *     - A null Object pointer (AST__NULL) will be returned if this
3962 c     function is invoked with the AST error status set, or if it
3963 f     function is invoked with STATUS set to an error value, or if it
3964 *     should fail for any reason.
3965 
3966 *  Status Handling:
3967 *     The protected interface to this function includes an extra
3968 *     parameter at the end of the parameter list descirbed above. This
3969 *     parameter is a pointer to the integer inherited status
3970 *     variable: "int *status".
3971 
3972 *--
3973 */
3974 
3975 /* Local Variables: */
3976    astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
3977    AstWinMap *new;              /* Pointer to new WinMap */
3978    va_list args;                /* Variable argument list */
3979 
3980 /* Get a pointer to the thread specific global data structure. */
3981    astGET_GLOBALS(NULL);
3982 
3983 /* Check the global status. */
3984    if ( !astOK ) return NULL;
3985 
3986 /* Initialise the WinMap, allocating memory and initialising the
3987    virtual function table as well if necessary. */
3988    new = astInitWinMap( NULL, sizeof( AstWinMap ), !class_init, &class_vtab,
3989                         "WinMap", ncoord, c1_in, c2_in, c1_out, c2_out );
3990 
3991 /* If successful, note that the virtual function table has been
3992    initialised. */
3993    if ( astOK ) {
3994       class_init = 1;
3995 
3996 /* Obtain the variable argument list and pass it along with the options string
3997    to the astVSet method to initialise the new WinMap's attributes. */
3998       va_start( args, status );
3999       astVSet( new, options, NULL, args );
4000       va_end( args );
4001 
4002 /* If an error occurred, clean up by deleting the new object. */
4003       if ( !astOK ) new = astDelete( new );
4004    }
4005 
4006 /* Return a pointer to the new WinMap. */
4007    return new;
4008 }
4009 
astWinMapId_(int ncoord,const double c1_in[],const double c2_in[],const double c1_out[],const double c2_out[],const char * options,...)4010 AstWinMap *astWinMapId_( int ncoord, const double c1_in[], const double c2_in[],
4011                          const double c1_out[], const double c2_out[],
4012                          const char *options, ... ) {
4013 /*
4014 *  Name:
4015 *     astWinMapId_
4016 
4017 *  Purpose:
4018 *     Create a WinMap.
4019 
4020 *  Type:
4021 *     Private function.
4022 
4023 *  Synopsis:
4024 *     #include "winmap.h"
4025 *     AstWinMap *astWinMapId_( int ncoord, const double c1_in[],
4026 *                              const double c2_in[], const double c1_out[],
4027 *                              const double c2_out[],
4028 *                              const char *options, ... )
4029 
4030 *  Class Membership:
4031 *     WinMap constructor.
4032 
4033 *  Description:
4034 *     This function implements the external (public) interface to the
4035 *     astWinMap constructor function. It returns an ID value (instead
4036 *     of a true C pointer) to external users, and must be provided
4037 *     because astWinMap_ has a variable argument list which cannot be
4038 *     encapsulated in a macro (where this conversion would otherwise
4039 *     occur).
4040 *
4041 *     The variable argument list also prevents this function from
4042 *     invoking astWinMap_ directly, so it must be a re-implementation
4043 *     of it in all respects, except for the final conversion of the
4044 *     result to an ID value.
4045 
4046 *  Parameters:
4047 *     As for astWinMap_.
4048 
4049 *  Returned Value:
4050 *     The ID value associated with the new WinMap.
4051 */
4052 
4053 /* Local Variables: */
4054    astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
4055    AstWinMap *new;              /* Pointer to new WinMap */
4056    va_list args;                /* Variable argument list */
4057    int *status;                 /* Pointer to inherited status value */
4058 
4059 /* Get a pointer to the inherited status value. */
4060    status = astGetStatusPtr;
4061 
4062 /* Get a pointer to the thread specific global data structure. */
4063    astGET_GLOBALS(NULL);
4064 
4065 /* Check the global status. */
4066    if ( !astOK ) return NULL;
4067 
4068 /* Initialise the WinMap, allocating memory and initialising the
4069    virtual function table as well if necessary. */
4070    new = astInitWinMap( NULL, sizeof( AstWinMap ), !class_init, &class_vtab,
4071                         "WinMap", ncoord, c1_in, c2_in, c1_out, c2_out );
4072 
4073 /* If successful, note that the virtual function table has been
4074    initialised. */
4075    if ( astOK ) {
4076       class_init = 1;
4077 
4078 /* Obtain the variable argument list and pass it along with the options string
4079    to the astVSet method to initialise the new WinMap's attributes. */
4080       va_start( args, options );
4081       astVSet( new, options, NULL, args );
4082       va_end( args );
4083 
4084 /* If an error occurred, clean up by deleting the new object. */
4085       if ( !astOK ) new = astDelete( new );
4086    }
4087 
4088 /* Return an ID value for the new WinMap. */
4089    return astMakeId( new );
4090 }
4091 
astInitWinMap_(void * mem,size_t size,int init,AstWinMapVtab * vtab,const char * name,int ncoord,const double * c1_in,const double * c2_in,const double * c1_out,const double * c2_out,int * status)4092 AstWinMap *astInitWinMap_( void *mem, size_t size, int init,
4093                            AstWinMapVtab *vtab, const char *name,
4094                            int ncoord, const double *c1_in,
4095                            const double *c2_in, const double *c1_out,
4096                            const double *c2_out, int *status ) {
4097 /*
4098 *+
4099 *  Name:
4100 *     astInitWinMap
4101 
4102 *  Purpose:
4103 *     Initialise a WinMap.
4104 
4105 *  Type:
4106 *     Protected function.
4107 
4108 *  Synopsis:
4109 *     #include "winmap.h"
4110 *     AstWinMap *astInitWinMap( void *mem, size_t size, int init,
4111 *                               AstWinMapVtab *vtab, const char *name,
4112 *                               int ncoord, const double *c1_in,
4113 *                               const double *c2_in,
4114 *                               const double *c1_out, const double *c2_out )
4115 
4116 *  Class Membership:
4117 *     WinMap initialiser.
4118 
4119 *  Description:
4120 *     This function is provided for use by class implementations to initialise
4121 *     a new WinMap object. It allocates memory (if necessary) to accommodate
4122 *     the WinMap plus any additional data associated with the derived class.
4123 *     It then initialises a WinMap structure at the start of this memory. If
4124 *     the "init" flag is set, it also initialises the contents of a virtual
4125 *     function table for a WinMap at the start of the memory passed via the
4126 *     "vtab" parameter.
4127 
4128 *  Parameters:
4129 *     mem
4130 *        A pointer to the memory in which the WinMap is to be initialised.
4131 *        This must be of sufficient size to accommodate the WinMap data
4132 *        (sizeof(WinMap)) plus any data used by the derived class. If a value
4133 *        of NULL is given, this function will allocate the memory itself using
4134 *        the "size" parameter to determine its size.
4135 *     size
4136 *        The amount of memory used by the WinMap (plus derived class data).
4137 *        This will be used to allocate memory if a value of NULL is given for
4138 *        the "mem" parameter. This value is also stored in the WinMap
4139 *        structure, so a valid value must be supplied even if not required for
4140 *        allocating memory.
4141 *     init
4142 *        A logical flag indicating if the WinMap's virtual function table is
4143 *        to be initialised. If this value is non-zero, the virtual function
4144 *        table will be initialised by this function.
4145 *     vtab
4146 *        Pointer to the start of the virtual function table to be associated
4147 *        with the new WinMap.
4148 *     name
4149 *        Pointer to a constant null-terminated character string which contains
4150 *        the name of the class to which the new object belongs (it is this
4151 *        pointer value that will subsequently be returned by the astGetClass
4152 *        method).
4153 *     ncoord
4154 *        The number of coordinate values per point.
4155 *     c1_in
4156 *        The input coordinates of corner C1 of the window.
4157 *     c2_in
4158 *        The input coordinates of corner C2 of the window.
4159 *     c1_out
4160 *        The output coordinates of corner C1 of the window.
4161 *     c2_out
4162 *        The output coordinates of corner C2 of the window.
4163 
4164 *  Returned Value:
4165 *     A pointer to the new WinMap.
4166 
4167 *  Notes:
4168 *     -  A null pointer will be returned if this function is invoked with the
4169 *     global error status set, or if it should fail for any reason.
4170 *-
4171 */
4172 
4173 /* Local Variables: */
4174    AstWinMap *new;              /* Pointer to new WinMap */
4175    double denom;                /* Denominotor */
4176    int axis;                    /* Axis index */
4177 
4178 /* Check the global status. */
4179    if ( !astOK ) return NULL;
4180 
4181 /* If necessary, initialise the virtual function table. */
4182    if ( init ) astInitWinMapVtab( vtab, name );
4183 
4184 /* Initialise. */
4185    new = NULL;
4186 
4187 /* Initialise a Mapping structure (the parent class) as the first component
4188    within the WinMap structure, allocating memory if necessary. Specify that
4189    the Mapping should be defined in both the forward and inverse directions. */
4190    new = (AstWinMap *) astInitMapping( mem, size, 0,
4191                                        (AstMappingVtab *) vtab, name,
4192                                        ncoord, ncoord, 1, 1 );
4193 
4194    if ( astOK ) {
4195 
4196 /* Initialise the WinMap data. */
4197 /* ---------------------------- */
4198 /* Allocate memory to hold the shift and scale for each axis. */
4199       new->a = (double *) astMalloc( sizeof(double)*(size_t)ncoord );
4200       new->b = (double *) astMalloc( sizeof(double)*(size_t)ncoord );
4201 
4202 /* Check the pointers can be used */
4203       if( astOK ){
4204 
4205 /* Calculater and store the shift and scale for each axis. */
4206          for( axis = 0; axis < ncoord; axis++ ){
4207 
4208 /* If any of the corners have not been provided, store bad values. */
4209             if( !c1_in || !c1_out || !c2_in || !c2_out ) {
4210                (new->b)[ axis ] = AST__BAD;
4211                (new->a)[ axis ] = AST__BAD;
4212 
4213 /* Otherwise, check the corners are good (not AST__BAD or NaN)... */
4214             } else if( astISGOOD(c2_in[ axis ]) && astISGOOD(c1_in[ axis ]) &&
4215                        astISGOOD(c2_out[ axis ]) && astISGOOD(c1_out[ axis ]) ){
4216 
4217                denom = c2_in[ axis ] - c1_in[ axis ];
4218                if( denom != 0.0 ){
4219                   (new->b)[ axis ] = ( c2_out[ axis ] - c1_out[ axis ] )/denom;
4220                   (new->a)[ axis ] = c1_out[ axis ] - (new->b)[ axis ]*c1_in[ axis ];
4221                } else {
4222                   (new->b)[ axis ] = AST__BAD;
4223                   (new->a)[ axis ] = AST__BAD;
4224                }
4225 
4226             } else {
4227                (new->b)[ axis ] = AST__BAD;
4228                (new->a)[ axis ] = AST__BAD;
4229             }
4230 
4231          }
4232 
4233       }
4234 
4235 /* If an error occurred, clean up by deleting the new WinMap. */
4236       if ( !astOK ) new = astDelete( new );
4237    }
4238 
4239 /* Return a pointer to the new WinMap. */
4240    return new;
4241 }
4242 
astLoadWinMap_(void * mem,size_t size,AstWinMapVtab * vtab,const char * name,AstChannel * channel,int * status)4243 AstWinMap *astLoadWinMap_( void *mem, size_t size,
4244                            AstWinMapVtab *vtab, const char *name,
4245                            AstChannel *channel, int *status ) {
4246 /*
4247 *+
4248 *  Name:
4249 *     astLoadWinMap
4250 
4251 *  Purpose:
4252 *     Load a WinMap.
4253 
4254 *  Type:
4255 *     Protected function.
4256 
4257 *  Synopsis:
4258 *     #include "winmap.h"
4259 *     AstWinMap *astLoadWinMap( void *mem, size_t size,
4260 *                               AstWinMapVtab *vtab, const char *name,
4261 *                               AstChannel *channel )
4262 
4263 *  Class Membership:
4264 *     WinMap loader.
4265 
4266 *  Description:
4267 *     This function is provided to load a new WinMap using data read
4268 *     from a Channel. It first loads the data used by the parent class
4269 *     (which allocates memory if necessary) and then initialises a
4270 *     WinMap structure in this memory, using data read from the input
4271 *     Channel.
4272 *
4273 *     If the "init" flag is set, it also initialises the contents of a
4274 *     virtual function table for a WinMap at the start of the memory
4275 *     passed via the "vtab" parameter.
4276 
4277 
4278 *  Parameters:
4279 *     mem
4280 *        A pointer to the memory into which the WinMap is to be
4281 *        loaded.  This must be of sufficient size to accommodate the
4282 *        WinMap data (sizeof(WinMap)) plus any data used by derived
4283 *        classes. If a value of NULL is given, this function will
4284 *        allocate the memory itself using the "size" parameter to
4285 *        determine its size.
4286 *     size
4287 *        The amount of memory used by the WinMap (plus derived class
4288 *        data).  This will be used to allocate memory if a value of
4289 *        NULL is given for the "mem" parameter. This value is also
4290 *        stored in the WinMap structure, so a valid value must be
4291 *        supplied even if not required for allocating memory.
4292 *
4293 *        If the "vtab" parameter is NULL, the "size" value is ignored
4294 *        and sizeof(AstWinMap) is used instead.
4295 *     vtab
4296 *        Pointer to the start of the virtual function table to be
4297 *        associated with the new WinMap. If this is NULL, a pointer
4298 *        to the (static) virtual function table for the WinMap class
4299 *        is used instead.
4300 *     name
4301 *        Pointer to a constant null-terminated character string which
4302 *        contains the name of the class to which the new object
4303 *        belongs (it is this pointer value that will subsequently be
4304 *        returned by the astGetClass method).
4305 *
4306 *        If the "vtab" parameter is NULL, the "name" value is ignored
4307 *        and a pointer to the string "WinMap" is used instead.
4308 
4309 *  Returned Value:
4310 *     A pointer to the new WinMap.
4311 
4312 *  Notes:
4313 *     - A null pointer will be returned if this function is invoked
4314 *     with the global error status set, or if it should fail for any
4315 *     reason.
4316 *-
4317 */
4318 
4319 /* Local Constants. */
4320    astDECLARE_GLOBALS            /* Pointer to thread-specific global data */
4321 #define KEY_LEN 50               /* Maximum length of a keyword */
4322 
4323 /* Local Variables: */
4324    AstWinMap *new;              /* Pointer to the new WinMap */
4325    char buff[ KEY_LEN + 1 ];    /* Buffer for keyword string */
4326    int axis;                    /* Axis index */
4327    int ncoord;                  /* The number of coordinate axes */
4328 
4329 /* Get a pointer to the thread specific global data structure. */
4330    astGET_GLOBALS(channel);
4331 
4332 /* Initialise. */
4333    new = NULL;
4334 
4335 /* Check the global error status. */
4336    if ( !astOK ) return new;
4337 
4338 /* If a NULL virtual function table has been supplied, then this is
4339    the first loader to be invoked for this WinMap. In this case the
4340    WinMap belongs to this class, so supply appropriate values to be
4341    passed to the parent class loader (and its parent, etc.). */
4342    if ( !vtab ) {
4343       size = sizeof( AstWinMap );
4344       vtab = &class_vtab;
4345       name = "WinMap";
4346 
4347 /* If required, initialise the virtual function table for this class. */
4348       if ( !class_init ) {
4349          astInitWinMapVtab( vtab, name );
4350          class_init = 1;
4351       }
4352    }
4353 
4354 /* Invoke the parent class loader to load data for all the ancestral
4355    classes of the current one, returning a pointer to the resulting
4356    partly-built WinMap. */
4357    new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name,
4358                          channel );
4359 
4360    if ( astOK ) {
4361 
4362 /* Get the number of axis for the mapping. */
4363       ncoord = astGetNin( (AstMapping *) new );
4364 
4365 /* Allocate memory to hold the scales and shifts. */
4366       new->a = (double *) astMalloc( sizeof(double)*(size_t)ncoord );
4367       new->b = (double *) astMalloc( sizeof(double)*(size_t)ncoord );
4368 
4369 /* Read input data. */
4370 /* ================ */
4371 /* Request the input Channel to read all the input data appropriate to
4372    this class into the internal "values list". */
4373       astReadClassData( channel, "WinMap" );
4374 
4375 /* Now read each individual data item from this list and use it to
4376    initialise the appropriate instance variable(s) for this class. */
4377 
4378 /* The scales and shifts. */
4379       for( axis = 0; axis < ncoord; axis++ ){
4380          (void) sprintf( buff, "sft%d", axis + 1 );
4381          (new->a)[ axis ] = astReadDouble( channel, buff, 0.0 );
4382          (void) sprintf( buff, "scl%d", axis + 1 );
4383          (new->b)[ axis ] = astReadDouble( channel, buff, 1.0 );
4384       }
4385    }
4386 
4387 /* If an error occurred, clean up by deleting the new WinMap. */
4388    if ( !astOK ) new = astDelete( new );
4389 
4390 /* Return the new WinMap pointer. */
4391    return new;
4392 
4393 /* Undefine macros local to this function. */
4394 #undef KEY_LEN
4395 }
4396 
4397 /* Virtual function interfaces. */
4398 /* ============================ */
4399 /* These provide the external interface to the virtual functions defined by
4400    this class. Each simply checks the global error status and then locates and
4401    executes the appropriate member function, using the function pointer stored
4402    in the object's virtual function table (this pointer is located using the
4403    astMEMBER macro defined in "object.h").
4404 
4405    Note that the member function may not be the one defined here, as it may
4406    have been over-ridden by a derived class. However, it should still have the
4407    same interface. */
4408 
astWinTerms_(AstWinMap * this,double ** scale,double ** shift,int * status)4409 int astWinTerms_( AstWinMap *this, double **scale, double **shift, int *status ){
4410    if( !astOK ) return 0;
4411    return (**astMEMBER(this,WinMap,WinTerms))( this, scale, shift, status );
4412 }
4413 
4414 
4415 
4416 
4417