diff options
Diffstat (limited to 'ast/matrixmap.c')
| -rw-r--r-- | ast/matrixmap.c | 5767 |
1 files changed, 5767 insertions, 0 deletions
diff --git a/ast/matrixmap.c b/ast/matrixmap.c new file mode 100644 index 0000000..ed9600b --- /dev/null +++ b/ast/matrixmap.c @@ -0,0 +1,5767 @@ +/* +*class++ +* Name: +* MatrixMap + +* Purpose: +* Map coordinates by multiplying by a matrix. + +* Constructor Function: +c astMatrixMap +f AST_MATRIXMAP + +* Description: +* A MatrixMap is form of Mapping which performs a general linear +* transformation. Each set of input coordinates, regarded as a +* column-vector, are pre-multiplied by a matrix (whose elements +* are specified when the MatrixMap is created) to give a new +* column-vector containing the output coordinates. If appropriate, +* the inverse transformation may also be performed. + +* Inheritance: +* The MatrixMap class inherits from the Mapping class. + +* Attributes: +* The MatrixMap class does not define any new attributes beyond +* those which are applicable to all Mappings. + +* Functions: +c The MatrixMap class does not define any new functions beyond those +f The MatrixMap class does not define any new routines beyond those +* which are applicable to all Mappings. + +* Copyright: +* Copyright (C) 1997-2006 Council for the Central Laboratory of the +* Research Councils +* Copyright (C) 2009 Science & Technology Facilities Council. +* All Rights Reserved. + +* Licence: +* This program is free software: you can redistribute it and/or +* modify it under the terms of the GNU Lesser General Public +* License as published by the Free Software Foundation, either +* version 3 of the License, or (at your option) any later +* version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU Lesser General Public License for more details. +* +* You should have received a copy of the GNU Lesser General +* License along with this program. If not, see +* <http://www.gnu.org/licenses/>. + +* Authors: +* DSB: D.S. Berry (Starlink) +* RFWS: R.F. Warren-Smith (Starlink) + +* History: +* 9-FEB-1996 (DSB): +* Original version. +* 13-NOV-1996 (DSB): +* Updated to support attributes, I/O and an external interface. +* 3-JUN-1997 (DSB): +* astMtrMult and astMtrRot made protected instead of public. +* 16-JUN-1997 (RFWS): +* Tidied public prologues. +* 24-JUN-1997 (DSB): +* Zero returned for coordinates which are indeterminate as a +* result of using an inverted, non-square, diagonal matrix. +* 10-OCT-1997 (DSB): +* o The inverse matrix is no longer dumped by the Dump function. +* Instead, it is re-calculated by the Load function. +* o The description of argument "form" in astMatrixMap corrected +* to indicate that a value of 2 produces a unit matrix. +* o String values used to represent choices externally, instead +* of integers. +* 24-NOV-1997 (DSB): +* Use of error code AST__OPT replaced by AST__RDERR. +* 28-JAN-1998 (DSB): +* Bug fix in astMtrMult: the matrix (forward or inverse) used for +* the "a" MatrixMap was determined by the Invert flag of the other +* ("this") MatrixMap. +* 14-APR-1998 (DSB): +* Bug fix in Dump. Previously, matrix elements with value AST__BAD +* were explicitly written out. Now they are not written out, since +* AST__BAD can have different values on different machines. Missing +* elements default to AST__BAD when read back in using astLoadMatrixMap. +* 20-APR-1998 (DSB): +* Bug fix in astLoadMatrixMap: initialise the pointer to the inverse +* matrix array to NULL if no inverse matrix is needed. +* 25-AUG-1998 (DSB): +* - Transform changed so that bad input axis values are not +* propagated to output axes which are independant of the input axis. +* - CompressMatrix changed to allow a tolerance of DBL_EPSILON when +* determining if a matrix is a unit matrix, or a diagonal matrix. +* - MapMerge changed to allow MatrixMaps to swap with PermMaps +* in order to move the MatrixMap closer to a Mapping with which it +* could merge. +* 22-FEB-1999 (DSB): +* Changed logic of MapMerge to avoid infinite looping. +* 5-MAY-1999 (DSB): +* More corrections to MapMerge: Cleared up errors in the use of the +* supplied invert flags, and corrected logic for deciding which +* neighbouring Mapping to swap with. +* 16-JUL-1999 (DSB): +* Fixed memory leaks in MatWin and MapMerge. +* 8-JAN-2003 (DSB): +* Changed private InitVtab method to protected astInitatrixMapVtab +* method. +* 11-SEP-2003 (DSB): +* Increased tolerance on checks for unit matrices within +* CompressMatrix. Now uses sqrt(DBL_EPSILON)*diag (previously was +* DBL_EPSILON*DIAG ). +* 10-NOV-2003 (DSB): +* Modified functions which swap a MatrixMap with another Mapping +* (e.g. MatSwapPerm, etc), to simplify the returned Mappings. +* 13-JAN-2003 (DSB): +* Modified the tolerance used by CompressMatrix when checking for +* zero matrix elements. Old system compared each element to thre +* size of the diagonal, but different scalings on different axes could +* cause this to trat as zero values which should nto be treated as +* zero. +* 23-APR-2004 (DSB): +* Changes to simplification algorithm. +* 8-JUL-2004 (DSB): +* astMtrMult - Report an error if either MatrixMap does not have a +* defined forward transformation. +* 1-SEP-2004 (DSB): +* Ensure do1 and do2 are initialised before use in MapMerge. +* 7-SEP-2005 (DSB): +* Take account of the Invert flag when using the zoom factor from +* a ZoomMap. +* 14-FEB-2006 (DSB): +* Correct row/col confusion in CompressMatrix. +* 15-MAR-2006 (DSB): +* Override astEqual. +* 15-MAR-2009 (DSB): +* MapSplit: Only create the returned Mapping if it would have some +* outputs. Also, do not create the returned Mapping if any output +* depends on a mixture of selected and unselected inputs. +* 16-JUL-2009 (DSB): +* MatPerm: Fix memory leak (mm2 was not being annulled). +* 2-OCT-2012 (DSB): +* - Check for Infs as well as NaNs. +* - In MapSplit do not split the MatrixMap if the resulting +* matrix would contain only bad elements. +* - Report an error if an attempt is made to create a MatrixMap +* containing only bad elements. +* 4-NOV-2013 (DSB): +* Allow a full form MatrixMap to be simplified to a diagonal form +* MatrixMap if all the off-diagonal values are zero. +* 23-APR-2015 (DSB): +* Improve MapMerge. If a MatrixMap can merge with its next-but-one +* neighbour, then swap the MatrixMap with its neighbour, so that +* it is then next its next-but-one neighbour, and then merge the +* two Mappings into a single Mapping. Previously, only the swap +* was performed - not the merger. And the swap was only performed +* if the intervening neighbour could not itself merge. This could +* result in an infinite simplification loop, which was detected by +* CmpMap and and aborted, resulting in no useful simplification. +* 15-JUN-2017 (DSB): +* A diagonal MatrixMap in which the diagonal elements are all zero +* cannot be simplified to a ZoomMap, since ZoomMaps cannot have +* zero zoom factor. +* 16-JUN-2017 (DSB): +* Fix error checking bug in MtrMult - it was checking for the +* inverse transformation of "this" instead of the forward +* transformation of "a". +* 7-NOW-2017 (DSB): +* Allow a diagonal MatrixMap to merge with a WinMap. +* 5-JUN-2018 (DSB): +* Include the inverse matrix in the dump of a MatrixMap. Previously, +* the inverse matrix was calculated afresh using function InvertMatrix +* when a MatrixMap was read from a dump. However this could introduce +* small round-trip errors if the inverse matrix in the original +* MatrixMap was created by astMtrRot etc, rather than the InvertMatrix +* function. +*class-- +*/ + +/* Module Macros. */ +/* ============== */ +/* Set the name of the class we are implementing. This indicates to + the header files that define class interfaces that they should make + "protected" symbols available. */ +#define astCLASS MatrixMap + +/* Define identifiers for the different forms of matrix storage. */ +#define FULL 0 +#define DIAGONAL 1 +#define UNIT 2 + +/* Include files. */ +/* ============== */ +/* Interface definitions. */ +/* ---------------------- */ + +#include "globals.h" /* Thread-safe global data access */ +#include "error.h" /* Error reporting facilities */ +#include "memory.h" /* Memory allocation facilities */ +#include "object.h" /* Base Object class */ +#include "pointset.h" /* Sets of points/coordinates */ +#include "mapping.h" /* Coordinate mappings (parent class) */ +#include "matrixmap.h" /* Interface definition for this class */ +#include "pal.h" /* SLALIB function definitions */ +#include "permmap.h" +#include "zoommap.h" +#include "unitmap.h" +#include "winmap.h" + +/* Error code definitions. */ +/* ----------------------- */ +#include "ast_err.h" /* AST error codes */ + +/* C header files. */ +/* --------------- */ +#include <ctype.h> +#include <math.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> + +/* Module Variables. */ +/* ================= */ + +/* Address of this static variable is used as a unique identifier for + member of this class. */ +static int class_check; +static const char *Form[3] = { "Full", "Diagonal", "Unit" }; /* Text values + used to represent storage form externally */ + +/* Pointers to parent class methods which are extended by this class. */ +static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * ); +static int *(* parent_mapsplit)( AstMapping *, int, const int *, AstMapping **, int * ); + + +#ifdef THREAD_SAFE +/* Define how to initialise thread-specific globals. */ +#define GLOBAL_inits \ + globals->Class_Init = 0; + +/* Create the function that initialises global data for this module. */ +astMAKE_INITGLOBALS(MatrixMap) + +/* Define macros for accessing each item of thread specific global data. */ +#define class_init astGLOBAL(MatrixMap,Class_Init) +#define class_vtab astGLOBAL(MatrixMap,Class_Vtab) + + +#include <pthread.h> + + +#else + + +/* Define the class virtual function table and its initialisation flag + as static variables. */ +static AstMatrixMapVtab class_vtab; /* Virtual function table */ +static int class_init = 0; /* Virtual function table initialised? */ + +#endif + +/* External Interface Function Prototypes. */ +/* ======================================= */ +/* The following functions have public prototypes only (i.e. no + protected prototypes), so we must provide local prototypes for use + within this module. */ +AstMatrixMap *astMatrixMapId_( int, int, int, const double [], const char *, ... ); + +/* Prototypes for Private Member Functions. */ +/* ======================================== */ +static AstMatrixMap *MatMat( AstMapping *, AstMapping *, int, int, int * ); +static AstMatrixMap *MatPerm( AstMatrixMap *, AstPermMap *, int, int, int, int * ); +static AstMatrixMap *MatZoom( AstMatrixMap *, AstZoomMap *, int, int, int * ); +static AstMatrixMap *MtrMult( AstMatrixMap *, AstMatrixMap *, int * ); +static AstMatrixMap *MtrRot( AstMatrixMap *, double, const double[], int * ); +static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * ); +static AstWinMap *MatWin2( AstMatrixMap *, AstWinMap *, int, int, int, int * ); +static double *InvertMatrix( int, int, int, double *, int * ); +static double Rate( AstMapping *, double *, int, int, int * ); +static int Equal( AstObject *, AstObject *, int * ); +static int FindString( int, const char *[], const char *, const char *, const char *, const char *, int * ); +static int Ustrcmp( const char *, const char *, int * ); +static int GetTranForward( AstMapping *, int * ); +static int GetIsLinear( AstMapping *, int * ); +static int GetTranInverse( AstMapping *, int * ); +static int CanSwap( AstMapping *, AstMapping *, int, int, int *, int * ); +static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * ); +static int PermOK( AstMapping *, int * ); +static int ScalingRowCol( AstMatrixMap *, int, int * ); +static void CompressMatrix( AstMatrixMap *, int * ); +static void Copy( const AstObject *, AstObject *, int * ); +static void Delete( AstObject *obj, int * ); +static void Dump( AstObject *, AstChannel *, int * ); +static void ExpandMatrix( AstMatrixMap *, int * ); +static void MatWin( AstMapping **, int *, int, int * ); +static void MatPermSwap( AstMapping **, int *, int, int * ); +static void PermGet( AstPermMap *, int **, int **, double **, int * ); +static void SMtrMult( int, int, int, const double *, double *, double*, int * ); +static int *MapSplit( AstMapping *, int, const int *, AstMapping **, int * ); + +/* Member functions. */ +/* ================= */ +static int CanSwap( AstMapping *map1, AstMapping *map2, int inv1, int inv2, + int *simpler, int *status ){ +/* +* Name: +* CanSwap + +* Purpose: +* Determine if two Mappings could be swapped. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* int CanSwap( AstMapping *map1, AstMapping *map2, int inv1, int inv2, +* int *simpler, int *status ) + +* Class Membership: +* MatrixMap member function + +* Description: +* This function returns a flag indicating if the pair of supplied +* Mappings could be replaced by an equivalent pair of Mappings from the +* same classes as the supplied pair, but in reversed order. Each pair +* of Mappings is considered to be compunded in series. The supplied +* Mapings are not changed in any way. + +* Parameters: +* map1 +* The Mapping to be applied first. +* map2 +* The Mapping to be applied second. +* inv1 +* The invert flag to use with map1. A value of zero causes the forward +* mapping to be used, and a non-zero value causes the inverse +* mapping to be used. +* inv2 +* The invert flag to use with map2. +* simpler +* Addresss of a location at which to return a flag indicating if +* the swapped Mappings would be intrinsically simpler than the +* original Mappings. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* 1 if the Mappings could be swapped, 0 otherwise. + +* Notes: +* - One of the supplied pair of Mappings must be a MatrixMap. +* - A value of 0 is returned if an error has already occurred, or if +* this function should fail for any reason. +*/ + +/* Local Variables: */ + AstMatrixMap *mat; /* Pointer to MatrixMap Mapping */ + AstMapping *nomat; /* Pointer to non-MatrixMap Mapping */ + const char *class1; /* Pointer to map1 class string */ + const char *class2; /* Pointer to map2 class string */ + const char *nomat_class; /* Pointer to non-MatrixMap class string */ + double *consts; /* Pointer to constants array */ + int *inperm; /* Pointer to input axis permutation array */ + int *outperm; /* Pointer to output axis permutation array */ + int i; /* Loop count */ + int invert[ 2 ]; /* Original invert flags */ + int nax; /* No. of in/out coordinates for the MatrixMap */ + int nin; /* No. of input coordinates for the PermMap */ + int nout; /* No. of output coordinates for the PermMap */ + int ret; /* Returned flag */ + +/* Check the global error status. */ + if ( !astOK ) return 0; + +/* Initialise */ + ret = 0; + *simpler = 0; + +/* Temporarily set the Invert attributes of both Mappings to the supplied + values. */ + invert[ 0 ] = astGetInvert( map1 ); + astSetInvert( map1, inv1 ); + + invert[ 1 ] = astGetInvert( map2 ); + astSetInvert( map2, inv2 ); + +/* Get the classes of the two mappings. */ + class1 = astGetClass( map1 ); + class2 = astGetClass( map2 ); + if( astOK ){ + +/* Get a pointer to the MatrixMap and non-MatrixMap Mappings. */ + if( !strcmp( class1, "MatrixMap" ) ){ + mat = (AstMatrixMap *) map1; + nomat = map2; + nomat_class = class2; + } else { + nomat = map1; + mat = (AstMatrixMap *) map2; + nomat_class = class1; + } + +/* Get the number of input axes for the MatrixMap. */ + nax = astGetNin( mat ); + +/* If it is a WinMap, the Mappings can be swapped. */ + if( !strcmp( nomat_class, "WinMap" ) ){ + ret = 1; + +/* If it is a PermMap, the Mappings can be swapped so long as: + 1) all links between input and output axes in the PermMap are + bi-directional. This does not preclude the existence of unconnected + axes, which do not have links (bi-directional or otherwise). + 2) The MatrixMap is square, and invertable. + 3) If the permMap is applied first, then each output of the PermMap + which is assigned a constant value must correspond to a "scaling" row + and column in the MatrixMap. I.e. if PermMap output axis "i" is + assigned a constant value, then row i and column i of the following + MatrixMap must contain only zeros, EXCEPT for the diagonal term (row + i, column i) which must be non-zero. If the Mappings are in the other + order, then the same applies to PermMap input axes assigned a constant + value. */ + +/* Check the other Mapping is a PermMap, and that the MatrixMap is square + and has an inverse. */ + } else if( !strcmp( nomat_class, "PermMap" ) && + nax == astGetNout( mat ) && ( mat->form == UNIT || + ( mat->i_matrix != NULL && + mat->f_matrix != NULL ) ) ) { + +/* Get the number of input and output coordinates for the PermMap. */ + nin = astGetNin( nomat ); + nout = astGetNout( nomat ); + +/* We need to know the axis permutation arrays and constants array for + the PermMap. */ + PermGet( (AstPermMap *) nomat, &outperm, &inperm, &consts, status ); + if( astOK ) { + +/* Indicate we can swap with the PermMap. */ + ret = 1; + +/* Check each output axis. If any links between axes are found which are + not bi-directional, indicate that we cannot swap with the PermMap. */ + for( i = 0; i < nout; i++ ){ + if( outperm[ i ] >= 0 && outperm[ i ] < nin ) { + if( inperm[ outperm[ i ] ] != i ) { + ret = 0; + break; + } + } + } + +/* Check each input axis. If any links between axes are found which are + not bi-directional, indicate that we cannot swap with the PermMap. */ + for( i = 0; i < nin; i++ ){ + if( inperm[ i ] >= 0 && inperm[ i ] < nout ) { + if( outperm[ inperm[ i ] ] != i ) { + ret = 0; + break; + } + } + } + +/* If the PermMap is suitable, check that any constant values fed from the + PermMap into the MatrixMap (in either forward or inverse direction) + are not changed by the MatrixMap. This requires the row and column for + each constant axis to be zeros, ecept for a value of 1.0 on the + diagonal. First deal with the cases where the PermMap is applied + first, so the outputs of the PermMap are fed into the MatrixMap in the + forward direction. */ + if( ret && ( nomat == map1 ) ) { + + if( nout != nax ){ + astError( AST__RDERR, "PermMap produces %d outputs, but the following" + "MatrixMap has %d inputs\n", status, nout, nax ); + ret = 0; + } + +/* Consider each output axis of the PermMap. */ + for( i = 0; i < nout && astOK ; i++ ) { + +/* If this PermMap output is assigned a constant... */ + if( outperm[ i ] < 0 || outperm[ i ] >= nin ) { + +/* Check the i'th row of the MatrixMap is all zero except for the i'th + column which must be non-zero. If not indicate that the MatrixMap cannot + swap with the PermMap and leave the loop. */ + if( !ScalingRowCol( mat, i, status ) ) { + ret = 0; + break; + } + } + } + } + +/* Now deal with the cases where the PermMap is applied second, so the inputs + of the PermMap are fed into the MatrixMap in the inverse direction. */ + if( ret && ( nomat == map2 ) ) { + + if( nin != nax ){ + astError( AST__RDERR, "Inverse PermMap produces %d inputs, but the " + "preceding MatrixMap has %d outputs\n", status, nin, nax ); + ret = 0; + } + +/* Consider each input axis of the PermMap. */ + for( i = 0; i < nin && astOK; i++ ){ + +/* If this PermMap input is assigned a constant (by the inverse Mapping)... */ + if( inperm[ i ] < 0 || inperm[ i ] >= nout ) { + +/* Check the i'th row of the MatrixMap is all zero except for the i'th + column which must be non-zero. If not indicate that the MatrixMap cannot + swap with the PermMap and leave the loop. */ + if( !ScalingRowCol( mat, i, status ) ) { + ret = 0; + break; + } + } + } + } + +/* If we can swap with the PermMap, the swapped Mappings may be + intrinsically simpler than the original mappings. */ + if( ret ) { + +/* If the PermMap precedes the WinMap, this will be the case if the PermMap + has more outputs than inputs. If the WinMap precedes the PermMap, this + will be the case if the PermMap has more inputs than outputs. */ + *simpler = ( nomat == map1 ) ? nout > nin : nin > nout; + } + +/* Free the axis permutation and constants arrays. */ + outperm = (int *) astFree( (void *) outperm ); + inperm = (int *) astFree( (void *) inperm ); + consts = (double *) astFree( (void *) consts ); + } + } + } + +/* Re-instate the original settings of the Invert attributes for the + supplied MatrixMaps. */ + astSetInvert( map1, invert[ 0 ] ); + astSetInvert( map2, invert[ 1 ] ); + +/* Return the answer. */ + return astOK ? ret : 0; +} + +static void CompressMatrix( AstMatrixMap *this, int *status ){ +/* +* Name: +* CompressMatrix + +* Purpose: +* If possible, reduce the amount of storage needed to store a MatrixMap. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* void CompressMatrix( AstMatrixMap *this, int *status ) + +* Class Membership: +* MatrixMap member function. + +* Description: +* The supplid MatrixMap is converted to its most compressed form +* (i.e no element values if it is a unit matrix, diagonal elements only +* if it is a diagonal matrix, or all elements otherwise). + +* Parameters: +* this +* A pointer to the MatrixMap to be compressed. +* status +* Pointer to the inherited status variable. + +*/ + +/* Local Variables: */ + double *a; /* Pointer to next element */ + double *colmax; /* Pointer to array holding column max values */ + double *fmat; /* Pointer to compressed forward matrix */ + double *rowmax; /* Pointer to array holding row max values */ + double mval; /* Matrix element value */ + int i; /* Loop count */ + int j; /* Loop count */ + int k; /* Loop count */ + int ncol; /* No. of columns in forward matrix */ + int ndiag; /* No. of diagonal elements in matrix */ + int new_form; /* Compressed storage form */ + int new_inv; /* New inverse requied? */ + int next_diag; /* Index of next diagonal element */ + int nrow; /* No. of rows in forward matrix */ + +/* Check the global error status. */ + if ( !astOK || !this ) return; + +/* Initialise variables to avoid "used of uninitialised variable" + messages from dumb compilers. */ + new_inv = 0; + +/* Get the dimensions of the forward matrix. */ + if( astGetInvert( this ) ){ + nrow = astGetNin( this ); + ncol = astGetNout( this ); + } else { + ncol = astGetNin( this ); + nrow = astGetNout( this ); + } + +/* Store the number of diagonal elements in the matrix. This is the + minimum of the number of rows and columns. */ + if( ncol < nrow ){ + ndiag = ncol; + } else { + ndiag = nrow; + } + +/* If the MatrixMap is already stored in UNIT form, it cannot be compressed + any further. */ + if( this->form == UNIT){ + return; + +/* Otherwise, if the MatrixMap is stored in DIAGONAL form, it could be + compressed into a UNIT MatrixMap if all the supplied element values are + one. */ + } else if( this->form == DIAGONAL ){ + new_form = UNIT; + for( i = 0; i < ndiag; i++ ){ + if( !astEQUAL( (this->f_matrix)[ i ], 1.0 ) ){ + new_form = DIAGONAL; + break; + } + } + +/* If it can be compressed, change the storage form and free the arrays + holding the diagonal element values. */ + if( new_form == UNIT ) { + this->f_matrix = (double *) astFree( (void *)( this->f_matrix ) ); + this->i_matrix = (double *) astFree( (void *)( this->i_matrix ) ); + this->form = UNIT; + } + +/* Otherwise, a full MatrixMap has been supplied, but this could be stored + in a unit or diagonal MatrixMap if the element values are appropriate. */ + } else { + new_form = FULL; + +/* Find the maximum absolute value in each column. Scale by + sqrt(DBL_EPSILON) to be come a lower limit for non-zero values. */ + colmax = astMalloc( ncol*sizeof( double ) ); + if( colmax ) { + for( j = 0; j < ncol; j++ ) { + colmax[ j ] = 0.0; + i = j; + for( k = 0; k < nrow; k++ ) { + mval = (this->f_matrix)[ i ]; + if( mval != AST__BAD ) { + mval = fabs( mval ); + if( mval > colmax[ j ] ) colmax[ j ] = mval; + } + i += ncol; + } + colmax[ j ] *= sqrt( DBL_EPSILON ); + } + } + +/* Find the maximum absolute value in each row. Scale by + sqrt(DBL_EPSILON) to be come a lower limit for non-zero values. */ + rowmax = astMalloc( nrow*sizeof( double ) ); + if( rowmax ) { + for( k = 0; k < nrow; k++ ) { + rowmax[ k ] = 0.0; + i = k*ncol; + for( j = 0; j < ncol; j++ ) { + mval = (this->f_matrix)[ i ]; + if( mval != AST__BAD ) { + mval = fabs( mval ); + if( mval > rowmax[ k ] ) rowmax[ k ] = mval; + } + i++; + } + rowmax[ k ] *= sqrt( DBL_EPSILON ); + } + } + +/* Check memory can be used */ + if( astOK ) { + +/* Initialise a flag indicating that the inverse matrix does not need to + be re-calculated. */ + new_inv = 0; + +/* Initially assume that the forward matrix is a unit matrix. */ + new_form = UNIT; + +/* Store a pointer to the next matrix element. */ + a = this->f_matrix; + +/* Loop through all the rows in the forward matrix array. */ + for( k = 0; k < nrow; k++ ) { + +/* Loop through all the elements in this column. */ + for( j = 0; j < ncol; j++, a++ ) { + +/* If this element is bad, use full form. */ + if( *a == AST__BAD ) { + new_form = FULL; + +/* Otherwise, if this is a diagonal term, check its value. If it is not one, + then the matrix cannot be a unit matrix, but it could still be a diagonal + matrix. */ + } else { + if( j == k ) { + if( *a != 1.0 && new_form == UNIT ) new_form = DIAGONAL; + +/* If this is not a diagonal element, and the element value is not zero, + then the matrix is not a diagonal matrix. Allow a tolerance of + SQRT(DBL_EPSILON) times the largest value in the same row or column as + the current matrix element. That is, an element must be insignificant + to both its row and its column to be considered as effectively zero. + Replace values less than this limit with zero. */ + } else { + mval = fabs( *a ); + if( mval <= rowmax[ k ] && + mval <= colmax[ j ] ) { + +/* If the element will change value, set a flag indicating that the inverse + matrix needs to be re-calculated. */ + if( *a != 0.0 ) new_inv = 1; + +/* Ensure this element value is zero. */ + *a = 0.0; + + } else { + new_form = FULL; + } + } + } + } + } + } + +/* Free memory. */ + colmax = astFree( colmax ); + rowmax = astFree( rowmax ); + +/* If it can be compressed into a UNIT MatrixMap, change the storage form and + free the arrays holding the element values. */ + if( new_form == UNIT ) { + this->f_matrix = (double *) astFree( (void *)( this->f_matrix ) ); + this->i_matrix = (double *) astFree( (void *)( this->i_matrix ) ); + this->form = UNIT; + +/* Otherwise, if it can be compressed into a DIAGONAL MatrixMap, copy the + diagonal elements from the full forward matrix into a newly allocated + array, use this array to replace the forward matrix array in the MatrixMap, + create a new inverse matrix, and change the storage form. */ + } else if( new_form == DIAGONAL ) { + fmat = astMalloc( sizeof(double)*(size_t)ndiag ); + if( fmat ){ + + next_diag = 0; + for( i = 0; i < ndiag; i++ ){ + fmat[ i ] = (this->f_matrix)[ next_diag ]; + next_diag += ncol + 1; + } + + (void) astFree( (void *) this->f_matrix ); + (void) astFree( (void *) this->i_matrix ); + + this->f_matrix = fmat; + this->i_matrix = InvertMatrix( DIAGONAL, nrow, ncol, fmat, status ); + this->form = DIAGONAL; + + } + +/* Calculate a new inverse matrix if necessary. */ + } else if( new_inv ) { + (void) astFree( (void *) this->i_matrix ); + this->i_matrix = InvertMatrix( FULL, nrow, ncol, this->f_matrix, status ); + } + } + + return; + +} + +static int Equal( AstObject *this_object, AstObject *that_object, int *status ) { +/* +* Name: +* Equal + +* Purpose: +* Test if two MatrixMaps are equivalent. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* int Equal( AstObject *this, AstObject *that, int *status ) + +* Class Membership: +* MatrixMap member function (over-rides the astEqual protected +* method inherited from the astMapping class). + +* Description: +* This function returns a boolean result (0 or 1) to indicate whether +* two MatrixMaps are equivalent. + +* Parameters: +* this +* Pointer to the first Object (a MatrixMap). +* that +* Pointer to the second Object. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* One if the MatrixMaps are equivalent, zero otherwise. + +* Notes: +* - A value of zero will be returned if this function is invoked +* with the global status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + AstMatrixMap *that; + AstMatrixMap *this; + double *that_matrix; + double *this_matrix; + int i; + int nin; + int nout; + int result; + +/* Initialise. */ + result = 0; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Obtain pointers to the two MatrixMap structures. */ + this = (AstMatrixMap *) this_object; + that = (AstMatrixMap *) that_object; + +/* Check the second object is a MatrixMap. We know the first is a + MatrixMap since we have arrived at this implementation of the virtual + function. */ + if( astIsAMatrixMap( that ) ) { + +/* Get the number of inputs and outputs and check they are the same for both. */ + nin = astGetNin( this ); + nout = astGetNout( this ); + if( astGetNout( that ) == nout && astGetNin( that ) == nin ) { + +/* Assume the MatrixMaps are equivalent. */ + result = 1; + +/* Ensure both MatrixMaps are stored in full form. */ + ExpandMatrix( this, status ); + ExpandMatrix( that, status ); + +/* Get pointers to the arrays holding the elements of the forward matrix + for both MatrixMaps. */ + if( astGetInvert( this ) ) { + this_matrix = this->i_matrix; + } else { + this_matrix = this->f_matrix; + } + + if( astGetInvert( that ) ) { + that_matrix = that->i_matrix; + } else { + that_matrix = that->f_matrix; + } + +/* If either of the above arrays is not available, try to get the inverse + matrix arrays. */ + if( !this_matrix || !that_matrix ) { + if( astGetInvert( this ) ) { + this_matrix = this->f_matrix; + } else { + this_matrix = this->i_matrix; + } + + if( astGetInvert( that ) ) { + that_matrix = that->f_matrix; + } else { + that_matrix = that->i_matrix; + } + } + +/* If both arrays are now available compare their elements. */ + if( this_matrix && that_matrix ) { + result = 1; + for( i = 0; i < nin*nout; i++ ) { + if( !astEQUAL( this_matrix[ i ], that_matrix[ i ] ) ){ + result = 0; + break; + } + } + } + +/* Ensure the supplied MatrixMaps are stored back in compressed form. */ + CompressMatrix( this, status ); + CompressMatrix( that, status ); + } + } + +/* If an error occurred, clear the result value. */ + if ( !astOK ) result = 0; + +/* Return the result, */ + return result; +} + +static void ExpandMatrix( AstMatrixMap *this, int *status ){ +/* +* Name: +* ExpandMatrix + +* Purpose: +* Ensure the MatrixMap is stored in full (non-compressed) form. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* void ExpandMatrix( AstMatrixMap *this, int *status ) + +* Class Membership: +* MatrixMap member function. + +* Description: +* If the supplid MatrixMap is stored in a compressed form (i.e no +* element values if it is a unit matrix, diagonal elements only +* if it is a diagonal matrix), it is expanded into a full MatrixMap +* in which all elements are stored. + +* Parameters: +* this +* A pointer to the MatrixMap to be expanded. +* status +* Pointer to the inherited status variable. + +*/ + +/* Local Variables: */ + double *fmat; /* Pointer to full forward matrix */ + double *imat; /* Pointer to full inverse matrix */ + int i; /* Loop count */ + int ncol; /* No. of columns in forward matrix */ + int ndiag; /* No. of diagonal elements in matrix */ + int nrow; /* No. of rows in forward matrix */ + +/* Check the global error status. Also return if the MatrixMap + pointer is null. */ + if ( !astOK || !this ) return; + +/* Return without action if the MatrixMap is already in full form. */ + if( this->form == FULL ) return; + +/* Get the dimensions of the forward matrix. */ + if( astGetInvert( this ) ){ + nrow = astGetNin( this ); + ncol = astGetNout( this ); + } else { + ncol = astGetNin( this ); + nrow = astGetNout( this ); + } + +/* Store the number of diagonal elements. */ + if( nrow > ncol ){ + ndiag = ncol; + } else { + ndiag = nrow; + } + +/* Allocate arrays to hold the full forward and inverse matrices. */ + fmat = (double *) astMalloc( sizeof( double )*(size_t)( nrow*ncol ) ); + imat = (double *) astMalloc( sizeof( double )*(size_t)( nrow*ncol ) ); + if( imat && fmat ){ + +/* Fill them both with zeros. */ + for( i = 0; i < nrow*ncol; i++ ) { + fmat[ i ] = 0.0; + imat[ i ] = 0.0; + } + +/* If a unit MatrixMap was supplied, put ones on the diagonals. */ + if( this->form == UNIT ){ + for( i = 0; i < ndiag; i++ ) { + fmat[ i*( ncol + 1 ) ] = 1.0; + imat[ i*( nrow + 1 ) ] = 1.0; + } + +/* If a diagonal MatrixMap was supplied, copy the diagonal terms from + the supplied MatrixMap. */ + } else if( this->form == DIAGONAL ){ + for( i = 0; i < ndiag; i++ ) { + fmat[ i*( ncol + 1 ) ] = (this->f_matrix)[ i ]; + imat[ i*( nrow + 1 ) ] = (this->i_matrix)[ i ]; + } + } + +/* Free any existing arrays in the MatrixMap and store the new ones. */ + (void) astFree( (void *) this->f_matrix ); + (void) astFree( (void *) this->i_matrix ); + + this->f_matrix = fmat; + this->i_matrix = imat; + +/* Update the storage form. */ + this->form = FULL; + +/* If either of the new matrices could not be allocated, ensure that + both have been freed. */ + } else { + fmat = (double *) astFree( (void *) fmat ); + imat = (double *) astFree( (void *) imat ); + } + + return; + +} + +static int FindString( int n, const char *list[], const char *test, + const char *text, const char *method, + const char *class, int *status ){ +/* +* Name: +* FindString + +* Purpose: +* Find a given string within an array of character strings. + +* Type: +* Private function. + +* Synopsis: +* #include "matrix.h" +* int FindString( int n, const char *list[], const char *test, +* const char *text, const char *method, const char *class, int *status ) + +* Class Membership: +* MatrixMap method. + +* Description: +* This function identifies a supplied string within a supplied +* array of valid strings, and returns the index of the string within +* the array. The test option may not be abbreviated, but case is +* insignificant. + +* Parameters: +* n +* The number of strings in the array pointed to be "list". +* list +* A pointer to an array of legal character strings. +* test +* A candidate string. +* text +* A string giving a description of the object, parameter, +* attribute, etc, to which the test value refers. +* This is only for use in constructing error messages. It should +* start with a lower case letter. +* method +* Pointer to a string holding the name of the calling method. +* This is only for use in constructing error messages. +* class +* Pointer to a string holding the name of the supplied object class. +* This is only for use in constructing error messages. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* The index of the identified string within the supplied array, starting +* at zero. + +* Notes: +* - A value of -1 is returned if an error has already occurred, or +* if this function should fail for any reason (for instance if the +* supplied option is not specified in the supplied list). + +*/ + +/* Local Variables: */ + int ret; /* The returned index */ + +/* Check global status. */ + if( !astOK ) return -1; + +/* Compare the test string with each element of the supplied list. Leave + the loop when a match is found. */ + for( ret = 0; ret < n; ret++ ) { + if( !Ustrcmp( test, list[ ret ], status ) ) break; + } + +/* Report an error if the supplied test string does not match any element + in the supplied list. */ + if( ret >= n ) { + astError( AST__RDERR, "%s(%s): Illegal value '%s' supplied for %s.", status, + method, class, test, text ); + ret = -1; + } + +/* Return the answer. */ + return ret; +} + +static int GetIsLinear( AstMapping *this_mapping, int *status ){ +/* +* Name: +* GetIsLinear + +* Purpose: +* Return the value of the IsLinear attribute for a MatrixMap. + +* Type: +* Private function. + +* Synopsis: +* #include "mapping.h" +* void GetIsLinear( AstMapping *this, int *status ) + +* Class Membership: +* MatrixMap member function (over-rides the protected astGetIsLinear +* method inherited from the Mapping class). + +* Description: +* This function returns the value of the IsLinear attribute for a +* Frame, which is always one. + +* Parameters: +* this +* Pointer to the MatrixMap. +* status +* Pointer to the inherited status variable. +*/ + return 1; +} + +static int Ustrcmp( const char *a, const char *b, int *status ){ +/* +* Name: +* Ustrncmp + +* Purpose: +* A case blind version of strcmp. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* int Ustrcmp( const char *a, const char *b ) + +* Class Membership: +* MatrixMap member function. + +* Description: +* Returns 0 if there are no differences between the two strings, and 1 +* otherwise. Comparisons are case blind. + +* Parameters: +* a +* Pointer to first string. +* b +* Pointer to second string. + +* Returned Value: +* Zero if the strings match, otherwise one. + +* Notes: +* - This function does not consider the sign of the difference between +* the two strings, whereas "strcmp" does. +* - This function attempts to execute even if an error has occurred. + +*/ + +/* Local Variables: */ + const char *aa; /* Pointer to next "a" character */ + const char *bb; /* Pointer to next "b" character */ + int ret; /* Returned value */ + +/* Initialise the returned value to indicate that the strings match. */ + ret = 0; + +/* Initialise pointers to the start of each string. */ + aa = a; + bb = b; + +/* Loop round each character. */ + while( 1 ){ + +/* We leave the loop if either of the strings has been exhausted. */ + if( !(*aa ) || !(*bb) ){ + +/* If one of the strings has not been exhausted, indicate that the + strings are different. */ + if( *aa || *bb ) ret = 1; + +/* Break out of the loop. */ + break; + +/* If neither string has been exhausted, convert the next characters to + upper case and compare them, incrementing the pointers to the next + characters at the same time. If they are different, break out of the + loop. */ + } else { + + if( toupper( (int) *(aa++) ) != toupper( (int) *(bb++) ) ){ + ret = 1; + break; + } + + } + + } + +/* Return the result. */ + return ret; + +} + +void astInitMatrixMapVtab_( AstMatrixMapVtab *vtab, const char *name, int *status ) { +/* +*+ +* Name: +* astInitMatrixMapVtab + +* Purpose: +* Initialise a virtual function table for a MatrixMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "matrixmap.h" +* void astInitMatrixMapVtab( AstMatrixMapVtab *vtab, const char *name ) + +* Class Membership: +* MatrixMap vtab initialiser. + +* Description: +* This function initialises the component of a virtual function +* table which is used by the MatrixMap class. + +* Parameters: +* vtab +* Pointer to the virtual function table. The components used by +* all ancestral classes will be initialised if they have not already +* been initialised. +* name +* Pointer to a constant null-terminated character string which contains +* the name of the class to which the virtual function table belongs (it +* is this pointer value that will subsequently be returned by the Object +* astClass function). +*- +*/ + +/* Local Variables: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstObjectVtab *object; /* Pointer to Object component of Vtab */ + AstMappingVtab *mapping; /* Pointer to Mapping component of Vtab */ + +/* Check the local error status. */ + if ( !astOK ) return; + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(NULL); + +/* Initialize the component of the virtual function table used by the + parent class. */ + astInitMappingVtab( (AstMappingVtab *) vtab, name ); + +/* Store a unique "magic" value in the virtual function table. This + will be used (by astIsAMatrixMap) to determine if an object belongs + to this class. We can conveniently use the address of the (static) + class_check variable to generate this unique value. */ + vtab->id.check = &class_check; + vtab->id.parent = &(((AstMappingVtab *) vtab)->id); + +/* Initialise member function pointers. */ +/* ------------------------------------ */ +/* Store pointers to the member functions (implemented here) that provide + virtual methods for this class. */ + vtab->MtrRot = MtrRot; + vtab->MtrMult = MtrMult; + +/* Save the inherited pointers to methods that will be extended, and + replace them with pointers to the new member functions. */ + object = (AstObjectVtab *) vtab; + mapping = (AstMappingVtab *) vtab; + + parent_transform = mapping->Transform; + mapping->Transform = Transform; + + parent_mapsplit = mapping->MapSplit; + mapping->MapSplit = MapSplit; + +/* Store replacement pointers for methods which will be over-ridden by + new member functions implemented here. */ + object->Equal = Equal; + mapping->GetIsLinear = GetIsLinear; + mapping->GetTranForward = GetTranForward; + mapping->GetTranInverse = GetTranInverse; + mapping->MapMerge = MapMerge; + mapping->Rate = Rate; + +/* Declare the destructor and copy constructor. */ + astSetDelete( (AstObjectVtab *) vtab, Delete ); + astSetCopy( (AstObjectVtab *) vtab, Copy ); + +/* Declare the class dump function. */ + astSetDump( vtab, Dump, "MatrixMap", "Matrix transformation" ); + +/* If we have just initialised the vtab for the current class, indicate + that the vtab is now initialised, and store a pointer to the class + identifier in the base "object" level of the vtab. */ + if( vtab == &class_vtab ) { + class_init = 1; + astSetVtabClassIdentifier( vtab, &(vtab->id) ); + } +} + + +static double *InvertMatrix( int form, int nrow, int ncol, double *matrix, int *status ){ +/* +* Name: +* InvertMatrix + +* Purpose: +* Invert a suplied matrix. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* double *InvertMatrix( int form, int nrow, int ncol, double *matrix, int *status ) + +* Class Membership: +* MatrixMap member function. + +* Description: +* This function returns a pointer to a matrix holding the inverse of +* the supplied matrix, or a NULL pointer if the inverse is not defined. +* The memory to store the inverse matrix is allocated internally, and +* should be freed using astFree when no longer required. +* +* The correspondence between a full matrix and its inverse is only +* unique if the matrix is square, and so a NULL pointer is returned if +* the supplied matrix is not square. + +* Parameters: +* form +* The form of the MatrixMap; UNIT, DIAGONAL or FULL. +* nrow +* Number of rows in the supplied matrix. +* ncol +* Number of columns in the supplied matrix. +* matrix +* A pointer to the input matrix. Elements should be stored in row +* order (i.e. (row 1,column 1 ), (row 1,column 2 )... (row 2,column 1), +* etc). +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Pointer to the output matrix. + +* Notes: +* - A NULL pointer is returned if a unit matrix is supplied. +* - A NULL pointer will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +* - No error is reported if the inverse is not defined. +*/ + +/* Local Variables: */ + double det; /* Determinant of supplied matrix */ + double mval; /* Matrix element value */ + double *out; /* Pointer to returned inverse matrix */ + double *vector; /* Pointer to vector used by palDmat */ + int i; /* Matrix element number */ + int *iw; /* Pointer to workspace used by palDmat */ + int nel; /* No. of elements in square matrix */ + int ndiag; /* No. of diagonal elements */ + int ok; /* Zero if any bad matrix values found */ + int sing; /* Zero if matrix is not singular */ + +/* Check the global error status. */ + if ( !astOK ) return NULL; + +/* Return a NULL pointer if the input matrix is NULL. */ + if( !matrix ) return NULL; + +/* If a unit matrix map has been supplied, return NULL. */ + if( form == UNIT ){ + return NULL; + +/* If a diagonal matrix has been supplied, allocate an array to hold + the diagonal terms of the inverse matrix. Store the reciprocal + of the input matrix diagonal terms in it. If any of the input diagonal + terms are zero or BAD, set the associated elements of the inverse matrix + BAD. */ + } else if( form == DIAGONAL ){ + if( nrow > ncol ) { + ndiag = ncol; + } else { + ndiag = nrow; + } + + out = (double *) astMalloc( sizeof( double )*(size_t)ndiag ); + + if( out ) { + for( i = 0; i < ndiag; i++ ) { + mval = matrix[ i ]; + if( mval != 0.0 && mval != AST__BAD ){ + out[ i ] = 1.0/mval; + } else { + out[ i ] = AST__BAD; + } + } + } + +/* If a full matrix has been supplied, initialise the returned pointer. */ + } else { + out = NULL; + +/* Check that the matrix is square. */ + if( nrow == ncol ){ + +/* Find the number of elements in the matrix. */ + nel = nrow*ncol; + +/* See if there are any bad values in the matrix. */ + ok = 1; + for ( i=0; i<nel; i++ ) { + if ( matrix[i] == AST__BAD ) { + ok = 0; + break; + } + } + +/* Only continue if there are no bad matrix values. */ + if( ok ) { + +/* Take a copy of the supplied matrix */ + out = (double *) astStore( NULL, (void *) matrix, + astSizeOf( (void *) matrix ) ); + +/* The SLALIB function which inverts the matrix also applies the inverse + matrix to a vector. We are not interested in the vector in this + instance, but we still have to provide one for SLALIB to use. Allocate + memory for the vector. */ + vector = (double *) astMalloc( sizeof(double)*(size_t) nrow ); + +/* If it was allocated succesfully, fill it with zeros. */ + if( astOK ){ + for ( i=0; i<nrow; i++ ) vector[i] = 0.0; + +/* Obtain work space and attempt to invert the matrix using SLALIB, then + free the work space. */ + iw = (int *) astMalloc( sizeof(int)*(size_t) nrow ); + if( astOK ) palDmat( nrow, out, vector, &det, &sing, iw ); + iw = (int *) astFree( (void *) iw ); + + } + +/* If the matrix could not be inverted, free the memory used to hold the + square matrix, and return the NULL pointer. */ + if ( !astOK || sing != 0 ){ + out = (double *) astFree( (void *) out ); + } + +/* Free the memory used to hold the vector. */ + vector = (double *) astFree( (void *) vector ); + } + } + } + +/* Return the pointer. */ + + return out; + +} + +static int MapMerge( AstMapping *this, int where, int series, int *nmap, + AstMapping ***map_list, int **invert_list, int *status ) { +/* +* Name: +* MapMerge + +* Purpose: +* Simplify a sequence of Mappings containing a MatrixMap. + +* Type: +* Private function. + +* Synopsis: +* #include "mapping.h" +* int MapMerge( AstMapping *this, int where, int series, int *nmap, +* AstMapping ***map_list, int **invert_list, int *status ) + +* Class Membership: +* MatrixMap method (over-rides the protected astMapMerge method +* inherited from the Mapping class). + +* Description: +* This function attempts to simplify a sequence of Mappings by +* merging a nominated MatrixMap in the sequence with its neighbours, +* so as to shorten the sequence if possible. +* +* In many cases, simplification will not be possible and the +* function will return -1 to indicate this, without further +* action. +* +* In most cases of interest, however, this function will either +* attempt to replace the nominated MatrixMap with a Mapping which it +* considers simpler, or to merge it with the Mappings which +* immediately precede it or follow it in the sequence (both will +* normally be considered). This is sufficient to ensure the +* eventual simplification of most Mapping sequences by repeated +* application of this function. +* +* In some cases, the function may attempt more elaborate +* simplification, involving any number of other Mappings in the +* sequence. It is not restricted in the type or scope of +* simplification it may perform, but will normally only attempt +* elaborate simplification in cases where a more straightforward +* approach is not adequate. + +* Parameters: +* this +* Pointer to the nominated MatrixMap which is to be merged with +* its neighbours. This should be a cloned copy of the MatrixMap +* pointer contained in the array element "(*map_list)[where]" +* (see below). This pointer will not be annulled, and the +* MatrixMap it identifies will not be modified by this function. +* where +* Index in the "*map_list" array (below) at which the pointer +* to the nominated MatrixMap resides. +* series +* A non-zero value indicates that the sequence of Mappings to +* be simplified will be applied in series (i.e. one after the +* other), whereas a zero value indicates that they will be +* applied in parallel (i.e. on successive sub-sets of the +* input/output coordinates). +* nmap +* Address of an int which counts the number of Mappings in the +* sequence. On entry this should be set to the initial number +* of Mappings. On exit it will be updated to record the number +* of Mappings remaining after simplification. +* map_list +* Address of a pointer to a dynamically allocated array of +* Mapping pointers (produced, for example, by the astMapList +* method) which identifies the sequence of Mappings. On entry, +* the initial sequence of Mappings to be simplified should be +* supplied. +* +* On exit, the contents of this array will be modified to +* reflect any simplification carried out. Any form of +* simplification may be performed. This may involve any of: (a) +* removing Mappings by annulling any of the pointers supplied, +* (b) replacing them with pointers to new Mappings, (c) +* inserting additional Mappings and (d) changing their order. +* +* The intention is to reduce the number of Mappings in the +* sequence, if possible, and any reduction will be reflected in +* the value of "*nmap" returned. However, simplifications which +* do not reduce the length of the sequence (but improve its +* execution time, for example) may also be performed, and the +* sequence might conceivably increase in length (but normally +* only in order to split up a Mapping into pieces that can be +* more easily merged with their neighbours on subsequent +* invocations of this function). +* +* If Mappings are removed from the sequence, any gaps that +* remain will be closed up, by moving subsequent Mapping +* pointers along in the array, so that vacated elements occur +* at the end. If the sequence increases in length, the array +* will be extended (and its pointer updated) if necessary to +* accommodate any new elements. +* +* Note that any (or all) of the Mapping pointers supplied in +* this array may be annulled by this function, but the Mappings +* to which they refer are not modified in any way (although +* they may, of course, be deleted if the annulled pointer is +* the final one). +* invert_list +* Address of a pointer to a dynamically allocated array which, +* on entry, should contain values to be assigned to the Invert +* attributes of the Mappings identified in the "*map_list" +* array before they are applied (this array might have been +* produced, for example, by the astMapList method). These +* values will be used by this function instead of the actual +* Invert attributes of the Mappings supplied, which are +* ignored. +* +* On exit, the contents of this array will be updated to +* correspond with the possibly modified contents of the +* "*map_list" array. If the Mapping sequence increases in +* length, the "*invert_list" array will be extended (and its +* pointer updated) if necessary to accommodate any new +* elements. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* If simplification was possible, the function returns the index +* in the "map_list" array of the first element which was +* modified. Otherwise, it returns -1 (and makes no changes to the +* arrays supplied). + +* Notes: +* - A value of -1 will be returned if this function is invoked +* with the global error status set, or if it should fail for any +* reason. +*/ + +/* Local Variables: */ + AstMapping **maplt; /* New mappings list pointer */ + AstMapping *map2; /* Pointer to replacement Mapping */ + AstMapping *mc[2]; /* Copies of supplied Mappings to swap */ + AstMapping *newmap; /* Pointer to replacement MatrixMap */ + AstMapping *smc0; /* Simplied Mapping */ + AstMapping *smc1; /* Simplied Mapping */ + AstMatrixMap *mm; /* Pointer to supplied MatrixMap */ + const char *class1; /* Pointer to first Mapping class string */ + const char *class2; /* Pointer to second Mapping class string */ + const char *nclass; /* Pointer to neighbouring Mapping class */ + double *b; /* Pointer to scale terms */ + double *new_mat; /* Pointer to elements of new MatrixMap */ + double factor; /* Zoom factor for new ZoomMap */ + int *invlt; /* New invert flags list pointer */ + int do1; /* Would a backward swap make a simplification? */ + int do2; /* Would a forward swap make a simplification? */ + int i1; /* Index of first MatrixMap to merge */ + int i2; /* Index of last MatrixMap to merge */ + int i; /* Loop counter */ + int ic[2]; /* Copies of supplied invert flags to swap */ + int invert; /* Should the inverted Mapping be used? */ + int j; /* Loop counter */ + int nin; /* Number of input coordinates for MatrixMap */ + int nmapt; /* No. of Mappings in list */ + int nout; /* Number of output coordinates for MatrixMap */ + int nstep1; /* No. of Mappings backwards to next mergable Mapping */ + int nstep2; /* No. of Mappings forward to next mergable Mapping */ + int result; /* Result value to return */ + int swaphi; /* Can MatrixMap be swapped with higher neighbour? */ + int swaplo; /* Can MatrixMap be swapped with lower neighbour? */ + int zoom; /* Can MatrixMap be replaced by a ZoomMap? */ + +/* Initialise. */ + result = -1; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Initialise variables to avoid "used of uninitialised variable" + messages from dumb compilers. */ + i1 = 0; + i2 = 0; + +/* Get the Invert attribute for the specified mapping. */ + invert = astGetInvert( ( *map_list )[ where ] ); + +/* Get the number of input and output axes for the MatrixMap. Swap these + if the supplied invert flag is not the same as the Invert attribute of + the Mapping. */ + if( ( invert && !( *invert_list )[ where ] ) || + ( !invert && ( *invert_list )[ where ] ) ) { + nout = astGetNin( ( *map_list )[ where ] ); + nin = astGetNout( ( *map_list )[ where ] ); + + } else { + nin = astGetNin( ( *map_list )[ where ] ); + nout = astGetNout( ( *map_list )[ where ] ); + } + +/* First of all, see if the MatrixMap can be replaced by a simpler Mapping, + without reference to the neighbouring Mappings in the list. */ +/* ======================================================================*/ + map2 = NULL; + mm = (AstMatrixMap *) ( *map_list )[ where ]; + +/* If the MatrixMap is a square unit matrix, it can be replaced by a + UnitMap. */ + if( mm->form == UNIT && nin == nout ){ + map2 = (AstMapping *) astUnitMap( nin, "", status ); + +/* If the MatrixMap is a square diagonal matrix with equal diagonal + terms, then it can be replaced by a ZoomMap, so long as the + diagonal elements are not all zero. */ + } else if( mm->form == DIAGONAL && nin == nout && + mm->f_matrix && mm->i_matrix && + (mm->f_matrix)[ 0 ] != AST__BAD ){ + zoom = 1; + b = mm->f_matrix + 1; + for( i = 1; i < nin; i++ ){ + if( !astEQUAL( *b, *( b - 1 ) ) ){ + zoom = 0; + break; + } + b++; + } + + if( zoom ){ + if( ( *invert_list )[ where ] ){ + factor = (mm->i_matrix)[ 0 ]; + } else { + factor = (mm->f_matrix)[ 0 ]; + } + + if( factor != 0.0 ){ + map2 = (AstMapping *) astZoomMap( nin, factor, "", status ); + } + } + +/* If the MatrixMap is a full matrix but all off-diagonal elements are + zero, it can be replaced by a diagonal MatrixMap. */ + } else if( mm->form == FULL && nin == nout && mm->f_matrix ){ + new_mat = astMalloc( sizeof( double )*nin ); + b = mm->f_matrix; + for( i = 0; i < nin && new_mat; i++ ){ + for( j = 0; j < nout; j++,b++ ){ + if( i == j ) { + new_mat[ i ] = *b; + } else if( *b != 0.0 ) { + new_mat = astFree( new_mat ); + break; + } + } + } + + if( new_mat ) { + map2 = (AstMapping *) astMatrixMap( nin, nout, 1, new_mat, "", + status ); + new_mat = astFree( new_mat ); + } + } + +/* If the MatrixMap can be replaced, annul the MatrixMap pointer in the + list and replace it with the new Mapping pointer, and indicate that the + forward transformation of the returned Mapping should be used. */ + if( map2 ){ + (void) astAnnul( ( *map_list )[ where ] ); + ( *map_list )[ where ] = map2; + ( *invert_list )[ where ] = 0; + +/* Return the index of the first modified element. */ + result = where; + +/* If the MatrixMap itself could not be simplified, see if it can be merged + with the Mappings on either side of it in the list. */ +/*==========================================================================*/ + } else { + +/* Store the classes of the neighbouring Mappings in the list. */ + class1 = ( where > 0 ) ? astGetClass( ( *map_list )[ where - 1 ] ) : NULL; + class2 = ( where < *nmap - 1 ) ? astGetClass( ( *map_list )[ where + 1 ] ) : NULL; + +/* In series. */ +/* ========== */ + if ( series ) { + +/* We first look to see if the MatrixMap can be merged with one of its + neighbours, resulting in a reduction of one in the number of Mappings + in the list. MatrixMaps can merge directly with another MatrixMap, a + ZoomMap, an invertable PermMap, or a UnitMap. */ + if( class1 && ( !strcmp( class1, "MatrixMap" ) || + !strcmp( class1, "ZoomMap" ) || + !strcmp( class1, "PermMap" ) || + !strcmp( class1, "UnitMap" ) ) ){ + nclass = class1; + i1 = where - 1; + i2 = where; + + } else if( class2 && ( !strcmp( class2, "MatrixMap" ) || + !strcmp( class2, "ZoomMap" ) || + !strcmp( class2, "PermMap" ) || + !strcmp( class2, "UnitMap" ) ) ){ + nclass = class2; + i1 = where; + i2 = where + 1; + + } else { + nclass = NULL; + } + +/* Only some PermMaps can be merged with (those which have consistent + forward and inverse mappings). If this is not one of them, set nclass + NULL to indicate this. */ + if( nclass && !strcmp( nclass, "PermMap" ) && + !PermOK( ( *map_list )[ (i1==where)?i2:i1 ], status ) ) nclass = NULL; + +/* If the MatrixMap is diagonal it can also merge with a WinMap. */ + if( !nclass && mm->form == DIAGONAL) { + if( class1 && ( !strcmp( class1, "WinMap" ) ) ){ + nclass = class1; + i1 = where - 1; + i2 = where; + + } else if( class2 && ( !strcmp( class2, "WinMap" ) ) ){ + nclass = class2; + i1 = where; + i2 = where + 1; + + } + } + +/* If the MatrixMap can merge with one of its neighbours, create the merged + Mapping. */ + if( nclass ){ + + if( !strcmp( nclass, "MatrixMap" ) ){ + newmap = (AstMapping *) MatMat( ( *map_list )[ i1 ], ( *map_list )[ i2 ], + ( *invert_list )[ i1 ], ( *invert_list )[ i2 ], status ); + invert = 0; + + } else if( !strcmp( nclass, "ZoomMap" ) ){ + if( i1 == where ){ + newmap = (AstMapping *) MatZoom( (AstMatrixMap *)( *map_list )[ i1 ], + (AstZoomMap *)( *map_list )[ i2 ], + ( *invert_list )[ i1 ], ( *invert_list )[ i2 ], status ); + } else { + newmap = (AstMapping *) MatZoom( (AstMatrixMap *)( *map_list )[ i2 ], + (AstZoomMap *)( *map_list )[ i1 ], + ( *invert_list )[ i2 ], ( *invert_list )[ i1 ], status ); + } + invert = 0; + + } else if( !strcmp( nclass, "PermMap" ) ){ + if( i1 == where ){ + newmap = (AstMapping *) MatPerm( (AstMatrixMap *)( *map_list )[ i1 ], + (AstPermMap *)( *map_list )[ i2 ], + ( *invert_list )[ i1 ], ( *invert_list )[ i2 ], 1, status ); + } else { + newmap = (AstMapping *) MatPerm( (AstMatrixMap *)( *map_list )[ i2 ], + (AstPermMap *)( *map_list )[ i1 ], + ( *invert_list )[ i2 ], ( *invert_list )[ i1 ], 0, status ); + } + invert = 0; + + } else if( !strcmp( nclass, "WinMap" ) ){ + if( i1 == where ){ + newmap = (AstMapping *) MatWin2( (AstMatrixMap *)( *map_list )[ i1 ], + (AstWinMap *)( *map_list )[ i2 ], + ( *invert_list )[ i1 ], ( *invert_list )[ i2 ], 1, status ); + } else { + newmap = (AstMapping *) MatWin2( (AstMatrixMap *)( *map_list )[ i2 ], + (AstWinMap *)( *map_list )[ i1 ], + ( *invert_list )[ i2 ], ( *invert_list )[ i1 ], 0, status ); + } + invert = 0; + + } else { + newmap = astClone( ( *map_list )[ where ] ); + invert = ( *invert_list )[ where ]; + } + +/* If succesfull... */ + if( astOK ){ + +/* Annul the first of the two Mappings, and replace it with the merged + MatrixMap. Also set the invert flag. */ + (void) astAnnul( ( *map_list )[ i1 ] ); + ( *map_list )[ i1 ] = newmap; + ( *invert_list )[ i1 ] = invert; + +/* Annul the second of the two Mappings, and shuffle down the rest of the + list to fill the gap. */ + (void) astAnnul( ( *map_list )[ i2 ] ); + for ( i = i2 + 1; i < *nmap; i++ ) { + ( *map_list )[ i - 1 ] = ( *map_list )[ i ]; + ( *invert_list )[ i - 1 ] = ( *invert_list )[ i ]; + } + +/* Clear the vacated element at the end. */ + ( *map_list )[ *nmap - 1 ] = NULL; + ( *invert_list )[ *nmap - 1 ] = 0; + +/* Decrement the Mapping count and return the index of the first + modified element. */ + ( *nmap )--; + result = i1; + + } + +/* If the MatrixMap could not merge directly with either of its neighbours, + we consider whether it would be worthwhile to swap the MatrixMap with + either of its neighbours. This can only be done for certain classes + of Mapping (WinMaps and some PermMaps), and will usually require both + Mappings to be modified (unless they are commutative). The advantage of + swapping the order of the Mappings is that it may result in the MatrixMap + being adjacent to a Mapping with which it can merge directly on the next + invocation of this function, thus reducing the number of Mappings + in the list. */ + } else { + +/* Set a flag if we could swap the MatrixMap with its higher neighbour. "do2" + is returned if swapping the Mappings would simplify either of the Mappings. */ + if( where + 1 < *nmap ){ + swaphi = CanSwap( ( *map_list )[ where ], + ( *map_list )[ where + 1 ], + ( *invert_list )[ where ], + ( *invert_list )[ where + 1 ], &do2, status ); + } else { + swaphi = 0; + do2 = 0; + } + +/* If so, step through each of the Mappings which follow the MatrixMap, + looking for a Mapping with which the MatrixMap could merge directly. Stop + when such a Mapping is found, or if a Mapping is found with which the + MatrixMap could definitely not swap. Note the number of Mappings which + separate the MatrixMap from the Mapping with which it could merge (if + any). */ + nstep2 = -1; + if( swaphi ){ + for( i2 = where + 1; i2 < *nmap; i2++ ){ + +/* See if we can merge with this Mapping. If so, note the number of steps + between the two Mappings and leave the loop. */ + nclass = astGetClass( ( *map_list )[ i2 ] ); + if( !strcmp( nclass, "MatrixMap" ) || + !strcmp( nclass, "ZoomMap" ) || + ( !strcmp( nclass, "PermMap" ) && PermOK( ( *map_list )[ i2 ], status ) ) || + !strcmp( nclass, "UnitMap" ) ) { + nstep2 = i2 - where - 1; + break; + } + +/* If there is no chance that we can swap with this Mapping, leave the loop + with -1 for the number of steps to indicate that no merging is possible. + MatrixMaps can swap with WinMaps and some permmaps. */ + if( strcmp( nclass, "WinMap" ) && + strcmp( nclass, "PermMap" ) ) { + break; + } + + } + + } + +/* Do the same working forward from the MatrixMap towards the start of the map + list. */ + if( where > 0 ){ + swaplo = CanSwap( ( *map_list )[ where - 1 ], + ( *map_list )[ where ], + ( *invert_list )[ where - 1 ], + ( *invert_list )[ where ], &do1, status ); + } else { + swaplo = 0; + do1 = 0; + } + + nstep1 = -1; + if( swaplo ){ + for( i1 = where - 1; i1 >= 0; i1-- ){ + + nclass = astGetClass( ( *map_list )[ i1 ] ); + if( !strcmp( nclass, "MatrixMap" ) || + ( !strcmp( nclass, "PermMap" ) && PermOK( ( *map_list )[ i1 ], status ) ) || + !strcmp( nclass, "ZoomMap" ) || + !strcmp( nclass, "UnitMap" ) ) { + nstep1 = where - 1 - i1; + break; + } + + if( strcmp( nclass, "WinMap" ) && + strcmp( nclass, "PermMap" ) ) { + break; + } + + } + + } + +/* Choose which neighbour to swap with so that the MatrixMap moves towards the + nearest Mapping with which it can merge. */ + if( do1 || ( + nstep1 != -1 && ( nstep2 == -1 || nstep2 > nstep1 ) ) ){ + nclass = class1; + i1 = where - 1; + i2 = where; + } else if( do2 || nstep2 != -1 ){ + nclass = class2; + i1 = where; + i2 = where + 1; + } else { + nclass = NULL; + } + +/* If there is a target Mapping in the list with which the MatrixMap could + merge, consider replacing the supplied Mappings with swapped Mappings to + bring the MatrixMap closer to the target Mapping. */ + if( nclass ){ + +/* Swap the Mappings. */ + if (!strcmp( nclass, "WinMap" ) ){ + MatWin( (*map_list) + i1, (*invert_list) + i1, where - i1, status ); + + } else if( !strcmp( nclass, "PermMap" ) ){ + MatPermSwap( (*map_list) + i1, (*invert_list) + i1, where - i1, status ); + } + +/* And then merge them. */ + if( where == i1 && where + 1 < *nmap ) { /* Merging upwards */ + map2 = astClone( (*map_list)[ where + 1 ] ); + nmapt = *nmap - where - 1; + maplt = *map_list + where + 1; + invlt = *invert_list + where + 1; + + (void) astMapMerge( map2, 0, series, &nmapt, &maplt, &invlt ); + map2 = astAnnul( map2 ); + *nmap = where + 1 + nmapt; + + } else if( where - 2 >= 0 ) { /* Merging downwards */ + map2 = astClone( (*map_list)[ where - 2 ] ); + nmapt = *nmap - where + 2; + maplt = *map_list + where - 2 ; + invlt = *invert_list + where - 2; + + (void) astMapMerge( map2, 0, series, &nmapt, &maplt, &invlt ); + map2 = astAnnul( map2 ); + *nmap = where - 2 + nmapt; + } + + result = i1; + +/* If there is no Mapping available for merging, it may still be + advantageous to swap with a neighbour because the swapped Mapping may + be simpler than the original Mappings. For instance, a PermMap may + strip rows of the MatrixMap leaving only a UnitMap. */ + } else if( swaphi || swaplo ) { + +/* Try swapping with each possible neighbour in turn. */ + for( i = 0; i < 2; i++ ) { + +/* Set up the class and pointers for the mappings to be swapped, first + the lower neighbour, then the upper neighbour. */ + if( i == 0 && swaplo ){ + nclass = class1; + i1 = where - 1; + i2 = where; + + } else if( i == 1 && swaphi ){ + nclass = class2; + i1 = where; + i2 = where + 1; + + } else { + nclass = NULL; + } + +/* If we have a Mapping to swap with... */ + if( nclass ) { + +/* Take copies of the Mapping and Invert flag arrays so we do not change + the supplied values. */ + mc[ 0 ] = (AstMapping *) astCopy( ( (*map_list) + i1 )[0] ); + mc[ 1 ] = (AstMapping *) astCopy( ( (*map_list) + i1 )[1] ); + ic[ 0 ] = ( (*invert_list) + i1 )[0]; + ic[ 1 ] = ( (*invert_list) + i1 )[1]; + +/* Swap these Mappings. */ + if( !strcmp( nclass, "WinMap" ) ){ + MatWin( mc, ic, where - i1, status ); + } else if( !strcmp( nclass, "PermMap" ) ){ + MatPermSwap( mc, ic, where - i1, status ); + } + +/* If neither of the swapped Mappings can be simplified further, then there + is no point in swapping the Mappings, so just annul the map copies. */ + smc0 = astSimplify( mc[0] ); + smc1 = astSimplify( mc[1] ); + + if( astGetClass( smc0 ) == astGetClass( mc[0] ) && + astGetClass( smc1 ) == astGetClass( mc[1] ) ) { + + mc[ 0 ] = (AstMapping *) astAnnul( mc[ 0 ] ); + mc[ 1 ] = (AstMapping *) astAnnul( mc[ 1 ] ); + +/* If one or both of the swapped Mappings could be simplified, then annul + the supplied Mappings and return the swapped mappings, storing the index + of the first modified Mapping. */ + } else { + (void ) astAnnul( ( (*map_list) + i1 )[0] ); + (void ) astAnnul( ( (*map_list) + i1 )[1] ); + + ( (*map_list) + i1 )[0] = mc[ 0 ]; + ( (*map_list) + i1 )[1] = mc[ 1 ]; + + ( (*invert_list) + i1 )[0] = ic[ 0 ]; + ( (*invert_list) + i1 )[1] = ic[ 1 ]; + + result = i1; + break; + } + +/* Annul the simplied Mappings */ + smc0 = astAnnul( smc0 ); + smc1 = astAnnul( smc1 ); + + } + } + } + } + } + } + +/* Return the result. */ + return result; +} + +static int *MapSplit( AstMapping *this_map, int nin, const int *in, AstMapping **map, int *status ){ +/* +* Name: +* MapSplit + +* Purpose: +* Create a Mapping representing a subset of the inputs of an existing +* MatrixMap. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* int *MapSplit( AstMapping *this, int nin, const int *in, AstMapping **map, int *status ) + +* Class Membership: +* MatrixMap method (over-rides the protected astMapSplit method +* inherited from the Mapping class). + +* Description: +* This function creates a new Mapping by picking specified inputs from +* an existing MatrixMap. This is only possible if the specified inputs +* correspond to some subset of the MatrixMap outputs. That is, there +* must exist a subset of the MatrixMap outputs for which each output +* depends only on the selected MatrixMap inputs, and not on any of the +* inputs which have not been selected. In addition, outputs that are +* not in this subset must not depend on any selected inputs. If these +* conditions are not met by the supplied MatrixMap, then a NULL Mapping +* is returned. + +* Parameters: +* this +* Pointer to the MatrixMap to be split (the MatrixMap is not actually +* modified by this function). +* nin +* The number of inputs to pick from "this". +* in +* Pointer to an array of indices (zero based) for the inputs which +* are to be picked. This array should have "nin" elements. If "Nin" +* is the number of inputs of the supplied MatrixMap, then each element +* should have a value in the range zero to Nin-1. +* map +* Address of a location at which to return a pointer to the new +* Mapping. This Mapping will have "nin" inputs (the number of +* outputs may be different to "nin"). A NULL pointer will be +* returned if the supplied MatrixMap has no subset of outputs which +* depend only on the selected inputs. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* A pointer to a dynamically allocated array of ints. The number of +* elements in this array will equal the number of outputs for the +* returned Mapping. Each element will hold the index of the +* corresponding output in the supplied MatrixMap. The array should be +* freed using astFree when no longer needed. A NULL pointer will +* be returned if no output Mapping can be created. + +* Notes: +* - If this function is invoked with the global error status set, +* or if it should fail for any reason, then NULL values will be +* returned as the function value and for the "map" pointer. +*/ + +/* Local Variables: */ + AstMatrixMap *this; /* Pointer to MatrixMap structure */ + double *mat; /* Pointer to matrix for supplied MatrixMap */ + double *pmat; /* Pointer to row start in returned matrix */ + double *prow; /* Pointer to row start in supplied matrix */ + double *rmat; /* Pointer to matrix for returned MatrixMap */ + double el; /* Next element value in supplied matrix */ + int *result; /* Pointer to returned array */ + int good; /* Would new matrix contain any good values/ */ + int i; /* Loop count */ + int icol; /* Column index within supplied MatrixMap */ + int iel; /* Index of next element from the input matrix */ + int irow; /* Row index within supplied MatrixMap */ + int isel; /* Does output depend on any selected inputs? */ + int ncol; /* Number of columns (inputs) in supplied MatrixMap */ + int nout; /* Number of outputs in returned MatrixMap */ + int nrow; /* Number of rows (outputs) in supplied MatrixMap */ + int ok; /* Are input indices OK? */ + int sel; /* Does any output depend on selected inputs? */ + int unsel; /* Does any output depend on unselected inputs? */ + +/* Initialise */ + result = NULL; + *map = NULL; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Invoke the parent astMapSplit method to see if it can do the job. */ + result = (*parent_mapsplit)( this_map, nin, in, map, status ); + +/* If not, we provide a special implementation here. */ + if( !result ) { + +/* Get a pointer to the MatrixMap structure. */ + this = (AstMatrixMap *) this_map; + +/* Get the number of inputs and outputs. */ + ncol = astGetNin( this ); + nrow = astGetNout( this ); + +/* Check the supplied input indices are usable. */ + ok = 1; + for( i = 0; i < nin; i++ ) { + if( in[ i ] < 0 || in[ i ] >= ncol ) { + ok = 0; + break; + } + } + + if( ok ) { + +/* Ensure the MatrixMap is stored in full form. */ + ExpandMatrix( this, status ); + +/* Allocate the largest array that could be necessary to hold the + returned array of Mapping outputs. */ + result = astMalloc( sizeof(int)*(size_t) nrow ); + +/* Allocate the largest array that could be necessary to hold the + matrix representing the returned MatrixMap. */ + rmat = astMalloc( sizeof(double)*(size_t) (nrow*ncol) ); + +/* Get the matrix which defines the current forward transformation. This + takes into account whether the MatrixMap has been inverted or not. */ + if( astGetInvert( this ) ) { + mat = this->i_matrix; + } else { + mat = this->f_matrix; + } + +/* We cannot create the require Mapping if the matrix is undefined. */ + if( !mat || !astOK ) { + ok = 0; + nout = 0; + good = 0; + +/* Otherwise, loop round all the rows in the matrix. */ + } else { + nout = 0; + good = 0; + pmat = rmat; + iel = 0; + for( irow = 0; irow < nrow; irow++ ) { + +/* Indicate that this output (i.e. row of the matrix) depends on neither + selected nor unselected inputs as yet. */ + sel = 0; + unsel = 0; + +/* Save a pointer to the first element of this row in the MatrixMap + matrix. */ + prow = mat + iel; + +/* Loop round all the elements in the current row of the matrix. */ + for( icol = 0; icol < ncol; icol++ ) { + +/* If this element is non-zero and non-bad, then output "irow" depends on + input "icol". */ + el = mat[ iel++ ]; + if( el != 0.0 && el != AST__BAD ) { + +/* Is input "icol" one of the selected inputs? */ + isel = 0; + for( i = 0; i < nin; i++ ) { + if( in[ i ] == icol ) { + isel = 1; + break; + } + } + +/* If so, note that this output depends on selected inputs. Otherwise note + it depends on unselected inputs. */ + if( isel ) { + sel = 1; + } else { + unsel = 1; + } + } + } + +/* If this output depends only on selected inputs, we can include it in + the returned Mapping.*/ + if( sel && !unsel ) { + +/* Store the index of the output within the original MatrixMap. */ + result[ nout ] = irow; + +/* Increment the number of outputs in the returned Mapping. */ + nout++; + +/* Copy the elements of the current matrix row which correspond to the + selected inputs into the new matrix. */ + for( i = 0; i < nin; i++ ) { + if( astISGOOD( prow[ in[ i ] ] ) ) { + *(pmat++) = prow[ in[ i ] ]; + good = 1; + } + } + } + +/* If this output depends on a selected input, but also depends on an + unselected input, we cannot split the MatrixMap. */ + if( sel && unsel ) { + ok = 0; + break; + } + } + } + + +/* If the returned Mapping can be created, create it. */ + if( ok && nout > 0 && good ) { + *map = (AstMapping *) astMatrixMap( nin, nout, 0, rmat, "", status ); + +/* Otherwise, free the returned array. */ + } else { + result = astFree( result ); + } + +/* Free resources. */ + rmat = astFree( rmat ); + +/* Re-compress the supplied MatrixMap. */ + CompressMatrix( this, status ); + } + } + +/* Free returned resources if an error has occurred. */ + if( !astOK ) { + result = astFree( result ); + *map = astAnnul( *map ); + } + +/* Return the list of output indices. */ + return result; +} + +static AstMatrixMap *MatMat( AstMapping *map1, AstMapping *map2, int inv1, + int inv2, int *status ){ +/* +* Name: +* MatMat + +* Purpose: +* Create a merged MatrixMap from two supplied MatrixMaps. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* AstMatrixMap *MatMat( AstMapping *map1, AstMapping *map2, int inv1, +* int inv2, int *status ) + +* Class Membership: +* MatrixMap member function + +* Description: +* This function creates a new MatrixMap which performs a mapping +* equivalent to applying the two supplied MatrixMaps in series, in the +* directions specified by the "invert" flags (the Invert attributes of +* the supplied MatrixMaps are ignored). + +* Parameters: +* map1 +* A pointer to the MatrixMap to apply first. +* map2 +* A pointer to the MatrixMap to apply second. +* inv1 +* The invert flag to use with map1. A value of zero causes the forward +* mapping to be used, and a non-zero value causes the inverse +* mapping to be used. +* inv2 +* The invert flag to use with map2. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Pointer to the new MatrixMap. + +* Notes: +* - The forward direction of the returned MatrixMap is equivalent to the +* combined effect of the two supplied MatrixMap, operating in the +* directions specified by "inv1" and "inv2". +* - A null pointer will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + AstMatrixMap *result; /* Pointer to output MatrixMap */ + int invert[ 2 ]; /* Original invert flags */ + +/* Check the global error status. */ + if ( !astOK ) return NULL; + +/* Initialise the returned pointer. */ + result = NULL; + +/* Temporarily set their Invert attributes to the supplied values. */ + invert[ 0 ] = astGetInvert( map1 ); + astSetInvert( map1, inv1 ); + + invert[ 1 ] = astGetInvert( map2 ); + astSetInvert( map2, inv2 ); + +/* Create a new MatrixMap by multiplying them together. */ + result = astMtrMult( (AstMatrixMap *) map1, (AstMatrixMap *) map2 ); + +/* Re-instate the original settings of the Invert attributes for the + supplied MatrixMaps. */ + astSetInvert( map1, invert[ 0 ] ); + astSetInvert( map2, invert[ 1 ] ); + +/* If an error has occurred, annull the returned MatrixMap. */ + if( !astOK ) result = astAnnul( result ); + +/* Return a pointer to the output MatrixMap. */ + return result; +} + +static AstMatrixMap *MatPerm( AstMatrixMap *mm, AstPermMap *pm, int minv, + int pinv, int mat1, int *status ){ +/* +* Name: +* MatPerm + +* Purpose: +* Create a MatrixMap by merging a MatrixMap and a PermMap. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* AstMatrixMap *MatPerm( AstMatrixMap *mm, AstPermMap *pm, int minv, +* int pinv, int mat1, int *status ) + +* Class Membership: +* MatrixMap member function + +* Description: +* This function creates a new MatrixMap which performs a mapping +* equivalent to applying the two supplied Mappings in series in the +* directions specified by the "invert" flags (the Invert attributes of +* the supplied MatrixMaps are ignored). + +* Parameters: +* mm +* A pointer to the MatrixMap. +* pm +* A pointer to the PermMap. +* minv +* The invert flag to use with mm. A value of zero causes the forward +* mapping to be used, and a non-zero value causes the inverse +* mapping to be used. +* pinv +* The invert flag to use with pm. +* mat1 +* If non-zero, then "mm" is applied first followed by "pm". Otherwise, +* "pm" is applied first followed by "mm". +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Pointer to the new MatrixMap. + +* Notes: +* - The forward direction of the returned MatrixMap is equivalent to the +* combined effect of the two supplied Mappings, operating in the +* directions specified by "pinv" and "minv". +* - A null pointer will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + AstMatrixMap *mm2; /* Pointer to an intermediate MatrixMap */ + AstMatrixMap *result; /* Pointer to output MatrixMap */ + AstPointSet *pset1; /* Pointer to a PointSet holding unpermuted unit vectors */ + AstPointSet *pset2; /* Pointer to a PointSet holding permuted unit vectors */ + double *matrix; /* Pointer to a matrix representing the PermMap */ + double *p; /* Pointer to next matrix element */ + double **ptr1; /* Pointer to the data in pset1 */ + double **ptr2; /* Pointer to the data in pset2 */ + int i; /* Axis index */ + int j; /* Point index */ + int nax; /* No. of axes in the PermMap */ + int old_minv; /* Original setting of MatrixMap Invert attribute */ + int old_pinv; /* Original setting of PermMap Invert attribute */ + +/* Check the global error status. */ + if ( !astOK ) return NULL; + +/* Initialise the returned pointer. */ + result = NULL; + +/* Temporarily set the Invert attributes of both Mappings to the supplied + values. */ + old_minv = astGetInvert( mm ); + astSetInvert( mm, minv ); + + old_pinv = astGetInvert( pm ); + astSetInvert( pm, pinv ); + +/* Get the number of axes in the PermMap. The PermMap will have the same + number of input and output axes because a check has already been made on + it to ensure that this is so (in function PermOK). */ + nax = astGetNin( pm ); + +/* We first represent the PermMap as a MatrixMap containing elements with + values zero or one. Each row of this matrix is obtained by transforming a + unit vector along each axis using the inverse PermMap. Allocate memory + to hold the matrix array, and create a PointSet holding the unit + vectors. */ + matrix = (double *) astMalloc( sizeof( double )*(size_t)( nax*nax ) ); + + pset1 = astPointSet( nax, nax, "", status ); + ptr1 = astGetPoints( pset1 ); + + pset2 = astPointSet( nax, nax, "", status ); + ptr2 = astGetPoints( pset2 ); + + if( astOK ){ + for( i = 0; i < nax; i++ ){ + for( j = 0; j < nax; j++ ) ptr1[ i ][ j ] = 0.0; + ptr1[ i ][ i ] = 1.0; + } + +/* Transform these unit vectors using the inverse PermMap. */ + (void) astTransform( pm, pset1, 0, pset2 ); + +/* Copy the transformed vectors into the matrix array. */ + p = matrix; + for( j = 0; j < nax; j++ ){ + for( i = 0; i < nax; i++ ) *(p++) = ptr2[ i ][ j ]; + } + +/* Create a MatrixMap holding this array. */ + mm2 = astMatrixMap( nax, nax, 0, matrix, "", status ); + +/* Create a new MatrixMap equal to the product of the supplied MatrixMap + and the MatrixMap just created from the PermMap. */ + if( mat1 ){ + result = astMtrMult( mm, mm2 ); + } else { + result = astMtrMult( mm2, mm ); + } + +/* Free everything. */ + mm2 = astAnnul( mm2 ) ; + } + + pset2 = astAnnul( pset2 ); + pset1 = astAnnul( pset1 ); + matrix = (double *) astFree( (void *) matrix ); + +/* Re-instate the original settings of the Invert attribute for the + supplied Mappings. */ + astSetInvert( mm, old_minv ); + astSetInvert( pm, old_pinv ); + +/* If an error has occurred, annull the returned MatrixMap. */ + if( !astOK ) result = astAnnul( result ); + +/* Return a pointer to the output MatrixMap. */ + return result; +} + +static void MatPermSwap( AstMapping **maps, int *inverts, int imm, int *status ){ +/* +* Name: +* MatPermSwap + +* Purpose: +* Swap a PermMap and a MatrixMap. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* void MatPermSwap( AstMapping **maps, int *inverts, int imm ) + +* Class Membership: +* MatrixMap member function + +* Description: +* A list of two Mappings is supplied containing a PermMap and a +* MatrixMap. These Mappings are annulled, and replaced with +* another pair of Mappings consisting of a PermMap and a MatrixMap +* in the opposite order. These Mappings are chosen so that their +* combined effect is the same as the original pair of Mappings. + +* Parameters: +* maps +* A pointer to an array of two Mapping pointers. +* inverts +* A pointer to an array of two invert flags. +* imm +* The index within "maps" of the MatrixMap. + +* Notes: +* - There are restictions on the sorts of PermMaps which can be +* swapped with a MatrixMap -- see function CanSwap. It is assumed +* that the supplied MatrixMap and PermMap satisfy these requirements. + +*/ + +/* Local Variables: */ + AstMatrixMap *mm; /* Pointer to the supplied MatrixMap */ + AstMatrixMap *mmnew; /* Pointer to new MatrixMap */ + AstMatrixMap *smmnew; /* Pointer to new simplified MatrixMap */ + AstPermMap *pm; /* Pointer to the supplied PermMap */ + AstPermMap *pmnew; /* Pointer to new PermMap */ + AstPermMap *spmnew; /* Pointer to new simplified PermMap */ + double *consts; /* Pointer to constants array */ + double *matrix; /* Supplied array of matrix elements */ + double *out_el; /* Pointer to next element of new MatrixMap */ + double *out_mat; /* Matrix elements for new MatrixMap */ + double c; /* Constant */ + double matel; /* Matrix element */ + int *inperm; /* Pointer to input axis permutation array */ + int *outperm; /* Pointer to output axis permutation array */ + int col; /* Index of matrix column */ + int i; /* Axis count */ + int k; /* Axis count */ + int nin; /* No. of axes in supplied PermMap */ + int nout; /* No. of axes in returned PermMap */ + int old_pinv; /* Invert value for the supplied PermMap */ + int row; /* Index of matrix row */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Initialise variables to avoid "used of uninitialised variable" + messages from dumb compilers. */ + mmnew = NULL; + pmnew = NULL; + +/* Store pointers to the supplied PermMap and the MatrixMap. */ + pm = (AstPermMap *) maps[ 1 - imm ]; + mm = (AstMatrixMap *) maps[ imm ]; + +/* Temporarily set the Invert attribute of the supplied PermMap to the + supplied value. */ + old_pinv = astGetInvert( pm ); + astSetInvert( pm, inverts[ 1 - imm ] ); + +/* Ensure the MatrixMap is stored in full form. */ + ExpandMatrix( mm, status ); + +/* Store a pointer to the required array of matrix elements. */ + if( inverts[ imm ] ) { + matrix = mm->i_matrix; + } else { + matrix = mm->f_matrix; + } + +/* Get the number of input and output axes of the PermMap. */ + nin = astGetNin( pm ); + nout = astGetNout( pm ); + +/* Allocate memory to hold the matrix elements for the swapped MatrixMap. + The number of rows and olumns in the new matrix must equal the number of + input or output axes for the PermMap, depending on whether the PermMap + or MatrixMap is applied first. */ + if( imm == 0 ) { + out_mat = (double *) astMalloc( sizeof( double )*(size_t)( nout*nout ) ); + } else { + out_mat = (double *) astMalloc( sizeof( double )*(size_t)( nin*nin ) ); + } + +/* We need to know the axis permutation arrays and constants array for + the PermMap. */ + PermGet( pm, &outperm, &inperm, &consts, status ); + if( astOK ) { + +/* First deal with cases where the MatrixMap is applied first. */ + if( imm == 0 ) { + +/* Consider each output axis of the PermMap. */ + for( i = 0; i < nout; i++ ) { + +/* If this output is connected to one of the input axes... */ + row = outperm[ i ]; + if( row >= 0 && row < nin ) { + +/* Permute the row of the supplied matrix which feeds the corresponding + PermMap input axis (i.e. axis outperm[k] ) using the forward PermMap. + Store zeros for any output axes which are assigned constants. This forms + row i of the new MatrixMap. */ + out_el = out_mat + nout*i; + for( k = 0; k < nout; k++ ){ + col = outperm[ k ]; + if( col >= 0 && col < nin ) { + *(out_el++) = *( matrix + nin*row + col ); + } else { + *(out_el++) = 0.0; + } + } + +/* If this output is asigned a constant value, use a "diagonal" vector for + row i of the new MatrixMap (i.e. all zeros except for a 1.0 in column + i ). */ + } else { + out_el = out_mat + nout*i; + for( k = 0; k < nout; k++ ) { + if( k != i ) { + *(out_el++) = 0.0; + } else { + *(out_el++) = 1.0; + } + } + } + } + +/* Create the new MatrixMap. */ + mmnew = astMatrixMap( nout, nout, 0, out_mat, "", status ); + +/* Any PermMap inputs which are assigned a constant value need to be + changed now, since they will no longer be scaled by the inverse + MatrixMap. CanSwap ensures that the inverse MatrixMap produces a + simple scaling for constant axes, so we change the PermMap constant + to be the constant AFTER scaling by the inverse MatrixMap. + + Consider each input axis of the PermMap. */ + for( i = 0; i < nin; i++ ) { + +/* If this input is assigned a constant value... */ + if( inperm[ i ] < 0 ) { + +/* Divide the supplied constant value by the corresponding diagonal term + in the supplied MatrixMap. */ + c = consts[ -inperm[ i ] - 1 ]; + if( c != AST__BAD ) { + matel = matrix[ i*( nin + 1 ) ]; + if( matel != 0.0 && matel != AST__BAD ) { + consts[ -inperm[ i ] - 1 ] /= matel; + } else { + consts[ -inperm[ i ] - 1 ] = AST__BAD; + } + } + } + } + +/* Now deal with cases where the PermMap is applied first. */ + } else { + +/* Consider each input axis of the PermMap. */ + for( i = 0; i < nin; i++ ) { + +/* If this input is connected to one of the output axes... */ + row = inperm[ i ]; + if( row >= 0 && row < nout ) { + +/* Permute the row of the supplied matrix which feeds the corresponding + PermMap output axis (i.e. axis inperm[k] ) using the inverse PermMap. + Store zeros for any input axes which are assigned constants. This forms + row i of the new MatrixMap. */ + out_el = out_mat + nin*i; + for( k = 0; k < nin; k++ ){ + col = inperm[ k ]; + if( col >= 0 && col < nout ) { + *(out_el++) = *( matrix + nout*row + col ); + } else { + *(out_el++) = 0.0; + } + } + +/* If this input is asigned a constant value, use a "diagonal" vector for + row i of the new MatrixMap (i.e. all zeros except for a 1.0 in column + i ). */ + } else { + out_el = out_mat + nin*i; + for( k = 0; k < nin; k++ ) { + if( k != i ) { + *(out_el++) = 0.0; + } else { + *(out_el++) = 1.0; + } + } + } + } + +/* Create the new MatrixMap. */ + mmnew = astMatrixMap( nin, nin, 0, out_mat, "", status ); + +/* Any PermMap outputs which are assigned a constant value need to be + changed now, since they will no longer be scaled by the forward + MatrixMap. CanSwap ensures that the forward MatrixMap produces a + simple scaling for constant axes, so we change the PermMap constant + to be the constant AFTER scaling by the forward MatrixMap. + + Consider each output axis of the PermMap. */ + for( i = 0; i < nout; i++ ) { + +/* If this output is assigned a constant value... */ + if( outperm[ i ] < 0 ) { + +/* Multiple the supplied constant value by the corresponding diagonal term in + the supplied MatrixMap. */ + c = consts[ -outperm[ i ] - 1 ]; + if( c != AST__BAD ) { + matel = matrix[ i*( nout + 1 ) ]; + if( matel != AST__BAD ) { + consts[ -outperm[ i ] - 1 ] *= matel; + } else { + consts[ -outperm[ i ] - 1 ] = AST__BAD; + } + } + } + } + } + +/* Create a new PermMap (since the constants may have changed). */ + pmnew = astPermMap( nin, inperm, nout, outperm, consts, "", status ); + +/* Free the axis permutation and constants arrays. */ + outperm = (int *) astFree( (void *) outperm ); + inperm = (int *) astFree( (void *) inperm ); + consts = (double *) astFree( (void *) consts ); + } + +/* Free the memory used to hold the new matrix elements. */ + out_mat = (double *) astFree( (void *) out_mat ); + +/* Ensure the supplied MatrixMap is stored back in compressed form. */ + CompressMatrix( mm, status ); + +/* Re-instate the original value of the Invert attribute of the supplied + PermMap. */ + astSetInvert( pm, old_pinv ); + + if( astOK ) { + +/* Annul the supplied PermMap. */ + (void) astAnnul( pm ); + +/* Simplify the returned Mappings. */ + spmnew = astSimplify( pmnew ); + pmnew = astAnnul( pmnew ); + + smmnew = astSimplify( mmnew ); + mmnew = astAnnul( mmnew ); + +/* Store a pointer to the new PermMap in place of the supplied MatrixMap. This + PermMap should be used in its forward direction. */ + maps[ imm ] = (AstMapping *) spmnew; + inverts[ imm ] = astGetInvert( spmnew ); + +/* Annul the supplied matrixMap. */ + (void) astAnnul( mm ); + +/* Store a pointer to the new MatrixMap. This MatrixMap should be used in + its forward direction. */ + maps[ 1 - imm ] = (AstMapping *) smmnew; + inverts[ 1 - imm ] = astGetInvert( smmnew ); + } + +/* Return. */ + return; +} + +static void MatWin( AstMapping **maps, int *inverts, int imm, int *status ){ +/* +* Name: +* MatWin + +* Purpose: +* Swap a WinMap and a MatrixMap. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* void MatWin( AstMapping **maps, int *inverts, int imm, int *status ) + +* Class Membership: +* WinMap member function + +* Description: +* A list of two Mappings is supplied containing a WinMap and a +* MatrixMap. These Mappings are annulled, and replaced with +* another pair of Mappings consisting of a WinMap and a MatrixMap +* in the opposite order. These Mappings are chosen so that their +* combined effect is the same as the original pair of Mappings. +* The scale factors in the returned WinMap are always unity (i.e. +* the differences in scaling get absorbed into the returned +* MatrixMap). + +* Parameters: +* maps +* A pointer to an array of two Mapping pointers. +* inverts +* A pointer to an array of two invert flags. +* imm +* The index within "maps" of the MatrixMap. +* status +* Pointer to the inherited status variable. + +*/ + +/* Local Variables: */ + AstMatrixMap *m1; /* Pointer to Diagonal scale factor MatrixMap */ + AstMatrixMap *m2; /* Pointer to returned MatrixMap */ + AstMatrixMap *sm2; /* Pointer to simplified returned MatrixMap */ + AstMatrixMap *mm; /* Pointer to the supplied MatrixMap */ + AstPointSet *pset1; /* Shift terms from supplied WinMap */ + AstPointSet *pset2; /* Shift terms for returned WinMap */ + AstWinMap *w1; /* Pointer to the returned WinMap */ + AstWinMap *sw1; /* Pointer to the simplified returned WinMap */ + AstWinMap *wm; /* Pointer to the supplied WinMap */ + double **ptr1; /* Pointer to pset1 data */ + double **ptr2; /* Pointer to pset2 data */ + double *a; /* Array of shift terms from supplied WinMap */ + double *aa; /* Pointer to next shift term */ + double *b; /* Array of scale terms from supplied WinMap */ + double *bb; /* Pointer to next scale term */ + int i; /* Axis count */ + int nin; /* No. of axes in supplied WinMap */ + int nout; /* No. of axes in returned WinMap */ + int old_minv; /* Invert value for the supplied MatrixMap */ + int old_winv; /* Invert value for the supplied WinMap */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Store pointers to the supplied WinMap and the MatrixMap. */ + wm = (AstWinMap *) maps[ 1 - imm ]; + mm = (AstMatrixMap *) maps[ imm ]; + +/* Temporarily set the Invert attribute of the supplied Mappings to the + supplied values. */ + old_winv = astGetInvert( wm ); + astSetInvert( wm, inverts[ 1 - imm ] ); + + old_minv = astGetInvert( mm ); + astSetInvert( mm, inverts[ imm ] ); + +/* Get copies of the shift and scale terms used by the WinMap. This + also returns the number of axes in the WinMap. */ + nin = astWinTerms( wm, &a, &b ); + +/* Create a diagonal MatrixMap holding the scale factors from the + supplied WinMap. */ + m1 = astMatrixMap( nin, nin, 1, b, "", status ); + +/* Create a PointSet holding a single position given by the shift terms + in the supplied WinMap. */ + pset1 = astPointSet( 1, nin, "", status ); + ptr1 = astGetPoints( pset1 ); + if( astOK ){ + aa = a; + for( i = 0; i < nin; i++ ) ptr1[ i ][ 0 ] = *(aa++); + } + +/* First deal with cases when the WinMap is applied first, followed by + the MatrixMap. */ + if( imm == 1 ){ + +/* Multiply the diagonal matrix holding the WinMap scale factors by the + supplied matrix. The resulting MatrixMap is the one to return in the + map list. */ + m2 = astMtrMult( m1, mm ); + +/* Transform the position given by the shift terms from the supplied + WinMap using the supplied MatrixMap to get the shift terms for + the returned WinMap. */ + pset2 = astTransform( mm, pset1, 1, NULL ); + +/* Now deal with cases when the MatrixMap is applied first, followed by + the WinMap. */ + } else { + +/* Multiply the supplied MatrixMap by the diagonal matrix holding scale + factors from the supplied WinMap. The resulting MatrixMap is the one to + return in the map list. */ + m2 = astMtrMult( mm, m1 ); + +/* Transform the position given by the shift terms from the supplied + WinMap using the inverse of the returned MatrixMap to get the shift + terms for the returned WinMap. */ + pset2 = astTransform( m2, pset1, 0, NULL ); + + } + +/* Re-instate the original value of the Invert attributes of the supplied + Mappings. */ + astSetInvert( wm, old_winv ); + astSetInvert( mm, old_minv ); + +/* Get pointers to the shift terms for the returned WinMap. */ + ptr2 = astGetPoints( pset2 ); + +/* Create the returned WinMap, initially with undefined corners. The number of + axes in the WinMap must equal the number of shift terms. */ + nout = astGetNcoord( pset2 ); + w1 = astWinMap( nout, NULL, NULL, NULL, NULL, "", status ); + +/* If succesful, store the scale and shift terms in the WinMap. The scale + terms are always unity. */ + if( astOK ){ + bb = w1->b; + aa = w1->a; + for( i = 0; i < nout; i++ ) { + *(bb++) = 1.0; + *(aa++) = ptr2[ i ][ 0 ]; + } + +/* Replace the supplied Mappings and invert flags with the ones found + above. Remember that the order of the Mappings is now swapped */ + (void) astAnnul( maps[ 0 ] ); + (void) astAnnul( maps[ 1 ] ); + + sw1 = astSimplify( w1 ); + w1 = astAnnul( w1 ); + + maps[ imm ] = (AstMapping *) sw1; + inverts[ imm ] = astGetInvert( sw1 ); + + sm2 = astSimplify( m2 ); + m2 = astAnnul( m2 ); + + maps[ 1 - imm ] = (AstMapping *) sm2; + inverts[ 1 - imm ] = astGetInvert( sm2 ); + + } + +/* Annul the MatrixMap and PointSet holding the scale and shift terms from the + supplied WinMap. */ + m1 = astAnnul( m1 ); + pset1 = astAnnul( pset1 ); + pset2 = astAnnul( pset2 ); + +/* Free the copies of the scale and shift terms from the supplied WinMap. */ + b = (double *) astFree( (void *) b ); + a = (double *) astFree( (void *) a ); + +/* Return. */ + return; +} + +static AstWinMap *MatWin2( AstMatrixMap *mm, AstWinMap *wm, int minv, + int winv, int mat1, int *status ){ +/* +* Name: +* MatWin2 + +* Purpose: +* Create a WinMap by merging a diagonal MatrixMap and a WinMap. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* AstWinMap *MatWin2( AstMatrixMap *mm, AstWinMap *wm, int minv, +* int winv, int mat1, int *status ) + +* Class Membership: +* MatrixMap member function + +* Description: +* This function creates a new WinMap which performs a mapping +* equivalent to applying the two supplied Mappings in series in the +* directions specified by the "invert" flags (the Invert attributes of +* the supplied MatrixMaps are ignored), in the order specified by +* "mat1". + +* Parameters: +* mm +* A pointer to the MatrixMap. Assumed to be diagonal. +* wm +* A pointer to the WinMap. +* minv +* The invert flag to use with mm. A value of zero causes the forward +* mapping to be used, and a non-zero value causes the inverse +* mapping to be used. +* winv +* The invert flag to use with wm. +* mat1 +* If non-zero, then "mm" is applied first followed by "wm". Otherwise, +* "wm" is applied first followed by "mm". +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Pointer to the new MatrixMap. + +* Notes: +* - The forward direction of the returned MatrixMap is equivalent to the +* combined effect of the two supplied Mappings, operating in the +* directions specified by "winv" and "minv". +* - A null pointer will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + AstWinMap *result; /* Pointer to output WinMap */ + double *ina; /* Input corner A in new WinMap */ + double *inb; /* Input corner B in new WinMap */ + double *newscales; /* Scales for new WinMap */ + double *newshifts; /* Shifts for new WinMap */ + double *outa; /* Output corner A in new WinMap */ + double *outb; /* Output corner B in new WinMap */ + double *scales2; /* Pointer to extended WinMap scales array */ + double *scales; /* Pointer to WinMap scales array */ + double *shifts; /* Pointer to WinMap shifts array */ + int i; /* Axis index */ + int ncol; /* No. of columns in the MatrixMap */ + int nrow; /* No. of rows in the MatrixMap */ + int nt; /* Number of axes in WinMap */ + int old_minv; /* Original setting of MatrixMap Invert attribute */ + int old_winv; /* Original setting of WinMap Invert attribute */ + +/* Check the global error status. */ + if ( !astOK ) return NULL; + +/* Initialise the returned pointer. */ + result = NULL; + +/* Temporarily set the Invert attributes of both Mappings to the supplied + values. */ + old_minv = astGetInvert( mm ); + astSetInvert( mm, minv ); + + old_winv = astGetInvert( wm ); + astSetInvert( wm, winv ); + +/* Get the number of inputs (columns) and outputs (rows) for the MatrixMap. */ + ncol = astGetNin( mm ); + nrow = astGetNout( mm ); + +/* Get the scales and shifts implemented by the WinMap. These take into + account the current Invert attribute of the WinMap. */ + nt = astWinTerms( wm, &shifts, &scales ); + +/* First deal with cases where the MatrixMap is applied first. */ + if( mat1 ){ + +/* Sanity check. */ + if( nt != nrow ) { + if( astOK ) astError( AST__INTER, "astMapMerge(%s): WinMap has %d axes, " + "but MatrixMap has %d rows (internal AST programming " + "error).", status, astGetClass(mm), nt, nrow ); + + } else { + +/* Allocate the array to hold the scale terms for the new WinMap. */ + newscales = astMalloc( nrow*sizeof(double) ); + +/* Ensure that the original scales array is padded with sufficient zeros + to allow it to be transformed using the matrixmap. */ + scales2 = astCalloc( ncol, sizeof(double) ); + if( astOK ) memcpy( scales2, scales, + (ncol<nrow?ncol:nrow)*sizeof(double) ); + +/* Use the MatrixMap to transform the scale terms from the WinMap. */ + astTranN( mm, 1, ncol, 1, scales2, 1, nrow, 1, newscales ); + +/* Free resources. */ + scales2 = astFree( scales2 ); + +/* The shifts are unchanged. */ + newshifts = shifts; + } + +/* Now deal with cases where the WinMap is applied first. */ + } else { + +/* Sanity check. */ + if( nt != ncol ) { + if( astOK ) astError( AST__INTER, "astMapMerge(%s): WinMap has %d axes, " + "but MatrixMap has %d columns (internal AST programming " + "error).", status, astGetClass(mm), nt, ncol ); + + } else { + +/* Allocate the array to hold the scale and shift terms for the new WinMap. */ + newscales = astMalloc( nrow*sizeof(double) ); + newshifts = astMalloc( nrow*sizeof(double) ); + +/* Use the MatrixMap to transform the scale terms from the WinMap. */ + astTranN( mm, 1, ncol, 1, scales, 1, nrow, 1, newscales ); + +/* Use the MatrixMap to transform the shift terms from the WinMap. */ + astTranN( mm, 1, ncol, 1, shifts, 1, nrow, 1, newshifts ); + + } + } + +/* Create the new WinMap. */ + ina = astMalloc( nt*sizeof(double) ); + inb = astMalloc( nt*sizeof(double) ); + outa = astMalloc( nt*sizeof(double) ); + outb = astMalloc( nt*sizeof(double) ); + if( astOK ) { + for( i = 0; i < nt; i++ ) { + ina[ i ] = 0.0; + inb[ i ] = 1.0; + outa[ i ] = newshifts[ i ]; + outb[ i ] = newscales[ i ] + newshifts[ i ]; + } + result = astWinMap( nt, ina, inb, outa, outb, "", status ); + } + +/* Re-instate the original settings of the Invert attribute for the + supplied Mappings. */ + astSetInvert( mm, old_minv ); + astSetInvert( wm, old_winv ); + +/* Free resources. */ + ina = astFree( ina ); + inb = astFree( inb ); + outa = astFree( outa ); + outb = astFree( outb ); + if( newscales != scales ) newscales = astFree( newscales ); + if( newshifts != shifts ) newshifts = astFree( newshifts ); + scales = astFree( scales ); + shifts = astFree( shifts ); + +/* If an error has occurred, annull the returned MatrixMap. */ + if( !astOK ) result = astAnnul( result ); + +/* Return a pointer to the output MatrixMap. */ + return result; +} + +static AstMatrixMap *MatZoom( AstMatrixMap *mm, AstZoomMap *zm, int minv, + int zinv, int *status ){ +/* +* Name: +* MatZoom + +* Purpose: +* Create a MatrixMap by merging a MatrixMap and a ZoomMap. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* AstMatrixMap *MatZoom( AstMatrixMap *mm, AstZoomMap *zm, int minv, +* int zinv, int *status ) + +* Class Membership: +* MatrixMap member function + +* Description: +* This function creates a new MatrixMap which performs a mapping +* equivalent to applying the two supplied Mappings in series in the +* directions specified by the "invert" flags (the Invert attributes of +* the supplied MatrixMaps are ignored). + +* Parameters: +* mm +* A pointer to the MatrixMap. +* zm +* A pointer to the ZoomMap. +* minv +* The invert flag to use with mm. A value of zero causes the forward +* mapping to be used, and a non-zero value causes the inverse +* mapping to be used. +* zinv +* The invert flag to use with zm. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Pointer to the new MatrixMap. + +* Notes: +* - The forward direction of the returned MatrixMap is equivalent to the +* combined effect of the two supplied Mappings, operating in the +* directions specified by "zinv" and "minv". +* - A null pointer will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + AstMatrixMap *mm2; /* Pointer to intermediate MatrixMap */ + AstMatrixMap *result; /* Pointer to output MatrixMap */ + double *matrix; /* Pointer to diagonal matrix elements array */ + double zfac; /* Zoom factor */ + int i; /* Axis index */ + int nrow; /* No. of rows in the MatrixMap */ + int old_minv; /* Original setting of MatrixMap Invert attribute */ + int old_zinv; /* Original setting of ZoomMap Invert attribute */ + +/* Check the global error status. */ + if ( !astOK ) return NULL; + +/* Initialise the returned pointer. */ + result = NULL; + +/* Temporarily set the Invert attributes of both Mappings to the supplied + values. */ + old_minv = astGetInvert( mm ); + astSetInvert( mm, minv ); + + old_zinv = astGetInvert( zm ); + astSetInvert( zm, zinv ); + +/* Get the number of rows in the MatrixMap (i.e. the number of output + axes). */ + nrow = astGetNout( mm ); + +/* Get the zoom factor implemented by the ZoomMap. Invert it if necessary + since astGetZoom does not take account of the Invert setting. */ + zfac = astGetZoom( zm ); + if( zinv ) zfac = 1.0 / zfac; + +/* Create a diagonal matrix map in which each diagonal element is equal + to the zoom factor. */ + matrix = (double *) astMalloc( sizeof( double )*(size_t) nrow ); + if( astOK ) { + for( i = 0; i < nrow; i++ ) matrix[ i ] = zfac; + } + mm2 = astMatrixMap( nrow, nrow, 1, matrix, "", status ); + matrix = (double *) astFree( (void *) matrix ); + +/* Create a new MatrixMap holding the product of the supplied MatrixMap + and the diagonal MatrixMap just created. */ + result = astMtrMult( mm, mm2 ); + mm2 = astAnnul( mm2 ); + +/* Re-instate the original settings of the Invert attribute for the + supplied Mappings. */ + astSetInvert( mm, old_minv ); + astSetInvert( zm, old_zinv ); + +/* If an error has occurred, annull the returned MatrixMap. */ + if( !astOK ) result = astAnnul( result ); + +/* Return a pointer to the output MatrixMap. */ + return result; +} + +/* Functions which access class attributes. */ +/* ---------------------------------------- */ +/* Implement member functions to access the attributes associated with + this class using the macros defined for this purpose in the + "object.h" file. For a description of each attribute, see the class + interface (in the associated .h file). */ + +static AstMatrixMap *MtrMult( AstMatrixMap *this, AstMatrixMap *a, int *status ){ +/* +*+ +* Name: +* astMtrMult + +* Purpose: +* Multiply a MatrixMap by another MatrixMap. + +* Type: +* Protected virtual function. + +* Synopsis: +* #include "matrixmap.h" +* AstMatrixMap *MtrMult( astMatrixMap *this, astMatrixMap *a ) + +* Class Membership: +* MatrixMap method + +* Description: +* This function multiples the matrices given by "this" and "a", returning +* a pointer to a new MatrixMap holding the product "a x this". +* +* The number of columns in the "a" matrix must match the number of +* rows in the "this" matrix. The number of rows in the returned +* MatrixMap is equal to the number of rows in "a", and the number of +* columns is the same as the number of rows in "this". + +* Parameters: +* this +* Pointer to the first MatrixMap. +* a +* Pointer to a second MatrixMap. + +* Returned Value: +* A pointer to the product MatrixMap. + +* Notes: +* - An error is reported if the two MatrixMaps have incompatible +* shapes, or if either MatrixMap does not have a defined forward +* transformation. +* - A null Object pointer will also be returned if this function +* is invoked with the AST error status set, or if it should fail +* for any reason. +*- +*/ + +/* Local variables. */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstMatrixMap *new; /* New MatrixMap holding the product matrix */ + double *a_matrix; /* Pointer to the forward "a" matrix */ + double *a_row; /* Pointer to start of current row in "a" */ + double a_val; /* Current element value from "a" */ + double factor; /* Diagonal matrix term */ + double *new_matrix; /* Pointer to the new forward "this" matrix */ + double *new_val; /* Pointer to current output element value */ + double sum; /* Dot product value */ + double *this_col; /* Pointer to start of current column in "this" */ + double *this_matrix; /* Pointer to the forward "this" matrix */ + double this_val; /* Current element value from "this" */ + int col; /* Current output column number */ + int i; /* Loop count */ + int minrow; /* Min. number of rows in "a" or "this" */ + int ncol_a; /* No. of columns in the "a" MatrixMap */ + int ncol_this; /* No. of columns in the "this" MatrixMap */ + int nrow_a; /* No. of rows in the "a" MatrixMap */ + int nrow_this; /* No. of rows in the "this" MatrixMap */ + int row; /* Current output row number */ + +/* Return a NULL pointer if an error has already occurred. */ + if ( !astOK ) return NULL; + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(NULL); + +/* Initialise */ + new = NULL; + +/* Report an error if eitherof the MatrixMaps doe snot have a defined + forward transformation.*/ + if( !astGetTranForward( this ) ){ + astError( AST__MTRML, "astMtrMult(%s): Cannot find the product of 2 " + "MatrixMaps- the first MatrixMap has no forward transformation.", status, + astClass(this) ); + return NULL; + } + + if( !astGetTranForward( a ) ){ + astError( AST__MTRML, "astMtrMult(%s): Cannot find the product of 2 " + "MatrixMaps- the second MatrixMap has no forward transformation.", status, + astClass(this) ); + return NULL; + } + +/* Report an error if the shapes of the two matrices are incompatible. */ + nrow_a = astGetNout( a ); + ncol_a = astGetNin( a ); + nrow_this = astGetNout( this ); + ncol_this = astGetNin( this ); + + if( ncol_a != nrow_this && astOK ){ + astError( AST__MTRML, "astMtrMult(%s): Number of rows in the first " + "MatrixMap (%d) does not equal number of columns in the " + "second MatrixMap (%d).", status, astClass(this), nrow_this, ncol_a ); + return NULL; + } + +/* Store the minimum number of rows in either matrix for later use. */ + if( nrow_a < nrow_this ){ + minrow = nrow_a; + } else { + minrow = nrow_this; + } + +/* Ensure that "this" is stored in FULL form (i.e. with all elements + stored explicitly, even if the matrix is a unit or diagonal matrix). */ + ExpandMatrix( this, status ); + +/* Store pointers to the current forward matrices (taking into + account the current states of the Mapping inversion flags ). */ + this_matrix = astGetInvert( this ) ? this->i_matrix : this->f_matrix; + a_matrix = astGetInvert( a ) ? a->i_matrix : a->f_matrix; + +/* Get memory to hold the product matrix in full form. */ + new_matrix = (double *) astMalloc( sizeof( double )* + (size_t)( nrow_a*ncol_this ) ); + if( astOK ){ + +/* First deal with cases where the "a" MatrixMap represents a unit + matrix. */ + if( a->form == UNIT ){ + +/* Copy the required number of rows from "this" to "new". */ + (void) memcpy( (void *) new_matrix, (const void *) this_matrix, + sizeof(double)*(size_t)( minrow*ncol_this ) ); + +/* If there are insufficient rows in "this", append some zero-filled rows. */ + if( minrow < nrow_a ){ + for( i = minrow*ncol_this; i < nrow_a*ncol_this; i++ ){ + new_matrix[ i ] = 0.0; + } + } + +/* Now deal with cases where the "a" MatrixMap represents a diagonal + matrix. */ + } else if( a->form == DIAGONAL ){ + +/* Scale the required number of rows from "this" storing them in "new", + and checking for bad values. */ + i = 0; + + for( row = 0; row < minrow; row++ ){ + factor = a_matrix[ row ]; + + if( factor != AST__BAD ){ + + for( col = 0; col < ncol_this; col++ ){ + this_val = this_matrix[ i ]; + if( this_val != AST__BAD ){ + new_matrix[ i ] = this_val*factor; + } else { + new_matrix[ i ] = AST__BAD; + } + i++; + } + + } else { + + for( col = 0; col < ncol_this; col++ ){ + new_matrix[ i++ ] = AST__BAD; + } + + } + } + +/* If there are insufficient rows in "this", append some zero-filled rows. */ + if( minrow < nrow_a ){ + for( i = minrow*ncol_this; i < nrow_a*ncol_this; i++ ){ + new_matrix[ i ] = 0.0; + } + } + + +/* Now deal with cases where the "a" MatrixMap represents a full, non-diagonal + matrix. */ + } else { + +/* Initialise a pointer to the next element in the product matrix. */ + new_val = new_matrix; + +/* Get a pointer to the start of each row of the "a" matrix. */ + for( row = 0; row < nrow_a; row++ ){ + a_row = a_matrix + ncol_a*row; + +/* Get a pointer to the start of each column of the "this" matrix. */ + for( col = 0; col < ncol_this; col++ ){ + this_col = this_matrix + col; + +/* Form the dot product of the current row from "a", and the current + column from "this", checking for bad values. */ + sum = 0.0; + for( i = 0; i < ncol_a; i++ ){ + a_val = a_row[ i ]; + this_val = this_col[ i*ncol_this ]; + if( a_val != AST__BAD && this_val != AST__BAD ){ + sum += a_val*this_val; + } else { + sum = AST__BAD; + break; + } + } + +/* Store the output matrix element value. */ + *(new_val++) = sum; + + } + } + } + +/* Create the new MatrixMap. */ + new = astInitMatrixMap( NULL, sizeof( AstMatrixMap ), !class_init, + &class_vtab, "MatrixMap", ncol_this, nrow_a, + FULL, new_matrix ); + +/* If possible, compress the new MatrixMap by removing off-diagonal zero + elements. */ + CompressMatrix( new, status ); + +/* Re-compress the original "this" MatrixMap. */ + CompressMatrix( this, status ); + + } + +/* Free the memory used to hold the product matrix in full form. */ + new_matrix = (double *) astFree( (void *) new_matrix ); + + return new; + +} + +static AstMatrixMap *MtrRot( AstMatrixMap *this, double theta, + const double axis[], int *status ){ +/* +*+ +* Name: +* astMtrRot + +* Purpose: +* Multiply a MatrixMap by a rotation matrix. + +* Type: +* Protected virtual function. + +* Synopsis: +* #include "matrixmap.h" +* AstMatrixMap *astMtrRot( astMatrixMap *this, double theta, +* const double axis[] ) + +* Class Membership: +* MatrixMap method. + +* Description: +* This function creates a new MatrixMap which is a copy of "this", +* rotated by a specified angle. It can only be used on MatrixMaps which +* have either 2 or 3 output coordinates. In the 3-D case, the rotation +* is about an arbitrary axis passing through the origin. + +* Parameters: +* this +* Pointer to the MatrixMap. +* theta +* The angle by which to rotate the matrix, in radians. If the matrix +* is applied to a 2-D vector position, the resulting vector is +* rotated clockwise about the origin (i.e. from the positive direction +* of the second axis to the positive direction of the first axis). If +* the vector positions are three dimensional, the rotation is clockwise +* when looking along the vector given by "axis". Note, "theta" measures +f when looking along the vector given by AXIS. Note, THETA measures +* the movemement of the vectors relative to a fixed reference frame. +* Alternatively, the reference frame can be thought of as rotating by +* (-theta) relative to the fixed vectors. +* axis +* A 3-D vector specifying the axis of rotation. This parameter is +* ignored if the output from MatrixMap is 2-dimensional. + +* Returned Value: +* A pointer to the rotated MatrixMap. + +* Notes: +* - A null Object pointer will also be returned if this function +* is invoked with the AST error status set, or if it should fail +* for any reason. +*- +*/ + +/* Local variables. */ + AstMatrixMap *new; /* New MatrixMap holding the rotated matrix */ + double as,a,b,c,d,e,f,g; /* Intermediate quantities */ + double axlen; /* Length of axis vector */ + double axlen2; /* Squared length of axis vector */ + double costh; /* Cos(rotation angle) */ + double sinth; /* Sin(rotation angle) */ + double rotmat[9]; /* Rotation matrix */ + double work[3]; /* Work space for matrix multiplication */ + int ncol; /* No. of columns in the MatrixMap */ + int nrow; /* No. of rows in the MatrixMap */ + +/* Return a NULL pointer if an error has already occurred. */ + if ( !astOK ) return NULL; + +/* Initialise the returned MarixMap to be a copy of the supplied MatrixMap. */ + new = astCopy( this ); + +/* Save the cos and sin of the rotation angle for future use. */ + costh = cos( theta ); + sinth = sin( theta ); + +/* Return without changing the MatrixMap if the rotation angle is a + multiple of 360 degrees. */ + if ( costh == 1.0 ) return new; + +/* Get the dimensions of the MatrixMap. */ + nrow = astGetNout( new ); + ncol = astGetNin( new ); + +/* First do rotation of a plane about the origin. */ + if( nrow == 2 ){ + +/* Ensure that the MatrixMap is stored in full form rather than + compressed form. */ + ExpandMatrix( new, status ); + +/* Form the 2x2 forward rotation matrix. Theta is the clockwise angle + of rotation. */ + rotmat[0] = costh; + rotmat[1] = sinth; + rotmat[2] = -sinth; + rotmat[3] = costh; + +/* Post-multiply the current forward matrix (depending on whether or not + the MatrixMap has been inverted) by the forward rotation matrix. */ + if( !astGetInvert( new ) ){ + SMtrMult( 1, 2, ncol, rotmat, new->f_matrix, work, status ); + } else { + SMtrMult( 1, 2, ncol, rotmat, new->i_matrix, work, status ); + } + +/* Now form the 2x2 inverse rotation matrix (the diagonal elements + don't change). */ + rotmat[1] = -sinth; + rotmat[2] = sinth; + +/* Pre-multiply the current inverse matrix (depending on whether or + not the MatrixMap has been inverted) by the inverse rotation matrix. */ + if( !astGetInvert( new ) ){ + SMtrMult( 0, ncol, 2, rotmat, new->i_matrix, work, status ); + } else { + SMtrMult( 0, ncol, 2, rotmat, new->f_matrix, work, status ); + } + +/* See if the matrix can be stored as a UNIT or DIAGONAL matrix. */ + CompressMatrix( new, status ); + +/* Now do rotation of a volume about an axis passing through the origin. */ + } else if( nrow == 3 ){ + +/* Find the length of the axis vector. Report an error if it has zero + length or has not been supplied. */ + if( axis ) { + axlen2 = axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]; + } else { + axlen2 = 0.0; + } + if( axlen2 <= 0.0 ) { + astError( AST__MTRAX, "astMtrRot(%s): NULL or zero length " + "axis vector supplied.", status, astClass(new) ); + } + axlen = sqrt( axlen2 ); + +/* Ensure that the MatrixMap is stored in full form rather than + compressed form. */ + ExpandMatrix( new, status ); + +/* Form commonly used terms in the rotation matrix. */ + as = sinth/axlen; + a = (1.0 - costh)/axlen2; + b = a*axis[0]*axis[1]; + c = as*axis[2]; + d = a*axis[0]*axis[2]; + e = as*axis[1]; + f = a*axis[1]*axis[2]; + g = as*axis[0]; + +/* Form the 3x3 forward rotation matrix. Theta is the clockwise angle + of rotation looking in the direction of the axis vector. */ + rotmat[0] = a*axis[0]*axis[0] + costh; + rotmat[1] = b - c; + rotmat[2] = d + e; + rotmat[3] = b + c; + rotmat[4] = a*axis[1]*axis[1] + costh; + rotmat[5] = f - g; + rotmat[6] = d - e; + rotmat[7] = f + g; + rotmat[8] = a*axis[2]*axis[2] + costh; + +/* Post-multiply the current forward matrix (depending on whether or not + the MatrixMap has been inverted) by the forward rotation matrix. */ + if( !astGetInvert( new ) ){ + SMtrMult( 1, 3, ncol, rotmat, new->f_matrix, work, status ); + } else { + SMtrMult( 1, 3, ncol, rotmat, new->i_matrix, work, status ); + } + +/* Now form the 3x3 inverse rotation matrix (the diagonal elements + don't change). */ + rotmat[1] = b + c; + rotmat[2] = d - e; + rotmat[3] = b - c; + rotmat[5] = f + g; + rotmat[6] = d + e; + rotmat[7] = f - g; + +/* Pre-multiply the current inverse matrix (depending on whether or + not the MatrixMap has been inverted) by the inverse rotation matrix. */ + if( !astGetInvert( new ) ){ + SMtrMult( 0, ncol, 3, rotmat, new->i_matrix, work, status ); + } else { + SMtrMult( 0, ncol, 3, rotmat, new->f_matrix, work, status ); + } + +/* See if the matrix can be stored as a UNIT or DIAGONAL matrix. */ + CompressMatrix( new, status ); + +/* Report an error if the matrix is not suitable for rotation. */ + } else { + astError( AST__MTR23, "astMtrRot(%s): Cannot rotate a %dx%d" + " MatrixMap.", status, astClass(new), nrow, ncol ); + } + +/* Delete the new MatrixMap if an error has occurred. */ + if( !astOK ) new = astDelete( new ); + + return new; + +} + +static void PermGet( AstPermMap *map, int **outperm, int **inperm, + double **consts, int *status ){ +/* +* Name: +* PermGet + +* Purpose: +* Get the axis permutation and constants array for a PermMap. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* void PermGet( AstPermMap *map, int **outperm, int **inperm, +* double **const, int *status ) + +* Class Membership: +* MatrixMap member function + +* Description: +* This function returns axis permutation and constants arrays which can +* be used to create a PermMap which is equivalent to the supplied PermMap. + +* Parameters: +* map +* The PermMap. +* outperm +* An address at which to return a popinter to an array of ints +* holding the output axis permutation array. The array should be +* released using astFree when no longer needed. +* inperm +* An address at which to return a popinter to an array of ints +* holding the input axis permutation array. The array should be +* released using astFree when no longer needed. +* consts +* An address at which to return a popinter to an array of doubles +* holding the constants array. The array should be released using +* astFree when no longer needed. +* status +* Pointer to the inherited status variable. + +* Notes: +* - NULL pointers are returned if an error has already occurred, or if +* this function should fail for any reason. +*/ + +/* Local Variables: */ + AstPointSet *pset1; /* PointSet holding input positions for PermMap */ + AstPointSet *pset2; /* PointSet holding output positions for PermMap */ + double **ptr1; /* Pointer to pset1 data */ + double **ptr2; /* Pointer to pset2 data */ + double *cnst; /* Pointer to constants array */ + double cn; /* Potential new constant value */ + double ip; /* Potential output axis index */ + double op; /* Potential input axis index */ + int *inprm; /* Pointer to input axis permutation array */ + int *outprm; /* Pointer to output axis permutation array */ + int i; /* Axis count */ + int nc; /* Number of constants stored so far */ + int nin; /* No. of input coordinates for the PermMap */ + int nout; /* No. of output coordinates for the PermMap */ + +/* Initialise. */ + if( outperm ) *outperm = NULL; + if( inperm ) *inperm = NULL; + if( consts ) *consts = NULL; + +/* Check the global error status and the supplied pointers. */ + if ( !astOK || !outperm || !inperm || !consts ) return; + +/* Get the number of input and output axes for the supplied PermMap. */ + nin = astGetNin( map ); + nout = astGetNout( map ); + +/* Allocate the memory for the returned arrays. */ + outprm = (int *) astMalloc( sizeof( int )* (size_t) nout ); + inprm = (int *) astMalloc( sizeof( int )* (size_t) nin ); + cnst = (double *) astMalloc( sizeof( double )* (size_t) ( nout + nin ) ); + +/* Returned the pointers to these arrays.*/ + *outperm = outprm; + *inperm = inprm; + *consts = cnst; + +/* Create two PointSets, each holding two points, which can be used for + input and output positions with the PermMap. */ + pset1 = astPointSet( 2, nin, "", status ); + pset2 = astPointSet( 2, nout, "", status ); + +/* Set up the two input positions to be [0,1,2...] and [-1,-1,-1,...]. The + first position is used to enumerate the axes, and the second is used to + check for constant axis values. */ + ptr1 = astGetPoints( pset1 ); + if( astOK ){ + for( i = 0; i < nin; i++ ){ + ptr1[ i ][ 0 ] = ( double ) i; + ptr1[ i ][ 1 ] = -1.0; + } + } + +/* Use the PermMap to transform these positions in the forward direction. */ + (void) astTransform( map, pset1, 1, pset2 ); + +/* No constant axis valeus found yet. */ + nc = 0; + +/* Look at the mapped positions to determine the output axis permutation + array. */ + ptr2 = astGetPoints( pset2 ); + if( astOK ){ + +/* Do each output axis. */ + for( i = 0; i < nout; i++ ){ + +/* If the output axis value is copied from an input axis value, the index + of the appropriate input axis will be in the mapped first position. */ + op = ptr2[ i ][ 0 ]; + +/* If the output axis value is assigned a constant value, the result of + mapping the two different input axis values will be the same. */ + cn = ptr2[ i ][ 1 ]; + if( op == cn ) { + +/* We have found another constant. Store it in the constants array, and + store the index of the constant in the output axis permutation array. */ + cnst[ nc ] = cn; + outprm[ i ] = -( nc + 1 ); + nc++; + +/* If the output axis values are different, then the output axis value + must be copied from the input axis value. */ + } else { + outprm[ i ] = (int) ( op + 0.5 ); + } + } + } + +/* Now do the same thing to determine the input permutation array. */ + if( astOK ){ + for( i = 0; i < nout; i++ ){ + ptr2[ i ][ 0 ] = ( double ) i; + ptr2[ i ][ 1 ] = -1.0; + } + } + + (void) astTransform( map, pset2, 0, pset1 ); + + if( astOK ){ + + for( i = 0; i < nin; i++ ){ + + ip = ptr1[ i ][ 0 ]; + cn = ptr1[ i ][ 1 ]; + if( ip == cn ) { + + cnst[ nc ] = cn; + inprm[ i ] = -( nc + 1 ); + nc++; + + } else { + inprm[ i ] = (int) ( ip + 0.5 ); + } + } + } + +/* Annul the PointSets. */ + pset1 = astAnnul( pset1 ); + pset2 = astAnnul( pset2 ); + +/* If an error has occurred, attempt to free the returned arrays. */ + if( !astOK ) { + *outperm = (int *) astFree( (void *) *outperm ); + *inperm = (int *) astFree( (void *) *inperm ); + *consts = (double *) astFree( (void *) *consts ); + } + +/* Return. */ + return; +} + +static int PermOK( AstMapping *pm, int *status ){ +/* +* Name: +* PermOK + +* Purpose: +* Determine if a PermMap can be merged with a MatrixMap. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* int PermOK( AstMapping *pm, int *status ) + +* Class Membership: +* PermMap member function + +* Description: +* This function returns a flag indicating if the supplied PermMap +* could be merged with a MatrixMap. For thios to be possible, the +* PermMap must have the same number of input and output axes, and the +* inverse and forward mappings must be consistent. + +* Parameters: +* pm +* The PermMap. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* 1 if the PermMap can be merged, 0 otherwise. + +* Notes: +* - A value of 0 is returned if an error has already occurred, or if +* this function should fail for any reason. +*/ + +/* Local Variables: */ + AstPointSet *pset1; /* PointSet holding input positions for PermMap */ + AstPointSet *pset2; /* PointSet holding output positions for PermMap */ + double **ptr1; /* Pointer to pset1 data */ + int i; /* Loop count */ + int nin; /* No. of input coordinates for the PermMap */ + int nout; /* No. of output coordinates for the PermMap */ + int ret; /* Returned flag */ + +/* Check the global error status. */ + if ( !astOK ) return 0; + +/* Initialise */ + ret = 0; + +/* The PermMap must have the same number of input and output coordinates. */ + nin = astGetNin( pm ); + nout = astGetNout( pm ); + if( nin == nout ){ + +/* Create two PointSets, each holding two points, which can be used for + the input and output positions with the PermMap. */ + pset1 = astPointSet( 2, nin, "", status ); + pset2 = astPointSet( 2, nout, "", status ); + +/* Set up the two input positions to be [1,2,3...] and [0,-1,-2,...] */ + ptr1 = astGetPoints( pset1 ); + if( astOK ){ + for( i = 0; i < nin; i++ ){ + ptr1[ i ][ 0 ] = ( double )( i + 1 ); + ptr1[ i ][ 1 ] = ( double )( -i ); + } + + } + +/* Use the PermMap to transform these positions in the forward direction. */ + (void) astTransform( pm, pset1, 1, pset2 ); + +/* Now transform the results back again using the inverse PermMap. */ + (void) astTransform( pm, pset2, 0, pset1 ); + +/* See if the input positions have changed. If they have, then the PermMap + does not have a consistent pair of transformations. If they have not, + then the transformations must be consistent because we used two + different input positions and only one could come out unchanged by + chance. */ + if( astOK ){ + ret = 1; + for( i = 0; i < nin; i++ ){ + if( ptr1[ i ][ 0 ] != ( double )( i + 1 ) || + ptr1[ i ][ 1 ] != ( double )( -i ) ){ + ret = 0; + break; + } + } + } + +/* Annul the PointSets. */ + pset1 = astAnnul( pset1 ); + pset2 = astAnnul( pset2 ); + } + +/* Return the answer. */ + return astOK ? ret : 0; +} + +static double Rate( AstMapping *this, double *at, int ax1, int ax2, int *status ){ +/* +* Name: +* Rate + +* Purpose: +* Calculate the rate of change of a Mapping output. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* result = Rate( AstMapping *this, double *at, int ax1, int ax2, int *status ) + +* Class Membership: +* MatrixMap member function (overrides the astRate method inherited +* from the Mapping class ). + +* Description: +* This function returns the rate of change of a specified output of +* the supplied Mapping with respect to a specified input, at a +* specified input position. + +* Parameters: +* this +* Pointer to the Mapping to be applied. +* at +* The address of an array holding the axis values at the position +* at which the rate of change is to be evaluated. The number of +* elements in this array should equal the number of inputs to the +* Mapping. +* ax1 +* The index of the Mapping output for which the rate of change is to +* be found (output numbering starts at 0 for the first output). +* ax2 +* The index of the Mapping input which is to be varied in order to +* find the rate of change (input numbering starts at 0 for the first +* input). +* status +* Pointer to the inherited status variable. + +* Returned Value: +* The rate of change of Mapping output "ax1" with respect to input +* "ax2", evaluated at "at", or AST__BAD if the value cannot be +* calculated. + +*/ + +/* Local Variables: */ + AstMatrixMap *map; + double *matrix; + double result; + +/* Check inherited status */ + if( !astOK ) return AST__BAD; + +/* Get a pointer to the MatrixMap structure. */ + map = (AstMatrixMap *) this; + +/* Get a pointer to the array holding the required matrix elements, according + to whether the MatrixMap has been inverted. */ + if( !astGetInvert( this ) ) { + matrix = map->f_matrix; + } else { + matrix = map->i_matrix; + } + +/* First deal with full MatrixMaps in which all matrix elements are stored. */ + if( map->form == FULL ){ + result = matrix[ ax1*astGetNin( this ) + ax2 ]; + +/* For unit matrices, the rate is unity if the input and output axes are + equal, and zero otherwise. */ + } else if( map->form == UNIT ){ + result = (ax1 == ax2 ) ? 1.0 : 0.0; + +/* For diagonal matrices, the rate is zero for off diagonal elements and + the matrix array stored the on-diagonal rates. */ + } else if( ax1 == ax2 ) { + result = matrix[ ax1 ]; + + } else { + result = 0.0; + } + +/* Return the result. */ + return result; +} + +static void SMtrMult( int post, int m, int n, const double *mat1, + double *mat2, double *work, int *status ){ +/* +* Name: +* SMtrMult + +* Purpose: +* Multiply a square matrix and a non-square matrix. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* void SMtrMult( int post, int m, int n, const double *mat1, +* double *mat2, double *work, int *status ) + +* Class Membership: +* MatrixMap member function. + +* Description: +* The matrix pointed to by "mat2" is modified by multiplying it by +* the square matrix pointed to by "mat1". If "post" is 1, then: +* +* mat2 -> mat1*mat2 (mat1 is mxm and mat2 is mxn) +* +* If "post" is zero, then: +* +* mat2 -> mat2*mat1 (mat1 is nxn and mat2 is mxn) +* +* The restriction that "mat1" must be square is imposed so that the +* returned matrix will have the same shape as the supplied matrix (mat1). + +* Parameters: +* post +* Specifies whether to post- or pre- multiply mat2 by mat1. +* m +* The number of rows in mat2. +* n +* The number of columns in mat2. +* mat1 +* The multiplier matrix. It must be square of size m or n, depending +* on "post". +* mat2 +* The multiplicand matrix. +* work +* Pointer to work space containing room for m doubles (if "post" +* is 1), or n doubles (if "post" is 0). +* status +* Pointer to the inherited status variable. + +* Notes: +* - No error is reported if "mat2" is supplied NULL. In this case +* it will also be returned NULL. +*/ + +/* Local Variables */ + double dot; /* Output matrix element value */ + const double *mat1_col; /* Pointer to start of current column of mat1 */ + const double *mat1_row; /* Pointer to start of current row of mat1 */ + double *mat2_col; /* Pointer to start of current column of mat2 */ + double *mat2_row; /* Pointer to start of current row of mat2 */ + double cel; /* Column element value */ + double rel; /* Row element value */ + int i; /* Index of current output row */ + int j; /* Index of current output column */ + int k; /* Dot product index */ + +/* Do nothing if mat2 is NULL */ + if ( mat2 ){ + +/* First deal with cases where the supplied matrix is post-multiplied + (i.e. the returned matrix is mat1*mat2). */ + if( post ){ + +/* Loop round each column of the output matrix, storing a pointer to + the start of the corresponding column of mat2. */ + for( j=0; j<n; j++ ){ + mat2_col = mat2 + j; + +/* Loop round each row of the output matrix, storing a pointer to + the start of the corresponding row of mat1. */ + for( i=0; i<m; i++ ){ + mat1_row = mat1 + i*m; + +/* Get the dot product of the corresponding row from mat1 and the + corresponding column from mat2 and store it in the work array. */ + dot = 0.0; + for( k=0; k<m; k++ ) { + rel = mat1_row[ k ]; + cel = mat2_col[ k*n ]; + if( rel != AST__BAD && cel != AST__BAD ){ + dot += rel*cel; + } else { + dot = AST__BAD; + break; + } + } + work[ i ] = dot; + } + +/* Copy the values stored in the work array to the current column of + the output matrix. */ + for( i=0; i<m; i++ ) mat2_col[ i*n ] = work[ i ]; + } + +/* Now deal with cases where the supplied matrix is pre-multiplied + (i.e. the returned matrix is mat2*mat1). */ + } else { + +/* Loop round each row of the output matrix, storing a pointer to + the start of the corresponding row of mat2. */ + for( i=0; i<m; i++ ){ + mat2_row = mat2 + i*n; + +/* Loop round each column of the output matrix, storing a pointer to + the start of the corresponding column of mat1. */ + for( j=0; j<n; j++ ){ + mat1_col = mat1 + j; + +/* Get the dot product of the corresponding row from mat2 and the + corresponding column from mat1 and store it in the work array. */ + dot = 0.0; + for( k=0; k<n; k++ ) { + rel = mat2_row[ k ]; + cel = mat1_col[ k*n ]; + if( rel != AST__BAD && cel != AST__BAD ){ + dot += rel*cel; + } else { + dot = AST__BAD; + break; + } + } + work[ j ] = dot; + } + +/* Copy the values stored in the work array to the current row of + the output matrix. */ + for( j=0; j<n; j++ ) mat2_row[ j ] = work[ j ]; + } + } + } + + return; + +} + +static int GetTranForward( AstMapping *this, int *status ) { +/* +* +* Name: +* GetTranForward + +* Purpose: +* Determine if a MatrixMap defines a forward coordinate transformation. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* int GetTranForward( AstMapping *this, int *status ) + +* Class Membership: +* MatrixMap member function (over-rides the astGetTranForward method +* inherited from the Mapping class). + +* Description: +* This function returns a value indicating if the MatrixMap is able +* to perform a forward coordinate transformation. + +* Parameters: +* this +* Pointer to the MatrixMap. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Zero if the forward coordinate transformation is not defined, or 1 if it +* is. + +* Notes: +* - A value of zero will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + AstMatrixMap *map; /* Pointer to MatrixMap to be queried */ + int invert; /* Has the mapping been inverted? */ + int result; /* The returned value */ + +/* Initialise. */ + result = 0; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Obtain a pointer to the MatrixMap. */ + map = (AstMatrixMap *) this; + +/* All unit MatrixMaps are defined in both directions. */ + if( map->form == UNIT ) { + result = 1; + +/* Otherwise, check that the appropriate array is defined in the + MatrixMap. */ + } else { + +/* Determine if the Mapping has been inverted. */ + invert = astGetInvert( this ); + +/* If OK, obtain the result. */ + if ( astOK ) { + + if( invert ){ + result = ( map->i_matrix != NULL ); + } else { + result = ( map->f_matrix != NULL ); + } + + } + + } + +/* Return the result. */ + return result; + +} + +static int GetTranInverse( AstMapping *this, int *status ) { +/* +* +* Name: +* GetTranInverse + +* Purpose: +* Determine if a MatrixMap defines an inverse coordinate transformation. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* int GetTranInverse( AstMapping *this, int *status ) + +* Class Membership: +* MatrixMap member function (over-rides the astGetTranInverse method +* inherited from the Mapping class). + +* Description: +* This function returns a value indicating if the MatrixMap is able +* to perform an inverse coordinate transformation. + +* Parameters: +* this +* Pointer to the MatrixMap. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Zero if the inverse coordinate transformation is not defined, or 1 if it +* is. + +* Notes: +* - A value of zero will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +*/ + +/* Local Variables: */ + AstMatrixMap *map; /* Pointer to MatrixMap to be queried */ + int invert; /* Has the mapping been inverted? */ + int result; /* The returned value */ + +/* Initialise. */ + result = 0; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Obtain a pointer to the MatrixMap. */ + map = (AstMatrixMap *) this; + +/* All unit MatrixMaps are defined in both directions. */ + if( map->form == UNIT ) { + result = 1; + +/* Otherwise, check that the appropriate array is defined in the + MatrixMap. */ + } else { + +/* Determine if the Mapping has been inverted. */ + invert = astGetInvert( this ); + +/* If OK, obtain the result. */ + if ( astOK ) { + + if( invert ){ + result = ( map->f_matrix != NULL ); + } else { + result = ( map->i_matrix != NULL ); + } + + } + + } + +/* Return the result. */ + return result; + +} + +static AstPointSet *Transform( AstMapping *this, AstPointSet *in, + int forward, AstPointSet *out, int *status ) { +/* +* Name: +* Transform + +* Purpose: +* Apply a MatrixMap to transform a set of points. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* AstPointSet *Transform( AstMapping *this, AstPointSet *in, +* int forward, AstPointSet *out, int *status ) + +* Class Membership: +* MatrixMap member function (over-rides the astTransform protected +* method inherited from the Mapping class). + +* Description: +* This function takes a MatrixMap and a set of points encapsulated in a +* PointSet and transforms the points by multiplying them by the matrix. + +* Parameters: +* this +* Pointer to the MatrixMap. +* in +* Pointer to the PointSet holding the input coordinate data. +* forward +* A non-zero value indicates that the forward coordinate transformation +* should be applied, while a zero value requests the inverse +* transformation. +* out +* Pointer to a PointSet which will hold the transformed (output) +* coordinate values. A NULL value may also be given, in which case a +* new PointSet will be created by this function. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Pointer to the output (possibly new) PointSet. + +* Notes: +* - A null pointer will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +* - The number of coordinate values per point in the input PointSet must +* match the number of columns in the MatrixMap being applied. +* - The number of coordinate values per point in the output PointSet will +* equal the number of rows in the MatrixMap being applied. +* - If an output PointSet is supplied, it must have space for sufficient +* number of points and coordinate values per point to accommodate the +* result. Any excess space will be ignored. +*/ + +/* Local Variables: */ + AstPointSet *result; /* Pointer to output PointSet */ + AstMatrixMap *map; /* Pointer to MatrixMap to be applied */ + double diag_term; /* Current diagonal element value */ + double *indata; /* Pointer to next input data value */ + double *matrix; /* Pointer to start of matrix element array */ + double *matrix_element; /* Pointer to current matrix element value */ + double *outdata; /* Pointer to next output data value */ + double **ptr_in; /* Pointer to input coordinate data */ + double **ptr_out; /* Pointer to output coordinate data */ + double sum; /* Partial output value */ + double val; /* Data value */ + int in_coord; /* Index of output coordinate */ + int nax; /* Output axes for which input axes exist */ + int ncoord_in; /* Number of coordinates per input point */ + int ncoord_out; /* Number of coordinates per output point */ + int npoint; /* Number of points */ + int out_coord; /* Index of output coordinate */ + int point; /* Loop counter for points */ + +/* Check the global error status. */ + if ( !astOK ) return NULL; + +/* Obtain a pointer to the MatrixMap. */ + map = (AstMatrixMap *) this; + +/* Apply the parent mapping using the stored pointer to the Transform member + function inherited from the parent Mapping class. This function validates + all arguments and generates an output PointSet if necessary, but does not + actually transform any coordinate values. */ + result = (*parent_transform)( this, in, forward, out, status ); + +/* We will now extend the parent astTransform method by performing the + calculations needed to generate the output coordinate values. */ + +/* Determine the numbers of points and coordinates per point from the input + and output PointSets and obtain pointers for accessing the input and + output coordinate values. */ + ncoord_in = astGetNcoord( in ); + ncoord_out = astGetNcoord( result ); + npoint = astGetNpoint( in ); + ptr_in = astGetPoints( in ); + ptr_out = astGetPoints( result ); + +/* Determine whether to apply the forward or inverse mapping, according to the + direction specified and whether the mapping has been inverted. */ + if ( astGetInvert( map ) ) forward = !forward; + +/* Get a pointer to the array holding the required matrix elements, according + to the direction of mapping required. */ + if ( forward ) { + matrix = map->f_matrix; + } else { + matrix = map->i_matrix; + } + +/* Perform coordinate arithmetic. */ +/* ------------------------------ */ + if ( astOK ) { + +/* First deal with full MatrixMaps in which all matrix elements are stored. */ + if( map->form == FULL ){ + +/* Loop to apply the matrix to each point in turn, checking for + (and propagating) bad values in the process. The matrix elements are + accessed sequentially in row order. The next matrix element to be + used is identified by a pointer which is initialised to point to the + first element of the matrix prior to processing each point. */ + for ( point = 0; point < npoint; point++ ) { + matrix_element = matrix; + +/* Each output co-ordinate value is created by summing the product of the + corresponding input co-ordinates and the elements of one row of the + matrix. */ + for ( out_coord = 0; out_coord < ncoord_out; out_coord++ ) { + sum = 0.0; + + for ( in_coord = 0; in_coord < ncoord_in; in_coord++ ) { + +/* Check the current input coordinate value and the current matrix element. + If the coordinate value is bad, then the output value will also be + bad unless the matrix element is zero. That is, a zero matrix element + results in the input coordinate value being ignored, even if it is bad. + This prevents bad input values being propagated to output axes which + are independant of the bad input axis. A bad matrix element always results + in the output value being bad. In either of these cases, break out of the + loop, remembering to advance the pointer to the next matrix element so + that it points to the start of the next row ready for doing the next + output coordinate. */ + if ( ( ptr_in[ in_coord ][ point ] == AST__BAD && + (*matrix_element) != 0.0 ) || + (*matrix_element) == AST__BAD ) { + sum = AST__BAD; + matrix_element += ncoord_in - in_coord; + break; + +/* If the input coordinate and the current matrix element are both + valid, increment the sum by their product, and step to the next matrix + element pointer If we arrive here with a bad input value, then the + matrix element must be zero, in which case the running sum is left + unchanged. */ + } else { + if ( ptr_in[ in_coord ][ point ] != AST__BAD ) { + sum += ptr_in[ in_coord ][ point ] * (*matrix_element); + } + matrix_element++; + } + } + +/* Store the output coordinate value. */ + ptr_out[ out_coord ][ point ] = sum; + + } + + } + +/* Now deal with unit and diagonal MatrixMaps. */ + } else { + +/* Find the number of output axes for which input data is available. */ + if( ncoord_in < ncoord_out ){ + nax = ncoord_in; + } else { + nax = ncoord_out; + } + +/* For unit matrices, copy the input axes to the corresponding output axes. */ + if( map->form == UNIT ){ + for( out_coord = 0; out_coord < nax; out_coord++ ) { + (void) memcpy( ptr_out[ out_coord ], + (const void *) ptr_in[ out_coord ], + sizeof( double )*(size_t)npoint ); + } + +/* For diagonal matrices, scale each input axis using the appropriate + diagonal element from the matrix, and store in the output. */ + } else { + for( out_coord = 0; out_coord < nax; out_coord++ ){ + diag_term = matrix[ out_coord ]; + outdata = ptr_out[ out_coord ]; + indata = ptr_in[ out_coord ]; + + if( diag_term != AST__BAD ){ + for( point = 0; point < npoint; point++ ){ + val = *(indata++); + if( val != AST__BAD ){ + *(outdata++) = diag_term*val; + } else { + *(outdata++) = AST__BAD; + } + } + + } else { + for( point = 0; point < npoint; point++ ){ + *(outdata++) = AST__BAD; + } + } + } + } + +/* If there are any remaining output axes, fill the first one with zeros. */ + if( nax < ncoord_out ){ + outdata = ptr_out[ nax ]; + for( point = 0; point < npoint; point++ ) *(outdata++) = 0.0; + +/* Copy this axis to any remaining output axes. */ + outdata = ptr_out[ nax ]; + for( out_coord = nax + 1; out_coord < ncoord_out; out_coord++ ) { + (void) memcpy( ptr_out[ out_coord ], (const void *) outdata, + sizeof( double )*(size_t)npoint ); + } + } + } + } + +/* Return a pointer to the output PointSet. */ + return result; +} + +static int ScalingRowCol( AstMatrixMap *map, int axis, int *status ){ +/* +* Name: +* ScalingRowCol + +* Purpose: +* Determine if a given row and column of a MatrixMap are zeros +* with a non-zero diagonal term. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* int ScalingRowCol( AstMatrixMap *map, int axis, int *status ) + +* Class Membership: +* MatrixMap member function + +* Description: +* This function returns a flag indicating if a MatrixMap presents a +* simple scaling for a given axis in both directions. The MatrixMap +* must be square. A value of one is returned if every element of the +* row and column corresponding to the given axis is zero, except for +* the diagonal term which must be non-zero. + +* Parameters: +* map +* The MatrixMap. +* axis +* The zero-based index of the axis to check. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* 1 if the row/column produces a simple scaling, 0 otherwise. + +*/ + +/* Local Variables: */ + double *el; /* Pointer to matrix element */ + int i; /* Element count */ + int ncol; /* No. of input coordinates */ + int ret; /* Returned flag */ + +/* Initialise */ + ret = 0; + +/* Check the global error status. */ + if ( !astOK ) return ret; + +/* If a unit or diagonal MatrixMap has been supplied, return 1. */ + if( map->form != FULL ){ + ret = 1; + +/* If a full matrix has been supplied... */ + } else { + +/* Assume the row/column gives a unit mapping. */ + ret = 1; + +/* Get the number of input axes for the MatrixMap. */ + ncol = astGetNin( map ); + +/* Check that all elements of the "axis"th row are effectively zero, except + for the "axis"th element which must be non-zero. */ + el = map->f_matrix + axis*ncol; + for( i = 0; i < ncol; i++ ) { + if( i == axis ) { + if( fabs( *el ) <= DBL_EPSILON ) { + ret = 0; + break; + } + } else if( fabs( *el ) > DBL_EPSILON ) { + ret = 0; + break; + } + el++; + } + +/* Check that all elements of the "axis"th column are effectively zero, except + for the "axis"th element which must be non-zero. */ + if( ret ) { + el = map->f_matrix + axis; + for( i = 0; i < ncol; i++ ) { + if( i == axis ) { + if( fabs( *el ) <= DBL_EPSILON ) { + ret = 0; + break; + } + } else if( fabs( *el ) > DBL_EPSILON ) { + ret = 0; + break; + } + el += ncol; + } + } + } + +/* Return the answer. */ + return astOK ? ret : 0; +} + +/* Functions which access class attributes. */ +/* ---------------------------------------- */ +/* Implement member functions to access the attributes associated with + this class using the macros defined for this purpose in the + "object.h" file. For a description of each attribute, see the class + interface (in the associated .h file). */ + +/* Copy constructor. */ +/* ----------------- */ +static void Copy( const AstObject *objin, AstObject *objout, int *status ) { +/* +* Name: +* Copy + +* Purpose: +* Copy constructor for MatrixMap objects. + +* Type: +* Private function. + +* Synopsis: +* void Copy( const AstObject *objin, AstObject *objout, int *status ) + +* Description: +* This function implements the copy constructor for MatrixMap objects. + +* Parameters: +* objin +* Pointer to the object to be copied. +* objout +* Pointer to the object being constructed. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* void + +* Notes: +* - This constructor makes a deep copy, including a copy of the matrix +* element values associated with the input MatrixMap. +*/ + + +/* Local Variables: */ + AstMatrixMap *in; /* Pointer to input MatrixMap */ + AstMatrixMap *out; /* Pointer to output MatrixMap */ + int nel; /* No. of elements in the matrix */ + int nin; /* No. of input coordinates */ + int nout; /* No. of output coordinates */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain pointers to the input and output MatrixMaps. */ + in = (AstMatrixMap *) objin; + out = (AstMatrixMap *) objout; + +/* Nullify the pointers stored in the output object since these will + currently be pointing at the input data (since the output is a simple + byte-for-byte copy of the input). Otherwise, the input data could be + freed by accidient if the output object is deleted due to an error + occuring in this function. */ + out->f_matrix = NULL; + out->i_matrix = NULL; + +/* If the input MatrixMap is a unit mapping, then no matrix elements are + stored with it, so do nothing in this case. */ + if( out->form != UNIT ){ + +/* Obtain the number of stored values in the MatrixMap. This is independant of + whether the Mapping has been inverted or not. If the MatrixMap is diagonal, + only the diagonal terms are stored. */ + nin = astGetNin( in ); + nout = astGetNout( in ); + + if( out->form == DIAGONAL ){ + if( nin < nout ){ + nel = nin; + } else { + nel = nout; + } + + } else { + nel = nin*nout; + } + +/* Store the forward matrix elements in the output MatrixMap. */ + out->f_matrix = (double *) astStore( NULL, (void *) in->f_matrix, + sizeof( double )*(size_t) nel ); + +/* Store the inverse matrix elements (if defined) in the output + MatrixMap. */ + if( in->i_matrix ){ + out->i_matrix = (double *) astStore( NULL, (void *) in->i_matrix, + sizeof( double )*(size_t) nel ); + } + +/* If an error has occurred, free the output MatrixMap arrays. */ + if( !astOK ) { + out->f_matrix = (double *) astFree( (void *) out->f_matrix ); + out->i_matrix = (double *) astFree( (void *) out->i_matrix ); + } + } + + return; + +} + +/* Destructor. */ +/* ----------- */ +static void Delete( AstObject *obj, int *status ) { +/* +* Name: +* Delete + +* Purpose: +* Destructor for MatrixMap objects. + +* Type: +* Private function. + +* Synopsis: +* void Delete( AstObject *obj, int *status ) + +* Description: +* This function implements the destructor for MatrixMap objects. + +* Parameters: +* obj +* Pointer to the object to be deleted. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* void + +* Notes: +* This function attempts to execute even if the global error status is +* set. +*/ + +/* Local Variables: */ + AstMatrixMap *this; /* Pointer to MatrixMap */ + +/* Obtain a pointer to the MatrixMap structure. */ + this = (AstMatrixMap *) obj; + +/* Free the arrays used to store element values for forward and inverse + matrices. */ + this->f_matrix = (double *) astFree( (void *) this->f_matrix ); + this->i_matrix = (double *) astFree( (void *) this->i_matrix ); +} + +/* Dump function. */ +/* -------------- */ +static void Dump( AstObject *this_object, AstChannel *channel, int *status ) { +/* +* Name: +* Dump + +* Purpose: +* Dump function for MatrixMap objects. + +* Type: +* Private function. + +* Synopsis: +* void Dump( AstObject *this, AstChannel *channel, int *status ) + +* Description: +* This function implements the Dump function which writes out data +* for the MatrixMap class to an output Channel. + +* Parameters: +* this +* Pointer to the MatrixMap whose data are being written. +* channel +* Pointer to the Channel to which the data are being written. +* status +* Pointer to the inherited status variable. +*/ + +#define KEY_LEN 50 /* Maximum length of a keyword */ + +/* Local Variables: */ + AstMatrixMap *this; /* Pointer to the MatrixMap structure */ + char buff[ KEY_LEN + 1 ]; /* Buffer for keyword string */ + int el; /* Element index */ + int nel; /* No. of elements in the matrix */ + int nin; /* No. of input coords */ + int nout; /* No. of output coords */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain a pointer to the MatrixMap structure. */ + this = (AstMatrixMap *) this_object; + +/* Find the number of elements stored for each matrix. */ + nin = astGetNin( this ); + nout = astGetNout( this ); + + if( this->form == FULL ){ + nel = nin*nout; + + } else if( this->form == DIAGONAL ){ + nel = astMIN( nin, nout ); + + } else { + nel = 0; + } + +/* Write out values representing the instance variables for the + MatrixMap class. */ + +/* The forward matrix. Note BAD values are not written out as the + AST__BAD value may differ on different machines. If a matrix element + is not found when reading the matrix back in again (in astLoadMatrixMap), + then it is assigned a default value of AST__BAD. */ + if( this->f_matrix ){ + for( el = 0; el < nel; el++ ){ + if( (this->f_matrix)[ el ] != AST__BAD ) { + (void) sprintf( buff, "M%d", el ); + astWriteDouble( channel, buff, 1, 1, (this->f_matrix)[ el ], + "Forward matrix value" ); + } + } + } + +/* The inverse matrix. */ + if( this->i_matrix ){ + for( el = 0; el < nel; el++ ){ + if( (this->i_matrix)[ el ] != AST__BAD ) { + (void) sprintf( buff, "IM%d", el ); + astWriteDouble( channel, buff, 1, 1, (this->i_matrix)[ el ], + "Inverse matrix value" ); + } + } + } + +/* The matrix storage form. */ + astWriteString( channel, "Form", 1, 1, Form[ this->form ], + "Matrix storage form" ); + +/* Undefine macros local to this function. */ +#undef KEY_LEN +} + +/* Standard class functions. */ +/* ========================= */ +/* Implement the astIsAMatrixMap and astCheckMatrixMap functions using the macros + defined for this purpose in the "object.h" header file. */ +astMAKE_ISA(MatrixMap,Mapping) +astMAKE_CHECK(MatrixMap) + +AstMatrixMap *astMatrixMap_( int nin, int nout, int form, + const double matrix[], const char *options, int *status, ...){ +/* +*++ +* Name: +c astMatrixMap +f AST_MATRIXMAP + +* Purpose: +* Create a MatrixMap. + +* Type: +* Public function. + +* Synopsis: +c #include "matrixmap.h" +c AstMatrixMap *astMatrixMap( int nin, int nout, int form, +c const double matrix[], +c const char *options, ... ) +f RESULT = AST_MATRIXMAP( NIN, NOUT, FORM, MATRIX, OPTIONS, STATUS ) + +* Class Membership: +* MatrixMap constructor. + +* Description: +* This function creates a new MatrixMap and optionally initialises +* its attributes. +* +* A MatrixMap is a form of Mapping which performs a general linear +* transformation. Each set of input coordinates, regarded as a +* column-vector, are pre-multiplied by a matrix (whose elements +* are specified when the MatrixMap is created) to give a new +* column-vector containing the output coordinates. If appropriate, +* the inverse transformation may also be performed. + +* Parameters: +c nin +f NIN = INTEGER (Given) +* The number of input coordinates, which determines the number +* of columns in the matrix. +c nout +f NOUT = INTEGER (Given) +* The number of output coordinates, which determines the number +* of rows in the matrix. +c form +f FORM = INTEGER (Given) +* An integer which indicates the form in which the matrix +* elements will be supplied. +* +c A value of zero indicates that a full "nout" x "nin" matrix +f A value of zero indicates that a full NOUT x NIN matrix +c of values will be supplied via the "matrix" parameter +f of values will be supplied via the MATRIX argument +* (below). In this case, the elements should be given in row +* order (the elements of the first row, followed by the +* elements of the second row, etc.). +* +* A value of 1 indicates that only the diagonal elements of the +* matrix will be supplied, and that all others should be +c zero. In this case, the elements of "matrix" should contain +f zero. In this case, the elements of MATRIX should contain +* only the diagonal elements, stored consecutively. +* +* A value of 2 indicates that a "unit" matrix is required, +* whose diagonal elements are set to unity (with all other +c elements zero). In this case, the "matrix" parameter is +c ignored and a NULL pointer may be supplied. +f elements zero). In this case, the MATRIX argument is not used. +c matrix +f MATRIX( * ) = DOUBLE PRECISION (Given) +* The array of matrix elements to be used, stored according to +c the value of "form". +f the value of FORM. +c options +f OPTIONS = CHARACTER * ( * ) (Given) +c Pointer to a null-terminated string containing an optional +c comma-separated list of attribute assignments to be used for +c initialising the new MatrixMap. The syntax used is identical to +c that for the astSet function and may include "printf" format +c specifiers identified by "%" symbols in the normal way. +f A character string containing an optional comma-separated +f list of attribute assignments to be used for initialising the +f new MatrixMap. The syntax used is identical to that for the +f AST_SET routine. +c ... +c If the "options" string contains "%" format specifiers, then +c an optional list of additional arguments may follow it in +c order to supply values to be substituted for these +c specifiers. The rules for supplying these are identical to +c those for the astSet function (and for the C "printf" +c function). +f STATUS = INTEGER (Given and Returned) +f The global status. + +* Returned Value: +c astMatrixMap() +f AST_MATRIXMAP = INTEGER +* A pointer to the new MatrixMap. + +* Notes: +* - In general, a MatrixMap's forward transformation will always +* be available (as indicated by its TranForward attribute), but +* its inverse transformation (TranInverse attribute) will only be +* available if the associated matrix is square and non-singular. +* - As an exception to this, the inverse transformation is always +* available if a unit or diagonal matrix is specified. In this +* case, if the matrix is not square, one or more of the input +* coordinate values may not be recoverable from a set of output +* coordinates. Any coordinates affected in this way will simply be +* set to the value zero. +* - A null Object pointer (AST__NULL) will be returned if this +c function is invoked with the AST error status set, or if it +f function is invoked with STATUS set to an error value, or if it +* should fail for any reason. + +* Status Handling: +* The protected interface to this function includes an extra +* parameter at the end of the parameter list descirbed above. This +* parameter is a pointer to the integer inherited status +* variable: "int *status". + +*-- +*/ + +/* Local Variables: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstMatrixMap *new; /* Pointer to new MatrixMap */ + va_list args; /* Variable argument list */ + +/* Check the global status. */ + if ( !astOK ) return NULL; + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(NULL); + +/* Initialise the MatrixMap, allocating memory and initialising the + virtual function table as well if necessary. */ + new = astInitMatrixMap( NULL, sizeof( AstMatrixMap ), !class_init, + &class_vtab, "MatrixMap", nin, nout, form, matrix); + +/* If successful, note that the virtual function table has been + initialised. */ + if ( astOK ) { + class_init = 1; + +/* Obtain the variable argument list and pass it along with the options string + to the astVSet method to initialise the new MatrixMap's attributes. */ + va_start( args, status ); + astVSet( new, options, NULL, args ); + va_end( args ); + +/* If an error occurred, clean up by deleting the new object. */ + if ( !astOK ) new = astDelete( new ); + } + +/* Return a pointer to the new MatrixMap. */ + return new; +} + +AstMatrixMap *astMatrixMapId_( int nin, int nout, int form, const double matrix[], + const char *options, ... ) { +/* +* Name: +* astMatrixMapId_ + +* Purpose: +* Create a MatrixMap. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* AstMatrixMap *astMatrixMapId_( int nin, int nout, int form, +* const double matrix[], const char *options, +* ... ) + +* Class Membership: +* MatrixMap constructor. + +* Description: +* This function implements the external (public) interface to the +* astMatrixMap constructor function. It returns an ID value (instead +* of a true C pointer) to external users, and must be provided +* because astMatrixMap_ has a variable argument list which cannot be +* encapsulated in a macro (where this conversion would otherwise +* occur). +* +* The variable argument list also prevents this function from +* invoking astMatrixMap_ directly, so it must be a re-implementation +* of it in all respects, except for the final conversion of the +* result to an ID value. + +* Parameters: +* As for astMatrixMap_. + +* Returned Value: +* The ID value associated with the new MatrixMap. +*/ + +/* Local Variables: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstMatrixMap *new; /* Pointer to new MatrixMap */ + va_list args; /* Variable argument list */ + int *status; /* Pointer to inherited status value */ + +/* Get a pointer to the inherited status value. */ + status = astGetStatusPtr; + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(NULL); + +/* Check the global status. */ + if ( !astOK ) return NULL; + +/* Initialise the MatrixMap, allocating memory and initialising the + virtual function table as well if necessary. */ + new = astInitMatrixMap( NULL, sizeof( AstMatrixMap ), !class_init, &class_vtab, + "MatrixMap", nin, nout, form, matrix ); + +/* If successful, note that the virtual function table has been + initialised. */ + if ( astOK ) { + class_init = 1; + +/* Obtain the variable argument list and pass it along with the options string + to the astVSet method to initialise the new MatrixMap's attributes. */ + va_start( args, options ); + astVSet( new, options, NULL, args ); + va_end( args ); + +/* If an error occurred, clean up by deleting the new object. */ + if ( !astOK ) new = astDelete( new ); + } + +/* Return an ID value for the new MatrixMap. */ + return astMakeId( new ); +} + +AstMatrixMap *astInitMatrixMap_( void *mem, size_t size, int init, + AstMatrixMapVtab *vtab, const char *name, + int nin, int nout, int form, + const double *matrix, int *status ) { +/* +*+ +* Name: +* astInitMatrixMap + +* Purpose: +* Initialise a MatrixMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "matrixmap.h" +* AstMatrixMap *astInitMatrixMap( void *mem, size_t size, int init, +* AstMatrixMapVtab *vtab, const char *name, +* int nin, int nout, int form, +* const double *matrix ) + +* Class Membership: +* MatrixMap initialiser. + +* Description: +* This function is provided for use by class implementations to initialise +* a new MatrixMap object. It allocates memory (if necessary) to accommodate +* the MatrixMap plus any additional data associated with the derived class. +* It then initialises a MatrixMap structure at the start of this memory. If +* the "init" flag is set, it also initialises the contents of a virtual +* function table for a MatrixMap at the start of the memory passed via the +* "vtab" parameter. + +* Parameters: +* mem +* A pointer to the memory in which the MatrixMap is to be initialised. +* This must be of sufficient size to accommodate the MatrixMap data +* (sizeof(MatrixMap)) plus any data used by the derived class. If a value +* of NULL is given, this function will allocate the memory itself using +* the "size" parameter to determine its size. +* size +* The amount of memory used by the MatrixMap (plus derived class data). +* This will be used to allocate memory if a value of NULL is given for +* the "mem" parameter. This value is also stored in the MatrixMap +* structure, so a valid value must be supplied even if not required for +* allocating memory. +* init +* A logical flag indicating if the MatrixMap's virtual function table is +* to be initialised. If this value is non-zero, the virtual function +* table will be initialised by this function. +* vtab +* Pointer to the start of the virtual function table to be associated +* with the new MatrixMap. +* name +* Pointer to a constant null-terminated character string which contains +* the name of the class to which the new object belongs (it is this +* pointer value that will subsequently be returned by the astGetClass +* method). +* nin +* The number of input coordinate values per point. This is the +* same as the number of columns in the matrix. +* nout +* The number of output coordinate values per point. This is the +* same as the number of rows in the matrix. +* form +* If "form" is 2 or larger, then a unit MatrixMap is created. In this +* case "matrix" is ignored and can be supplied as NULL. If "form" is +* 1, then a diagonal MatrixMap is created. In this case, the number of +* values in "matrix" should be equal to the minimum of nin and nout, +* and "matrix" should contain the corresponding diagonal terms, in row +* order. If "form" is 0 or less, then a full MatrixMap is created, and +* "matrix" should contain all nin*nout element values. +* matrix +* A pointer to an array of matrix element values. The values should be +* supplied in row order. The content of this array is determined by +* "form". If a full MatrixMap is to be created then the array starts +* with (row 1, column 1), then comes (row 1, column 2), up to (row 1, +* column nin), then (row 2, column 1), (row 2, column 2), and so on, +* finishing with (row nout, column nin) ). An error is reported if a +* NULL value is supplied unless "form" is 2 or more. + +* Returned Value: +* A pointer to the new MatrixMap. + +* Notes: +* - A null pointer will be returned if this function is invoked with the +* global error status set, or if it should fail for any reason. +*- +*/ + +/* Local Variables: */ + AstMatrixMap *new; /* Pointer to new MatrixMap */ + double *fmat; /* Pointer to the forward matrix */ + double *imat; /* Pointer to the inverse matrix */ + int i; /* Loop count */ + int nel; /* No. of elements in matrix array */ + int nuse; /* Number of usable matrix elements */ + int used_form; /* Form limited to 0, 1 or 2 */ + +/* Check the global status. */ + if ( !astOK ) return NULL; + +/* If necessary, initialise the virtual function table. */ + if ( init ) astInitMatrixMapVtab( vtab, name ); + +/* Initialise. */ + new = NULL; + +/* Report an error if a NULL matrix was supplied, unless a unit MatrixMap + has been requested. */ + if( form < 2 && !matrix ){ + astError( AST__MTRMT, "astInitMatrixMap(%s): NULL matrix supplied.", status, + name ); + + } else { + +/* Initialise a Mapping structure (the parent class) as the first component + within the MatrixMap structure, allocating memory if necessary. Specify that + the Mapping should be defined in both the forward and inverse directions. */ + new = (AstMatrixMap *) astInitMapping( mem, size, 0, + (AstMappingVtab *) vtab, name, + nin, nout, 1, 1 ); + + if ( astOK ) { + +/* Initialise the MatrixMap data. */ +/* ---------------------------- */ +/* If a unit MatrixMap is being created, then no additional storage is + required. */ + if( form > 1 ){ + nel = 0; + used_form = UNIT; + +/* If a diagonal MatrixMap is being created, then memory is needed to hold + the diagonal terms. */ + } else if( form == 1 ){ + if( nin < nout ){ + nel = nin; + } else { + nel = nout; + } + used_form = DIAGONAL; + +/* If a full MatrixMap is being created, then memory is needed to hold + all the terms. */ + } else { + nel = nin*nout ; + used_form = FULL; + } + +/* Allocate memory for the forward matrix, storing the supplied matrix + values in it. */ + fmat = (double *) astStore( NULL, (void *) matrix, + sizeof(double)*(size_t)nel ); + +/* Replace any NaNs by AST__BAD and count the number of usable values. */ + if( nel > 0 ) { + nuse = 0; + for( i = 0; i < nel; i++ ) { + if( !astISFINITE(fmat[ i ]) ) { + fmat[ i ] = AST__BAD; + } else if( fmat[ i ] != AST__BAD ) { + nuse++; + } + } + +/* Report an error if there are no usable values. */ + if( nuse == 0 && astOK ) { + astError( AST__MTRMT, "astInitMatrixMap(%s): Supplied matrix " + "contains only bad values.", status, name ); + } + } + +/* Create an inverse matrix if possible. */ + imat = InvertMatrix( used_form, nout, nin, fmat, status ); + +/* Store the matrix arrays. */ + new->form = used_form; + new->f_matrix = fmat; + new->i_matrix = imat; + +/* Attempt to compress the MatrixMap into DIAGONAL or UNIT form. */ + CompressMatrix( new, status ); + +/* If an error occurred, clean up by deleting the new MatrixMap. */ + if ( !astOK ) new = astDelete( new ); + } + } + +/* Return a pointer to the new MatrixMap. */ + return new; +} + +AstMatrixMap *astLoadMatrixMap_( void *mem, size_t size, + AstMatrixMapVtab *vtab, const char *name, + AstChannel *channel, int *status ) { +/* +*+ +* Name: +* astLoadMatrixMap + +* Purpose: +* Load a MatrixMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "matrixmap.h" +* AstMatrixMap *astLoadMatrixMap( void *mem, size_t size, +* AstMatrixMapVtab *vtab, const char *name, +* AstChannel *channel ) + +* Class Membership: +* MatrixMap loader. + +* Description: +* This function is provided to load a new MatrixMap using data read +* from a Channel. It first loads the data used by the parent class +* (which allocates memory if necessary) and then initialises a +* MatrixMap structure in this memory, using data read from the input +* Channel. +* +* If the "init" flag is set, it also initialises the contents of a +* virtual function table for a MatrixMap at the start of the memory +* passed via the "vtab" parameter. + + +* Parameters: +* mem +* A pointer to the memory into which the MatrixMap is to be +* loaded. This must be of sufficient size to accommodate the +* MatrixMap data (sizeof(MatrixMap)) plus any data used by derived +* classes. If a value of NULL is given, this function will +* allocate the memory itself using the "size" parameter to +* determine its size. +* size +* The amount of memory used by the MatrixMap (plus derived class +* data). This will be used to allocate memory if a value of +* NULL is given for the "mem" parameter. This value is also +* stored in the MatrixMap structure, so a valid value must be +* supplied even if not required for allocating memory. +* +* If the "vtab" parameter is NULL, the "size" value is ignored +* and sizeof(AstMatrixMap) is used instead. +* vtab +* Pointer to the start of the virtual function table to be +* associated with the new MatrixMap. If this is NULL, a pointer +* to the (static) virtual function table for the MatrixMap class +* is used instead. +* name +* Pointer to a constant null-terminated character string which +* contains the name of the class to which the new object +* belongs (it is this pointer value that will subsequently be +* returned by the astGetClass method). +* +* If the "vtab" parameter is NULL, the "name" value is ignored +* and a pointer to the string "MatrixMap" is used instead. + +* Returned Value: +* A pointer to the new MatrixMap. + +* Notes: +* - A null pointer will be returned if this function is invoked +* with the global error status set, or if it should fail for any +* reason. +*- +*/ + +#define KEY_LEN 50 /* Maximum length of a keyword */ + + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ +/* Local Variables: */ + AstMatrixMap *new; /* Pointer to the new MatrixMap */ + char buff[ KEY_LEN + 1 ]; /* Buffer for keyword string */ + const char *form; /* String form */ + int def; /* Is the matrix defined? */ + int el; /* Element index */ + int nel; /* No. of elements in the matrix */ + int nin; /* No. of input coords */ + int nout; /* No. of output coords */ + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(channel); + +/* Initialise. */ + new = NULL; + +/* Check the global error status. */ + if ( !astOK ) return new; + +/* If a NULL virtual function table has been supplied, then this is + the first loader to be invoked for this MatrixMap. In this case the + MatrixMap belongs to this class, so supply appropriate values to be + passed to the parent class loader (and its parent, etc.). */ + if ( !vtab ) { + size = sizeof( AstMatrixMap ); + vtab = &class_vtab; + name = "MatrixMap"; + +/* If required, initialise the virtual function table for this class. */ + if ( !class_init ) { + astInitMatrixMapVtab( vtab, name ); + class_init = 1; + } + } + +/* Invoke the parent class loader to load data for all the ancestral + classes of the current one, returning a pointer to the resulting + partly-built MatrixMap. */ + new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name, + channel ); + + if ( astOK ) { + +/* Read input data. */ +/* ================ */ +/* Request the input Channel to read all the input data appropriate to + this class into the internal "values list". */ + astReadClassData( channel, "MatrixMap" ); + +/* Now obtain the Matrix storage form from this list. */ + form = astReadString( channel, "form", Form[FULL] ); + new->form = FindString( 3, Form, form, "the MatrixMap component 'Form'", + "astRead", astGetClass( channel ), status ); + form = astFree( (void *) form ); + +/* Find the number of elements stored for each matrix. */ + nin = astGetNin( (AstMapping *) new ); + nout = astGetNout( (AstMapping *) new ); + + if( new->form == FULL ){ + nel = nin*nout; + + } else if( new->form == DIAGONAL ){ + nel = astMIN( nin, nout ); + + } else { + nel = 0; + } + +/* Allocate memory to hold the forward matrix. */ + new->f_matrix = (double *) astMalloc( sizeof(double)*(size_t)nel ); + +/* Now read the other data items from the list and use them to + initialise the appropriate instance variable(s) for this class. */ + +/* The forward matrix. */ + if( new->f_matrix ){ + def = 0; + + for( el = 0; el < nel; el++ ){ + (void) sprintf( buff, "m%d", el ); + (new->f_matrix)[ el ] = astReadDouble( channel, buff, AST__BAD ); + if( (new->f_matrix)[ el ] != AST__BAD ) def = 1; + } + +/* Store a NULL pointer if no elements of the matrix were found. */ + if( !def ) new->f_matrix = (double *) astFree( (void *) new->f_matrix ); + + } + +/* The inverse matrix. */ + new->i_matrix = (double *) astMalloc( sizeof(double)*(size_t)nel ); + if( new->i_matrix ){ + def = 0; + + for( el = 0; el < nel; el++ ){ + (void) sprintf( buff, "im%d", el ); + (new->i_matrix)[ el ] = astReadDouble( channel, buff, AST__BAD ); + if( (new->i_matrix)[ el ] != AST__BAD ) def = 1; + } + +/* If no elements of the matrix were found, create an inverse matrix if + possible, otherwise store a NULL pointer. Note, prior to AST 8.6.3, the + inverse matrix was not included in the dump - it was always recalculated + using InvertMatrix, but this led to small round-trip errors in cases, + where the original inverse matrix was not created using InvertMatrix + (e.g. was created by astMtrRot). */ + if( !def ) { + new->i_matrix = (double *) astFree( (void *) new->i_matrix ); + if( new->f_matrix ){ + new->i_matrix = InvertMatrix( new->form, nout, nin, new->f_matrix, status ); + } else { + new->i_matrix = NULL; + } + } + } + +/* If an error occurred, clean up by deleting the new MatrixMap. */ + if ( !astOK ) new = astDelete( new ); + } + +/* Return the new MatrixMap pointer. */ + return new; + +/* Undefine macros local to this function. */ +#undef KEY_LEN +} + +/* Virtual function interfaces. */ +/* ============================ */ +/* These provide the external interface to the virtual functions defined by + this class. Each simply checks the global error status and then locates and + executes the appropriate member function, using the function pointer stored + in the object's virtual function table (this pointer is located using the + astMEMBER macro defined in "object.h"). + + Note that the member function may not be the one defined here, as it may + have been over-ridden by a derived class. However, it should still have the + same interface. */ + +AstMatrixMap *astMtrRot_( AstMatrixMap *this, double theta, + const double axis[], int *status ){ + if( !astOK ) return NULL; + return (**astMEMBER(this,MatrixMap,MtrRot))( this, theta, axis, status ); +} + +AstMatrixMap *astMtrMult_( AstMatrixMap *this, AstMatrixMap *a, int *status ){ + if( !astOK ) return NULL; + return (**astMEMBER(this,MatrixMap,MtrMult))( this, a, status ); +} + + + + |