diff options
author | William Joye <wjoye@cfa.harvard.edu> | 2019-05-10 18:53:10 (GMT) |
---|---|---|
committer | William Joye <wjoye@cfa.harvard.edu> | 2019-05-10 18:53:10 (GMT) |
commit | a213dadabce2b2e02eca8376765fa425f01792f5 (patch) | |
tree | 9a21116c977447511fed092e80415ae664a92945 /ast/polymap.c | |
parent | 9ef925b353e71b6b8abd23b59ca1adef488e3a3b (diff) | |
download | blt-a213dadabce2b2e02eca8376765fa425f01792f5.zip blt-a213dadabce2b2e02eca8376765fa425f01792f5.tar.gz blt-a213dadabce2b2e02eca8376765fa425f01792f5.tar.bz2 |
upgrade ast 8.7.1
Diffstat (limited to 'ast/polymap.c')
-rw-r--r-- | ast/polymap.c | 6359 |
1 files changed, 6359 insertions, 0 deletions
diff --git a/ast/polymap.c b/ast/polymap.c new file mode 100644 index 0000000..f03a631 --- /dev/null +++ b/ast/polymap.c @@ -0,0 +1,6359 @@ +/* +*class++ +* Name: +* PolyMap + +* Purpose: +* Map coordinates using polynomial functions. + +* Constructor Function: +c astPolyMap +f AST_POLYMAP + +* Description: +* A PolyMap is a form of Mapping which performs a general polynomial +* transformation. Each output coordinate is a polynomial function of +* all the input coordinates. The coefficients are specified separately +* for each output coordinate. The forward and inverse transformations +* are defined independantly by separate sets of coefficients. If no +* inverse transformation is supplied, the default behaviour is to use +* an iterative method to evaluate the inverse based only on the forward +* transformation (see attribute IterInverse). + +* Inheritance: +* The PolyMap class inherits from the Mapping class. + +* Attributes: +* In addition to those attributes common to all Mappings, every +* PolyMap also has the following attributes: +* +* - IterInverse: Provide an iterative inverse transformation? +* - NiterInverse: Maximum number of iterations for iterative inverse +* - TolInverse: Target relative error for iterative inverse + +* Functions: +c In addition to those functions applicable to all Objects, the +c following functions may also be applied to all Mappings: +f In addition to those routines applicable to all Objects, the +f following routines may also be applied to all Mappings: +* +c - astPolyCoeffs: Retrieve the coefficients of a PolyMap transformation +c - astPolyTran: Fit a PolyMap inverse or forward transformation +f - AST_POLYCOEFFS: Retrieve the coefficients of a PolyMap transformation +f - AST_POLYTRAN: Fit a PolyMap inverse or forward transformation + +* Copyright: +* Copyright (C) 2017 East Asian Observatory. +* Copyright (C) 1997-2006 Council for the Central Laboratory of the +* Research Councils +* Copyright (C) 2011 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) + +* History: +* 27-SEP-2003 (DSB): +* Original version. +* 13-APR-2005 (DSB): +* Changed the keys used by the Dump/astLoadPolyMap functions. They +* used to exceed 8 characters and consequently caused problems for +* FitsChans. +* 20-MAY-2005 (DSB): +* Correct the indexing of keywords produced in the Dump function. +* 20-APR-2006 (DSB): +* Guard against undefined transformations in Copy. +* 10-MAY-2006 (DSB): +* Override astEqual. +* 4-JUL-2008 (DSB): +* Fixed loop indexing problems in Equal function. +* 27-MAY-2011 (DSB): +* Added public method astPolyTran. +* 18-JUL-2011 (DSB): +* - Added attributes IterInverse, NiterInverse and TolInverse. +* - Do not report an error if astPolyTran fails to fit an inverse. +* 15-OCT-2011 (DSB): +* Improve argument checking and error reporting in PolyTran +* 8-MAY-2014 (DSB): +* Move to using CMinPack for minimisations. +* 11-NOV-2016 (DSB): +* - Fix bug in MapMerge that could cause a seg fault. It did not +* check that the PolyMap had a defined transformation before accessing +* the transformation's coefficient array. +* - Fix similar bugs in Equal that could cause seg faults. +* 15-MAR-2017 (DSB): +* - Change the GetTranForward and GetTranInverse functions so that they +* take into account the state of the Invert attribute. +* - Improve docs for the IterInverse attribute to explain that the +* inverse transformation replaced is always the original inverse +* transformation, as defined by the arguments supplied to the PolyMap +* constructor, regardless of the state of the Invert attribute. +* 17-MAR-2017 (DSB): +* - Add the astPolyCoeffs method. +* 30-MAR-2017 (DSB): +* Modify the astPolyTran method so that it can be used by the +* ChebyMap class to determine new transformations implemented as +* Chebyshev polynomials. +* 27-JUN-2017 (DSB): +* In SamplePoly1D/2D ensure the final sample on each axis does not +* go above the supplied upper bound. This can happen due to rounding +* error. This is important for ChebyMaps since points outside the +* bounds are set bad when transformed using a ChebyMap, causing NaNs +* to be generated in lmder1 (cminpack minimisation function). +* 3-JUL-2017 (DSB): +* Within FitPoly1D and FitPoly2D, use an initial guess that represents +* a unit mapping between normalised input and output values, rather +* than unnormalised PolyMap values. This is in case the PolyMap being +* fitted includes a change of scale (e.g. the PolyMap input is in "mm" +* but the output is in "rads" and includes some large scaling factor +* to do the conversion). +* 9-JAN-2018 (DSB): +* Correct Transform to take account of AST__BAD coeffs correctly. +* Previously bad coeffs could generate NaN output values. +* 27-APR-2018 (DSB): +* When calculating the iterative inverse, use an initial guess based +* on the linear truncation of the PolyMap rather than a UnitMap. +* 3-MAY-2018 (DSB): +* Clear the IterInverse attribute in the PolyMap returned by +* astPolyTran (why would you want to use the iterative inverse if +* you have just created an inverse fit?). +* 4-MAY-2018 (DSB): +* Fix various places where there was confusion over whether the +* "forward transformation" refers to the forward transformation +* of the original uninverted PolyMap, or the current forward +* transformation of the PolyMap (i.e. taking the "Invert" flag into +* account). +*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 PolyMap + +/* 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 "matrixmap.h" /* Matrix multiplication mappings */ +#include "shiftmap.h" /* Shift of origin mappings */ +#include "mapping.h" /* Coordinate mappings (parent class) */ +#include "cmpmap.h" /* Compound mappings */ +#include "polymap.h" /* Interface definition for this class */ +#include "unitmap.h" /* Unit mappings */ +#include "cminpack/cminpack.h" /* Levenberg - Marquardt minimization */ +#include "pal.h" /* SLALIB function definitions */ + +/* 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> +#include <limits.h> +#include <float.h> + +/* Module Variables. */ +/* ================= */ + +/* Address of this static variable is used as a unique identifier for + member of this class. */ +static int class_check; + +/* Pointers to parent class methods which are extended by this class. */ +static AstPointSet *(* parent_transform)( AstMapping *, AstPointSet *, int, AstPointSet *, int * ); +static const char *(* parent_getattrib)( AstObject *, const char *, int * ); +static int (* parent_testattrib)( AstObject *, const char *, int * ); +static void (* parent_clearattrib)( AstObject *, const char *, int * ); +static void (* parent_setattrib)( AstObject *, const char *, int * ); +static int (* parent_getobjsize)( AstObject *, int * ); + +#if defined(THREAD_SAFE) +static int (* parent_managelock)( AstObject *, int, int, AstObject **, int * ); +#endif + + +#ifdef THREAD_SAFE +/* Define how to initialise thread-specific globals. */ +#define GLOBAL_inits \ + globals->Class_Init = 0; \ + globals->GetAttrib_Buff[ 0 ] = 0; + +/* Create the function that initialises global data for this module. */ +astMAKE_INITGLOBALS(PolyMap) + +/* Define macros for accessing each item of thread specific global data. */ +#define class_init astGLOBAL(PolyMap,Class_Init) +#define class_vtab astGLOBAL(PolyMap,Class_Vtab) +#define getattrib_buff astGLOBAL(LutMap,GetAttrib_Buff) + +#include <pthread.h> + + +#else + +static char getattrib_buff[ 101 ]; + +/* Define the class virtual function table and its initialisation flag + as static variables. */ +static AstPolyMapVtab 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. */ +AstPolyMap *astPolyMapId_( int, int, int, const double[], int, const double[], const char *, ... ); + + +/* Prototypes for Private Member Functions. */ +/* ======================================== */ +static AstMapping *LinearGuess( AstPolyMap *, int * ); +static AstPointSet *Transform( AstMapping *, AstPointSet *, int, AstPointSet *, int * ); +static AstPolyMap **GetJacobian( AstPolyMap *, int * ); +static AstPolyMap *PolyTran( AstPolyMap *, int, double, double, int, const double *, const double *, int * ); +static double **SamplePoly1D( AstPolyMap *, int, double **, double, double, int, int *, double[2], int * ); +static double **SamplePoly2D( AstPolyMap *, int, double **, const double *, const double *, int, int *, double[4], int * ); +static double *FitPoly1D( AstPolyMap *, int, int, double, int, double **, double[2], int *, double *, int * ); +static double *FitPoly2D( AstPolyMap *, int, int, double, int, double **, double[4], int *, double *, int * ); +static int Equal( AstObject *, AstObject *, int * ); +static int GetObjSize( AstObject *, int * ); +static int GetTranForward( AstMapping *, int * ); +static int GetTranInverse( AstMapping *, int * ); +static int MPFunc1D( void *, int, int, const double *, double *, double *, int, int ); +static int MPFunc2D( void *, int, int, const double *, double *, double *, int, int ); +static int MapMerge( AstMapping *, int, int, int *, AstMapping ***, int **, int * ); +static int ReplaceTransformation( AstPolyMap *, int, double, double, int, const double *, const double *, int * ); +static void Copy( const AstObject *, AstObject *, int * ); +static void Delete( AstObject *obj, int * ); +static void Dump( AstObject *, AstChannel *, int * ); +static void FitPoly1DInit( AstPolyMap *, int, double **, AstMinPackData *, double *, int *); +static void FitPoly2DInit( AstPolyMap *, int, double **, AstMinPackData *, double *, int *); +static void FreeArrays( AstPolyMap *, int, int * ); +static void IterInverse( AstPolyMap *, AstPointSet *, AstPointSet *, int * ); +static void LMFunc1D( const double *, double *, int, int, void * ); +static void LMFunc2D( const double *, double *, int, int, void * ); +static void LMJacob1D( const double *, double *, int, int, void * ); +static void LMJacob2D( const double *, double *, int, int, void * ); +static void PolyCoeffs( AstPolyMap *, int, int, double *, int *, int * ); +static void PolyPowers( AstPolyMap *, double **, int, const int *, double **, int, int, int * ); +static void StoreArrays( AstPolyMap *, int, int, const double *, int * ); + +#if defined(THREAD_SAFE) +static int ManageLock( AstObject *, int, int, AstObject **, int * ); +#endif + +static const char *GetAttrib( AstObject *, const char *, int * ); +static int TestAttrib( AstObject *, const char *, int * ); +static void ClearAttrib( AstObject *, const char *, int * ); +static void SetAttrib( AstObject *, const char *, int * ); + +static int GetIterInverse( AstPolyMap *, int * ); +static int TestIterInverse( AstPolyMap *, int * ); +static void ClearIterInverse( AstPolyMap *, int * ); +static void SetIterInverse( AstPolyMap *, int, int * ); + +static int GetNiterInverse( AstPolyMap *, int * ); +static int TestNiterInverse( AstPolyMap *, int * ); +static void ClearNiterInverse( AstPolyMap *, int * ); +static void SetNiterInverse( AstPolyMap *, int, int * ); + +static double GetTolInverse( AstPolyMap *, int * ); +static int TestTolInverse( AstPolyMap *, int * ); +static void ClearTolInverse( AstPolyMap *, int * ); +static void SetTolInverse( AstPolyMap *, double, int * ); + + +/* Member functions. */ +/* ================= */ + +static void ClearAttrib( AstObject *this_object, const char *attrib, int *status ) { +/* +* Name: +* ClearAttrib + +* Purpose: +* Clear an attribute value for a PolyMap. + +* Type: +* Private function. + +* Synopsis: +* #include "polymap.h" +* void ClearAttrib( AstObject *this, const char *attrib, int *status ) + +* Class Membership: +* PolyMap member function (over-rides the astClearAttrib protected +* method inherited from the Mapping class). + +* Description: +* This function clears the value of a specified attribute for a +* PolyMap, so that the default value will subsequently be used. + +* Parameters: +* this +* Pointer to the PolyMap. +* attrib +* Pointer to a null-terminated string specifying the attribute +* name. This should be in lower case with no surrounding white +* space. +* status +* Pointer to the inherited status variable. +*/ + +/* Local Variables: */ + AstPolyMap *this; /* Pointer to the PolyMap structure */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain a pointer to the PolyMap structure. */ + this = (AstPolyMap *) this_object; + +/* Check the attribute name and clear the appropriate attribute. */ + +/* IterInverse. */ +/* ------------ */ + if ( !strcmp( attrib, "iterinverse" ) ) { + astClearIterInverse( this ); + +/* NiterInverse. */ +/* ------------- */ + } else if ( !strcmp( attrib, "niterinverse" ) ) { + astClearNiterInverse( this ); + +/* TolInverse. */ +/* ----------- */ + } else if ( !strcmp( attrib, "tolinverse" ) ) { + astClearTolInverse( this ); + +/* If the attribute is still not recognised, pass it on to the parent + method for further interpretation. */ + } else { + (*parent_clearattrib)( this_object, attrib, status ); + } +} + +static int Equal( AstObject *this_object, AstObject *that_object, int *status ) { +/* +* Name: +* Equal + +* Purpose: +* Test if two PolyMaps are equivalent. + +* Type: +* Private function. + +* Synopsis: +* #include "polymap.h" +* int Equal( AstObject *this, AstObject *that, int *status ) + +* Class Membership: +* PolyMap 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 PolyMaps are equivalent. + +* Parameters: +* this +* Pointer to the first Object (a PolyMap). +* that +* Pointer to the second Object. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* One if the PolyMaps 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: */ + AstPolyMap *that; + AstPolyMap *this; + int i, j, k; + int nin; + int nout; + int result; + int tmp; + +/* Initialise. */ + result = 0; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Obtain pointers to the two PolyMap structures. */ + this = (AstPolyMap *) this_object; + that = (AstPolyMap *) that_object; + +/* Check the second object is a PolyMap. We know the first is a + PolyMap since we have arrived at this implementation of the virtual + function. */ + if( astIsAPolyMap( that ) ) { + +/* Get the number of inputs and outputs and check they are the same for both. */ + nin = astGetNin( this ); + nout = astGetNout( this ); + if( astGetNin( that ) == nin && astGetNout( that ) == nout ) { + +/* If the Invert flags for the two PolyMaps differ, it may still be possible + for them to be equivalent. First compare the PolyMaps if their Invert + flags are the same. In this case all the attributes of the two PolyMaps + must be identical. */ + if( astGetInvert( this ) == astGetInvert( that ) ) { + +/* Assume they are the same until we find a difference. */ + result = 1; + +/* The "_f" and "_i" suffixes in the PolyMap structure array names refer + to the forward and inverse transformations of the original uninverted + PolyMap. So we need to ensure the "nout" and "nin" values also refer + to the uninverted values. So if the PolyMaps are inverted, swap nout + and nin. */ + if( astGetInvert( this ) ) { + tmp = nin; + nin = nout; + nout = tmp; + } + +/* Check properties of the forward transformation. */ + if( this->ncoeff_f && that->ncoeff_f ) { + for( i = 0; i < nout && result; i++ ) { + if( this->ncoeff_f[ i ] != that->ncoeff_f[ i ] ){ + result = 0; + } + } + } else if( this->ncoeff_f || that->ncoeff_f ) { + result = 0; + } + + if( this->mxpow_f && that->mxpow_f ) { + for( i = 0; i < nout && result; i++ ) { + if( this->mxpow_f[ i ] != that->mxpow_f[ i ] ){ + result = 0; + } + } + } else if( this->mxpow_f || that->mxpow_f ) { + result = 0; + } + + if( this->coeff_f && that->coeff_f ) { + for( i = 0; i < nout && result; i++ ) { + for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) { + if( !astEQUAL( this->coeff_f[ i ][ j ], + that->coeff_f[ i ][ j ] ) ) { + result = 0; + } + } + } + } else if( this->coeff_f || that->coeff_f ) { + result = 0; + } + + if( this->power_f && that->power_f ) { + for( i = 0; i < nout && result; i++ ) { + for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) { + for( k = 0; k < nin && result; k++ ) { + if( this->power_f[ i ][ j ][ k ] != + that->power_f[ i ][ j ][ k ] ) { + result = 0; + } + } + } + } + } else if( this->power_f || that->power_f ) { + result = 0; + } + +/* Check properties of the inverse transformation. */ + if( this->ncoeff_i && that->ncoeff_i ) { + for( i = 0; i < nout && result; i++ ) { + if( this->ncoeff_i[ i ] != that->ncoeff_i[ i ] ){ + result = 0; + } + } + } else if( this->ncoeff_i || that->ncoeff_i ) { + result = 0; + } + + if( this->mxpow_i && that->mxpow_i ) { + for( i = 0; i < nout && result; i++ ) { + if( this->mxpow_i[ i ] != that->mxpow_i[ i ] ){ + result = 0; + } + } + } else if( this->mxpow_i || that->mxpow_i ) { + result = 0; + } + + if( this->coeff_f && that->coeff_f ) { + for( i = 0; i < nout && result; i++ ) { + for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) { + if( !astEQUAL( this->coeff_f[ i ][ j ], + that->coeff_f[ i ][ j ] ) ) { + result = 0; + } + } + } + } else if( this->coeff_f || that->coeff_f ) { + result = 0; + } + + if( this->power_f && that->power_f ) { + for( i = 0; i < nout && result; i++ ) { + for( j = 0; j < this->ncoeff_f[ i ] && result; j++ ) { + for( k = 0; k < nin && result; k++ ) { + if( this->power_f[ i ][ j ][ k ] != + that->power_f[ i ][ j ][ k ] ) { + result = 0; + } + } + } + } + } else if( this->power_f || that->power_f ) { + result = 0; + } + +/* If the Invert flags for the two PolyMaps differ, the attributes of the two + PolyMaps must be inversely related to each other. */ + } else { + +/* In the specific case of a PolyMap, Invert flags must be equal. */ + result = 0; + + } + } + } + +/* If an error occurred, clear the result value. */ + if ( !astOK ) result = 0; + +/* Return the result, */ + return result; +} + +static double *FitPoly1D( AstPolyMap *this, int forward, int nsamp, double acc, + int order, double **table, double scales[2], int *ncoeff, + double *racc, int *status ){ +/* +* Name: +* FitPoly1D + +* Purpose: +* Fit a (1-in,1-out) polynomial to a supplied set of data. + +* Type: +* Private function. + +* Synopsis: +* double *FitPoly1D( AstPolyMap *this, int forward, int nsamp, double acc, +* int order, double **table, double scales[2], int *ncoeff, +* double *racc, int *status ) + +* Description: +* This function fits a least squares 1D polynomial curve to the +* positions in a supplied table, and returns the coefficients of a +* PolyMap to describe the fit. For the purposes of this function, +* the input to the fit is refered to as x1 and the output as y1. So +* the returned coefficients describe a PolyMap with forward +* transformation: +* +* y1 = P1( x1 ) + +* Parameters: +* this +* Pointer to the PolyMap. +* forward +* Non-zero if the forward transformation of "this" is being +* replaced. Zero if the inverse transformation of "this" is being +* replaced. +* nsamp +* The number of (x1,y1) positions in the supplied table. +* acc +* The required accuracy, expressed as a geodesic distance within +* the polynomials output space (not the normalised tabular space). +* order +* The maximum power (plus one) of x1 within P1. So for instance, a +* value of 3 refers to a quadratic polynomial. +* table +* Pointer to an array of 2 pointers. Each of these pointers points +* to an array of "nsamp" doubles, being the normalised and sampled +* values for x1 and y1 in that order. +* scales +* Array holding the scaling factors used to produced the normalised +* values in the two columns of the table. Multiplying the normalised +* table values by the scale factor produces input or output axis +* values for the returned PolyMap +* ncoeff +* Pointer to an int in which to return the number of coefficients +* described by the returned array. +* racc +* Pointer to a double in which to return the achieved accuracy +* (which may be greater than "acc"). +* status +* Pointer to the inherited status variable. + +* Returned Value: +* A pointer to an array of doubles defining the polynomial in the +* form required by the PolyMap contructor. The number of coefficients +* is returned via "ncoeff". If the polynomial could not be found, +* then NULL is returned. The returned pointer should be freed using +* astFree when no longer needed. + +*/ + +/* Local Variables: */ + AstMinPackData data; + double *coeffs; + double *pc; + double *pr; + double *pxp1; + double *result; + double *work1; + double *work2; + double *work4; + double f1; + double f2; + double maxterm; + double term; + double tv; + int *work3; + int info; + int k; + int ncof; + int w1; + +/* Initialise returned value */ + result = NULL; + *ncoeff = 0; + *racc = 10*acc; + +/* Check inherited status */ + if( !astOK ) return result; + +/* Number of coefficients per poly. */ + ncof = order; + +/* Initialise the elements of the structure. */ + data.order = order; + data.nsamp = nsamp; + data.init_jac = 1; + data.xp1 = astMalloc( nsamp*order*sizeof( double ) ); + data.xp2 = NULL; + data.y[ 0 ] = table[ 1 ]; + data.y[ 1 ] = NULL; + +/* Work space to hold coefficients. */ + coeffs = astMalloc( ncof*sizeof( double ) ); + +/* Other work space. */ + work1 = astMalloc( nsamp*sizeof( double ) ); + work2 = astMalloc( ncof*nsamp*sizeof( double ) ); + work3 = astMalloc( ncof*sizeof( int ) ); + work4 = astMalloc( (5*ncof+nsamp)*sizeof( double ) ); + if( astOK ) { + +/* Find all the required powers of x1 and store them in the "xp1" + component of the data structure. The required initialisation is done + differently for different subclasses of PolyMap, so we need to wrap + it up in a virtual function. */ + astFitPoly1DInit( this, forward, table, &data, scales ); + +/* The initial guess at the coefficient values represents a unit + transformation in (normalised) tabulated (x,y) values. Using normalised + values means that we are, effectively, including a guess at the linear + scaling factor between input and output of the PolyMap (e.g. the PolyMap + may have inputs in mm and outputs in radians). */ + for( k = 0; k < ncof; k++ ) coeffs[ k ] = 0.0; + coeffs[ 1 ] = 1.0; + +/* Find the best coefficients */ + info = lmder1( MPFunc1D, &data, nsamp, ncof, coeffs, work1, work2, nsamp, + sqrt(DBL_EPSILON), work3, work4, (5*ncof+nsamp) ); + if( info == 0 ) astError( AST__MNPCK, "astPolyMap(PolyTran): Minpack error " + "detected (possible programming error).", status ); + +/* Return the achieved accuracy. The "work1" array holds the normalised Y + residuals at each tabulated point. */ + pr = work1; + tv = 0.0; + for( k = 0; k < nsamp; k++,pr++ ) tv += (*pr)*(*pr); + *racc = scales[ 1 ]*sqrt( tv/nsamp ); + +/* The best fitting polynomial coefficient found above relate to the + polynomial between the scaled positions stored in "table". These + scaled positions are related to PolyMap input/output axis values via + the scale factors supplied in "scales". Find the initial factor for the + current output. */ + f1 = scales[ 1 ]; + f2 = 1.0; + +/* Look at each coefficient. */ + pc = coeffs; + for( w1 = 0; w1 < order; w1++,pc++ ) { + +/* Get a pointer to the powers of X1 appropriate for the current coefficient, + at the first sample. */ + pxp1 = data.xp1 + w1; + +/* We find the contribution which this coefficient makes to the total + polynomial value. Find the maximum contribution made at any sample + points. */ + maxterm = 0.0; + for( k = 0; k < nsamp; k++ ) { + +/* Get the absolute value of the polynomial term that uses the current + coefficient. */ + term = fabs( ( *pc )*( *pxp1 ) ); + +/* Update the maximum term found at any sample. */ + if( term > maxterm ) maxterm = term; + +/* Increment the pointers to refer to the next sample. */ + pxp1 += order; + } + +/* If the maximum contribution made by this term is less than the + required accuracy, set the coefficient value to zero. */ + if( maxterm*f1 < acc ) { + *pc = 0.0; + +/* Scale the best fitting polynomial coefficient found above to take + account of the fact that the tabulated input and output positions in + "table" were are not actual PolyMap input and output axis values, but + are scaled by the factors stored in "scales". */ + } else { + *pc *= f1/f2; + } + + f2 *= scales[ 0 ]; + } + +/* Convert the array of coefficients into PolyMap form. */ + result = astMalloc( ncof*3*sizeof( double ) ); + + *ncoeff = 0; + pr = result; + pc = coeffs; + for( w1 = 0; w1 < order; w1++,pc++ ) { + if( *pc != 0.0 ) { + *(pr++) = *pc; + *(pr++) = 1; + *(pr++) = w1; + (*ncoeff)++; + } + } + +/* Truncate the returned array. */ + result = astRealloc( result, (*ncoeff)*3*sizeof( double ) ); + } + +/* Free resources. */ + coeffs = astFree( coeffs ); + data.xp1 = astFree( data.xp1 ); + work1 = astFree( work1 ); + work2 = astFree( work2 ); + work3 = astFree( work3 ); + work4 = astFree( work4 ); + +/* Return the coefficient array. */ + return result; + +} + +static void FitPoly1DInit( AstPolyMap *this, int forward, double **table, + AstMinPackData *data, double *scales, int *status ){ +/* +*+ +* Name: +* astFitPoly1DInit + +* Purpose: +* Perform initialisation needed for FitPoly1D + +* Type: +* Protected function. + +* Synopsis: +* #include "polymap.h" +* void astFitPoly1DInit( AstPolyMap *this, int forward, double **table, +* AstMinPackData *data, double *scales, +* int *status ) + +* Class Membership: +* PolyMap virtual function. + +* Description: +* This function performs initialisation needed for FitPoly1D. + +* Parameters: +* this +* Pointer to the PolyMap. +* forward +* Non-zero if the forward transformation of "this" is being +* replaced. Zero if the inverse transformation of "this" is being +* replaced. +* table +* Pointer to an array of 2 pointers. Each of these pointers points +* to an array of "nsamp" doubles, being the scaled and sampled values +* for x1 and y1 in that order. +* data +* Pointer to a structure holding information to pass the the +* service function invoked by the minimisation function. +* scales +* Array holding the scaling factors for the two columns of the table. +* Multiplying the table values by the scale factor produces PolyMap +* input or output axis values. +*- +*/ + +/* Local Variables; */ + double *px1; + double *pxp1; + double tv; + double x1; + int k; + int w1; + +/* Check the local error status. */ + if ( !astOK ) return; + +/* Get pointers to the supplied x1 values. */ + px1 = table[ 0 ]; + +/* Get pointers to the location for the next power of x1. */ + pxp1 = data->xp1; + +/* Loop round all samples. */ + for( k = 0; k < data->nsamp; k++ ) { + +/* Get the current x1 value. */ + x1 = *(px1++); + +/* Find all the required powers of x1 and store them in the "xp1" + component of the data structure. */ + tv = 1.0; + for( w1 = 0; w1 < data->order; w1++ ) { + *(pxp1++) = tv; + tv *= x1; + } + } +} + +static double *FitPoly2D( AstPolyMap *this, int forward, int nsamp, double acc, + int order, double **table, double scales[4], + int *ncoeff, double *racc, int *status ){ +/* +* Name: +* FitPoly2D + +* Purpose: +* Fit a (2-in,2-out) polynomial to a supplied set of data. + +* Type: +* Private function. + +* Synopsis: +* double *FitPoly2D( AstPolyMap *this, int forward, int nsamp, double acc, +* int order, double **table, double scales[4], +* int *ncoeff, double *racc, int *status ) + +* Description: +* This function fits a pair of least squares 2D polynomial surfaces +* to the positions in a supplied table, and returns the coefficients +* of a PolyMap to describe the fit. For the purposes of this function, +* the inputs to the fit is refered to as (x1,x2) and the output as +* (y1,y2). So the returned coefficients describe a PolyMap with forward +* transformations: +* +* y1 = P1( x1, x2 ) +* y2 = P2( x1, x2 ) +* +* P1 and P2 have the same maximum power on each input (specified by +* the "order" parameter). + +* Parameters: +* this +* Pointer to the PolyMap. +* forward +* Non-zero if the forward transformation of "this" is being +* replaced. Zero if the inverse transformation of "this" is being +* replaced. +* nsamp +* The number of (x1,x2,y1,y2) positions in the supplied table. +* acc +* The required accuracy, expressed as a geodesic distance within +* the polynomials output space (not the normalised tabular space). +* order +* The maximum power (plus one) of x1 or x2 within P1 and P2. So for +* instance, a value of 3 refers to a quadratic polynomial. Note, cross +* terms with total powers greater than or equal to "order" are not +* included in the fit. So the maximum number of terms in the fitted +* polynomial is order*(order+1)/2. +* table +* Pointer to an array of 4 pointers. Each of these pointers points +* to an array of "nsamp" doubles, being the normalised and sampled +* values for x1, x2, y1 or y2 in that order. +* scales +* Array holding the scaling factors used to produced the normalised +* values in the four columns of the table. Multiplying the normalised +* table values by the scale factor produces input or output axis +* values for the returned PolyMap +* ncoeff +* Pointer to an ant in which to return the number of coefficients +* described by the returned array. +* racc +* Pointer to a double in which to return the achieved accuracy +* (which may be greater than "acc"). +* status +* Pointer to the inherited status variable. + +* Returned Value: +* A pointer to an array of doubles defining the polynomial in the +* form required by the PolyMap contructor. The number of coefficients +* is returned via "ncoeff". If the polynomial could not be found with +* sufficient accuracy , then NULL is returned. The returned pointer +* should be freed using astFree when no longer needed. + +*/ + +/* Local Variables: */ + AstMinPackData data; + double *coeffs; + double *pc; + double *pr; + double *pxp1; + double *pxp2; + double *result; + double *work1; + double *work2; + double *work4; + double f1; + double f20; + double f2; + double f3; + double facc; + double maxterm; + double term; + double tv; + int *work3; + int info; + int iout; + int k; + int ncof; + int w12; + int w1; + int w2; + +/* Initialise returned value */ + result = NULL; + *ncoeff = 0; + *racc = 10*acc; + +/* Check inherited status */ + if( !astOK ) return result; + +/* Number of coefficients per poly. */ + ncof = order*( order + 1 )/2; + +/* Initialise the elements of the structure. */ + data.order = order; + data.nsamp = nsamp; + data.init_jac = 1; + data.xp1 = astMalloc( nsamp*order*sizeof( double ) ); + data.xp2 = astMalloc( nsamp*order*sizeof( double ) ); + data.y[ 0 ] = table[ 2 ]; + data.y[ 1 ] = table[ 3 ]; + +/* Work space to hold coefficients. */ + coeffs = astMalloc( 2*ncof*sizeof( double ) ); + +/* Other work space. */ + work1 = astMalloc( 2*nsamp*sizeof( double ) ); + work2 = astMalloc( 4*ncof*nsamp*sizeof( double ) ); + work3 = astMalloc( 2*ncof*sizeof( int ) ); + work4 = astMalloc( 2*(5*ncof+nsamp)*sizeof( double ) ); + if( astOK ) { + +/* Find all the required powers of x1 and x2 and store them in the "xp1" + and "xp2" components of the data structure. The required initialisation + is done differently for different subclasses of PolyMap, so we need to + wrap it up in a virtual function. */ + astFitPoly2DInit( this, forward, table, &data, scales ); + +/* The initial guess at the coefficient values represents a unit + transformation in (normalised) tabulated (x,y) values. Using normalised + values means that we are, effectively, including a guess at the linear + scaling factor between input and output of the PolyMap (e.g. the PolyMap + may have inputs in mm and outputs in radians). */ + for( k = 0; k < 2*ncof; k++ ) coeffs[ k ] = 0.0; + coeffs[ 1 ] = 1.0; + coeffs[ 5 ] = 1.0; + +/* Find the best coefficients */ + info = lmder1( MPFunc2D, &data, 2*nsamp, 2*ncof, coeffs, work1, work2, + 2*nsamp, sqrt(DBL_EPSILON), work3, work4, 2*(5*ncof+nsamp) ); + if( info == 0 ) astError( AST__MNPCK, "astPolyMap(PolyTran): Minpack error " + "detected (possible programming error).", status ); + +/* Return the achieved accuracy. */ + pr = work1; + tv = 0.0; + for( k = 0; k < 2*nsamp; k++,pr++ ) tv += (*pr)*(*pr); + facc = 1.0/(scales[2]*scales[2]) + 1.0/(scales[3]*scales[3]); + *racc = sqrt( tv/(2*nsamp*facc) ); + +/* Pointer to the first coefficient. */ + pc = coeffs; + +/* Look at coefficients for each output in turn. */ + for( iout = 0; iout < 2 && astOK; iout++ ) { + +/* The best fitting polynomial coefficient found above relate to the + polynomial between the scaled positions stored in "table". These + scaled positions are related to PolyMap input/output axis values via + the scale factors supplied in "scales". Find the initial factor for the + current output. */ + f1 = scales[ 2 + iout ]; + +/* Look at each coefficient for the current output. */ + f20 = 1.0; + for( w12 = 0; w12 < order; w12++ ) { + f3 = 1.0; + f2 = f20; + for( w2 = 0; w2 <= w12; w2++,pc++ ) { + w1 = w12 - w2; + +/* Get pointers to the powers of X1 and X2 appropriate for the current + coefficient, at the first sample. */ + pxp1 = data.xp1 + w1; + pxp2 = data.xp2 + w2; + +/* We find the contribution which this coefficient makes to the total + polynomial value. Find the maximum contribution made at any sample + points. */ + maxterm = 0.0; + for( k = 0; k < nsamp; k++ ) { + +/* Get the absolute value of the polynomial term that uses the current + coefficient. */ + term = fabs( ( *pc )*( *pxp1 )*( *pxp2 ) ); + +/* Update the maximum term found at any sample. */ + if( term > maxterm ) maxterm = term; + +/* Increment the pointers to refer to the next sample. */ + pxp1 += order; + pxp2 += order; + } + +/* If the maximum contribution made by this term is less than the + required accuracy, set the coefficient value to zero. */ + if( maxterm*f1 < acc ) { + *pc = 0.0; + +/* Scale the best fitting polynomial coefficient found above to take + account of the fact that the tabulated input and output positions in + "table" were are not actual PolyMap input and output axis values, but + are scaled by the factors stored in "scales". */ + } else { + *pc *= f1/( f2*f3 ); + } + + f2 /= scales[ 0 ]; + f3 *= scales[ 1 ]; + } + + f20 *= scales[ 0 ]; + } + } + +/* Convert the array of coefficients into PolyMap form. */ + result = astMalloc( 2*ncof*4*sizeof( double ) ); + + *ncoeff = 0; + pr = result; + pc = coeffs; + for( iout = 0; iout < 2 && astOK; iout++ ) { + for( w12 = 0; w12 < order; w12++ ) { + for( w2 = 0; w2 <= w12; w2++,pc++ ) { + w1 = w12 - w2; + if( *pc != 0.0 ) { + *(pr++) = *pc; + *(pr++) = iout + 1; + *(pr++) = w1; + *(pr++) = w2; + (*ncoeff)++; + } + } + } + } + +/* Truncate the returned array. */ + result = astRealloc( result, (*ncoeff)*4*sizeof( double ) ); + } + +/* Free resources. */ + coeffs = astFree( coeffs ); + data.xp1 = astFree( data.xp1 ); + data.xp2 = astFree( data.xp2 ); + work1 = astFree( work1 ); + work2 = astFree( work2 ); + work3 = astFree( work3 ); + work4 = astFree( work4 ); + +/* Return the coefficient array. */ + return result; + +} + +static void FitPoly2DInit( AstPolyMap *this, int forward, double **table, + AstMinPackData *data, double *scales, int *status ){ +/* +*+ +* Name: +* astFitPoly2DInit + +* Purpose: +* Perform initialisation needed for FitPoly2D + +* Type: +* Protected function. + +* Synopsis: +* #include "polymap.h" +* void astFitPoly2DInit( AstPolyMap *this, int forward, double **table, +* AstMinPackData *data, double *scales, +* int *status ) + +* Class Membership: +* PolyMap virtual function. + +* Description: +* This function performs initialisation needed for FitPoly2D. + +* Parameters: +* this +* Pointer to the PolyMap. +* forward +* Non-zero if the forward transformation of "this" is being +* replaced. Zero if the inverse transformation of "this" is being +* replaced. +* table +* Pointer to an array of 4 pointers. Each of these pointers points +* to an array of "nsamp" doubles, being the scaled and sampled values +* for x1, x2, y1 or y2 in that order. +* data +* Pointer to a structure holding information to pass the the +* service function invoked by the minimisation function. +* scales +* Array holding the scaling factors for the four columns of the table. +* Multiplying the table values by the scale factor produces PolyMap +* input or output axis values. +*- +*/ + +/* Local Variables; */ + double *px1; + double *px2; + double *pxp1; + double *pxp2; + double tv; + double x1; + double x2; + int k; + int w1; + int w2; + +/* Check the local error status. */ + if ( !astOK ) return; + +/* Get pointers to the supplied x1 and x2 values. */ + px1 = table[ 0 ]; + px2 = table[ 1 ]; + +/* Get pointers to the location for the next power of x1 and x2. */ + pxp1 = data->xp1; + pxp2 = data->xp2; + +/* Loop round all samples. */ + for( k = 0; k < data->nsamp; k++ ) { + +/* Get the current x1 and x2 values. */ + x1 = *(px1++); + x2 = *(px2++); + +/* Find all the required powers of x1 and store them in the "xp1" + component of the data structure. */ + tv = 1.0; + for( w1 = 0; w1 < data->order; w1++ ) { + *(pxp1++) = tv; + tv *= x1; + } + +/* Find all the required powers of x2 and store them in the "xp2" + comonent of the data structure. */ + tv = 1.0; + for( w2 = 0; w2 < data->order; w2++ ) { + *(pxp2++) = tv; + tv *= x2; + } + } +} + +static void FreeArrays( AstPolyMap *this, int forward, int *status ) { +/* +* Name: +* FreeArrays + +* Purpose: +* Free the dynamic arrays contained within a PolyMap + +* Type: +* Private function. + +* Synopsis: +* void FreeArrays( AstPolyMap *this, int forward, int *status ) + +* Description: +* This function frees all the dynamic arrays allocated as part of a +* PolyMap. + +* Parameters: +* this +* Pointer to the PolyMap. +* forward +* If non-zero, the arrays for the forward transformation are freed. +* Otherwise, the arrays for the inverse transformation are freed. +* 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: */ + int nin; /* Number of inputs */ + int nout; /* Number of outputs */ + int i; /* Loop count */ + int j; /* Loop count */ + +/* Get the number of inputs and outputs of the uninverted Mapping. */ + nin = ( (AstMapping *) this )->nin; + nout = ( (AstMapping *) this )->nout; + +/* Free the dynamic arrays for the forward transformation. */ + if( forward != astGetInvert( this ) ) { + + if( this->coeff_f ) { + for( i = 0; i < nout; i++ ) { + this->coeff_f[ i ] = astFree( this->coeff_f[ i ] ); + } + this->coeff_f = astFree( this->coeff_f ); + } + + if( this->power_f ) { + for( i = 0; i < nout; i++ ) { + if( this->ncoeff_f && this->power_f[ i ] ) { + for( j = 0; j < this->ncoeff_f[ i ]; j++ ) { + this->power_f[ i ][ j ] = astFree( this->power_f[ i ][ j ] ); + } + } + this->power_f[ i ] = astFree( this->power_f[ i ] ); + } + this->power_f = astFree( this->power_f ); + } + + this->ncoeff_f = astFree( this->ncoeff_f ); + this->mxpow_f = astFree( this->mxpow_f ); + +/* Free the dynamic arrays for the inverse transformation. */ + } else { + + if( this->coeff_i ) { + for( i = 0; i < nin; i++ ) { + this->coeff_i[ i ] = astFree( this->coeff_i[ i ] ); + } + this->coeff_i = astFree( this->coeff_i ); + } + + if(this->power_i ) { + for( i = 0; i < nin; i++ ) { + if( this->ncoeff_i && this->power_i[ i ] ) { + for( j = 0; j < this->ncoeff_i[ i ]; j++ ) { + this->power_i[ i ][ j ] = astFree( this->power_i[ i ][ j ] ); + } + } + this->power_i[ i ] = astFree( this->power_i[ i ] ); + } + this->power_i = astFree( this->power_i ); + } + + this->ncoeff_i = astFree( this->ncoeff_i ); + this->mxpow_i = astFree( this->mxpow_i ); + } +} + +static const char *GetAttrib( AstObject *this_object, const char *attrib, int *status ) { +/* +* Name: +* GetAttrib + +* Purpose: +* Get the value of a specified attribute for a PolyMap. + +* Type: +* Private function. + +* Synopsis: +* #include "polymap.h" +* const char *GetAttrib( AstObject *this, const char *attrib, int *status ) + +* Class Membership: +* PolyMap member function (over-rides the protected astGetAttrib +* method inherited from the Mapping class). + +* Description: +* This function returns a pointer to the value of a specified +* attribute for a PolyMap, formatted as a character string. + +* Parameters: +* this +* Pointer to the PolyMap. +* attrib +* Pointer to a null-terminated string containing the name of +* the attribute whose value is required. This name should be in +* lower case, with all white space removed. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* - Pointer to a null-terminated string containing the attribute +* value. + +* Notes: +* - The returned string pointer may point at memory allocated +* within the PolyMap, or at static memory. The contents of the +* string may be over-written or the pointer may become invalid +* following a further invocation of the same function or any +* modification of the PolyMap. A copy of the string should +* therefore be made if necessary. +* - 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: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstPolyMap *this; /* Pointer to the PolyMap structure */ + const char *result; /* Pointer value to return */ + double dval; /* Floating point attribute value */ + int ival; /* Integer attribute value */ + +/* Initialise. */ + result = NULL; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Get a pointer to the thread specific global data structure. */ + astGET_GLOBALS(this_object); + +/* Obtain a pointer to the PolyMap structure. */ + this = (AstPolyMap *) this_object; + +/* Compare "attrib" with each recognised attribute name in turn, + obtaining the value of the required attribute. If necessary, write + the value into "getattrib_buff" as a null-terminated string in an appropriate + format. Set "result" to point at the result string. */ + +/* IterInverse. */ +/* ------------ */ + if ( !strcmp( attrib, "iterinverse" ) ) { + ival = astGetIterInverse( this ); + if ( astOK ) { + (void) sprintf( getattrib_buff, "%d", ival ); + result = getattrib_buff; + } + +/* NiterInverse. */ +/* ------------- */ + } else if ( !strcmp( attrib, "niterinverse" ) ) { + ival = astGetNiterInverse( this ); + if ( astOK ) { + (void) sprintf( getattrib_buff, "%d", ival ); + result = getattrib_buff; + } + +/* TolInverse. */ +/* ----------- */ + } else if ( !strcmp( attrib, "tolinverse" ) ) { + dval = astGetTolInverse( this ); + if ( astOK ) { + (void) sprintf( getattrib_buff, "%.*g", AST__DBL_DIG, dval ); + result = getattrib_buff; + } + +/* If the attribute name was not recognised, pass it on to the parent + method for further interpretation. */ + } else { + result = (*parent_getattrib)( this_object, attrib, status ); + } + +/* Return the result. */ + return result; + +} + +static AstPolyMap **GetJacobian( AstPolyMap *this, int *status ){ +/* +* Name: +* GetJacobian + +* Purpose: +* Get a description of a Jacobian matrix for the original +* fwd transformation of a PolyMap. + +* Type: +* Private function. + +* Synopsis: +* AstPolyMap **GetJacobian( AstPolyMap *this, int *status ) + +* Description: +* This function returns a set of PolyMaps which define the Jacobian +* matrix of the original forward transformation of the supplied PolyMap +* (i.e. the Negated attribute is assumed to be zero). +* +* The Jacobian matrix has "nout" rows and "nin" columns, where "nin" +* and "nout" are the number of inputs and outputs of the original PolyMap. +* Row "i", column "j" of the matrix holds the rate of change of the +* i'th PolyMap output with respect to the j'th PolyMap input. +* +* Since the values in the Jacobian matrix vary across the input space +* of the PolyMap, the matrix is returned in the form of a set of new +* PolyMaps which generate the elements of the Jacobian for any given +* position in the input space. The "nout" values in a single column of +* the Jacobian matrix are generated by the "nout" outputs from a single +* new PolyMap. The whole matrix is described by "nin" PolyMaps. +* +* The returned PolyMaps are cached in the supplied PolyMap object in +* order to speed up subsequent calls to this function. + +* Parameters: +* this +* The PolyMap for which the Jacbian is required. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* A pointer to an array of "nin" PolyMap pointers, where "nin" is the number +* of inputs for the sipplied PolyMap. The returned array should not be changed +* in any way, and the PolyMaps should not be freed (they will be freed when +* the supplied PolyMap is deleted). + +*/ + +/* Local Variables: */ + double *coeffs; + double *pc; + int icof; + int icol; + int iin; + int irow; + int ncof; + int ncof_row; + int ncof_total; + int nin; + int nout; + int power; + +/* Check inherited status */ + if( !astOK ) return NULL; + +/* Ensure there is a Jacobian stored in the PolyMap. */ + if( !this->jacobian ) { + +/* Get the number of inputs and outputs of the original forward + transformation. */ + if( astGetInvert( this ) ){ + nout = astGetNin( this ); + nin = astGetNout( this ); + } else { + nin = astGetNin( this ); + nout = astGetNout( this ); + } + +/* Allocate memory to hold pointers to the PolyMaps used to describe the + Jacobian matrix. */ + this->jacobian = astCalloc( nin, sizeof(AstPolyMap *) ); + +/* Find the total number of coefficients used to describe the supplied + PolyMap (forward transformation) and allocate work space to hold the + coefficients for a single new PolyMap forward transformation. */ + ncof = 0; + for( irow = 0; irow < nout; irow++ ) { + ncof += this->ncoeff_f[ irow ]; + } + coeffs = astMalloc( ncof*( 2 + nin )*sizeof( double ) ); + +/* Check pointers can be used safely. */ + if( astOK ) { + +/* The Jacobian matrix has "nout" rows and "nin" columns. The "nout" values + in a single column of the Jacobian matrix corresponds to the "nout" outputs + from a single PolyMap. The whole matrix is described by "nin" PolyMaps. + Loop over each column of the Matrix, creating the corresponding PolyMap + for each. */ + for( icol = 0; icol < nin; icol++ ) { + +/* Initialise the total number of coefficients used to describe the + element of the PolyMap. */ + ncof_total = 0; + +/* Loop over each row of the Jacobian matrix (i.e. each PolyMap output). */ + pc = coeffs; + for( irow = 0; irow < nout; irow++ ) { + +/* Loop over each coefficient used in the polynomial that generates the + current PolyMap output. */ + ncof_row = this->ncoeff_f[ irow ]; + for( icof = 0; icof < ncof_row; icof++ ) { + +/* Get the power of input "icol" associated with the current coefficient. */ + power = (int)( this->power_f[ irow ][ icof ][ icol ] + 0.5 ); + +/* We can skip the coefficient if the power is zero. */ + if( power > 0 ) { + ncof_total++; + +/* Store the coefficient value, modified so that it describes a + polynomial that has been differentiated with respect to input "icol". */ + *(pc++) = this->coeff_f[ irow ][ icof ]*power; + +/* Store the output PolyMap to which the coeff relates. */ + *(pc++) = irow + 1; + +/* Store the powers of the inputs associated with the coeff. These are + the same as the original powers, except that the power of "icol" + (the input with respect to which the output has been differentiated) + is reduced by one. */ + for( iin = 0; iin < nin; iin++ ) { + if( iin != icol ) { + *(pc++) = this->power_f[ irow ][ icof ][ iin ]; + } else { + *(pc++) = this->power_f[ irow ][ icof ][ iin ] - 1; + } + } + } + } + } + +/* Create the PolyMap and store a pointer to it in the jacobian array in + the supplied PolyMap. */ + (this->jacobian)[ icol ] = astPolyMap( nin, nout, ncof_total, coeffs, + 0, NULL, " ", status ); + } + } + +/* Free resources */ + coeffs = astFree( coeffs ); + } + +/* Return the Jacobian. */ + return this->jacobian; +} + +static int GetObjSize( AstObject *this_object, int *status ) { +/* +* Name: +* GetObjSize + +* Purpose: +* Return the in-memory size of an Object. + +* Type: +* Private function. + +* Synopsis: +* #include "polymap.h" +* int GetObjSize( AstObject *this, int *status ) + +* Class Membership: +* PolyMap member function (over-rides the astGetObjSize protected +* method inherited from the parent class). + +* Description: +* This function returns the in-memory size of the supplied PolyMap, +* in bytes. + +* Parameters: +* this +* Pointer to the PolyMap. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* The Object size, in bytes. + +* 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: */ + AstPolyMap *this; + int ic; + int nc; + int result; + +/* Initialise. */ + result = 0; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Obtain a pointers to the PolyMap structure. */ + this = (AstPolyMap *) this_object; + +/* Invoke the GetObjSize method inherited from the parent class, and then + add on any components of the class structure defined by this class + which are stored in dynamically allocated memory. */ + result = (*parent_getobjsize)( this_object, status ); + + if( this->jacobian ) { + nc = astGetNin( this ); + for( ic = 0; ic < nc; ic++ ) { + result += astGetObjSize( (this->jacobian)[ ic ] ); + } + result += sizeof( AstPolyMap * )*nc; + } + +/* If an error occurred, clear the result value. */ + if ( !astOK ) result = 0; + +/* Return the result, */ + return result; +} + +static int GetTranForward( AstMapping *this, int *status ) { +/* +* +* Name: +* GetTranForward + +* Purpose: +* Determine if a PolyMap defines a forward coordinate transformation. + +* Type: +* Private function. + +* Synopsis: +* #include "mapping.h" +* int GetTranForward( AstMapping *this, int *status ) + +* Class Membership: +* PolyMap member function (over-rides the astGetTranForward method +* inherited from the Mapping class). + +* Description: +* This function returns a value indicating if the PolyMap is able +* to perform a forward coordinate transformation. + +* Parameters: +* this +* Pointer to the PolyMap. +* 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: */ + AstPolyMap *map; /* Pointer to PolyMap to be queried */ + int result; /* The returned value */ + +/* Check the global error status. */ + if ( !astOK ) return 0; + +/* Obtain a pointer to the PolyMap. */ + map = (AstPolyMap *) this; + +/* First deal with cases where the PolyMap has not been inverted. */ + if( !astGetInvert( this ) ) { + +/* The PolyMap has a defined forward transformation if one or more + coefficients values were supplied for the original forward + transformation. It is not possible to replace the original forward + transformation with an iterative algorithm. */ + result = map->ncoeff_f ? 1 : 0; + +/* Now deal with cases where the PolyMap has been inverted. */ + } else { + +/* The PolyMap has a defined forward transformation if one or more + coefficients values were supplied for the original inverse + transformation, or if the original inverse transformation is being + approximated using an iterative algorithm. */ + result = ( map->ncoeff_i || astGetIterInverse( map ) ) ? 1 : 0; + } + +/* Return the result. */ + return result; +} + +static int GetTranInverse( AstMapping *this, int *status ) { +/* +* +* Name: +* GetTranInverse + +* Purpose: +* Determine if a PolyMap defines an inverse coordinate transformation. + +* Type: +* Private function. + +* Synopsis: +* #include "matrixmap.h" +* int GetTranInverse( AstMapping *this, int *status ) + +* Class Membership: +* PolyMap member function (over-rides the astGetTranInverse method +* inherited from the Mapping class). + +* Description: +* This function returns a value indicating if the PolyMap is able +* to perform an inverse coordinate transformation. + +* Parameters: +* this +* Pointer to the PolyMap. +* 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: */ + AstPolyMap *map; /* Pointer to PolyMap to be queried */ + int result; /* The returned value */ + +/* Check the global error status. */ + if ( !astOK ) return 0; + +/* Obtain a pointer to the PolyMap. */ + map = (AstPolyMap *) this; + +/* First deal with cases where the PolyMap has not been inverted. */ + if( !astGetInvert( this ) ) { + +/* The PolyMap has a defined inverse transformation if one or more + coefficients values were supplied for the original inverse + transformation, or if the original inverse transformation is being + approximated using an iterative algorithm. */ + result = ( map->ncoeff_i || astGetIterInverse( map ) ) ? 1 : 0; + +/* Now deal with cases where the PolyMap has been inverted. */ + } else { + +/* The PolyMap has a defined inverse transformation if one or more + coefficients values were supplied for the original forward + transformation. It is not possible to replace the original forward + transformation with an iterative algorithm. */ + result = map->ncoeff_f ? 1 : 0; + } + +/* Return the result. */ + return result; +} + +void astInitPolyMapVtab_( AstPolyMapVtab *vtab, const char *name, int *status ) { +/* +*+ +* Name: +* astInitPolyMapVtab + +* Purpose: +* Initialise a virtual function table for a PolyMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "polymap.h" +* void astInitPolyMapVtab( AstPolyMapVtab *vtab, const char *name ) + +* Class Membership: +* PolyMap vtab initialiser. + +* Description: +* This function initialises the component of a virtual function +* table which is used by the PolyMap 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 astIsAPolyMap) 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->PolyPowers = PolyPowers; + vtab->FitPoly1DInit = FitPoly1DInit; + vtab->FitPoly2DInit = FitPoly2DInit; + vtab->PolyTran = PolyTran; + vtab->PolyCoeffs = PolyCoeffs; + + vtab->ClearIterInverse = ClearIterInverse; + vtab->GetIterInverse = GetIterInverse; + vtab->SetIterInverse = SetIterInverse; + vtab->TestIterInverse = TestIterInverse; + + vtab->ClearNiterInverse = ClearNiterInverse; + vtab->GetNiterInverse = GetNiterInverse; + vtab->SetNiterInverse = SetNiterInverse; + vtab->TestNiterInverse = TestNiterInverse; + + vtab->ClearTolInverse = ClearTolInverse; + vtab->GetTolInverse = GetTolInverse; + vtab->SetTolInverse = SetTolInverse; + vtab->TestTolInverse = TestTolInverse; + +/* 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_getobjsize = object->GetObjSize; + object->GetObjSize = GetObjSize; + +#if defined(THREAD_SAFE) + parent_managelock = object->ManageLock; + object->ManageLock = ManageLock; +#endif + + parent_clearattrib = object->ClearAttrib; + object->ClearAttrib = ClearAttrib; + parent_getattrib = object->GetAttrib; + object->GetAttrib = GetAttrib; + parent_setattrib = object->SetAttrib; + object->SetAttrib = SetAttrib; + parent_testattrib = object->TestAttrib; + object->TestAttrib = TestAttrib; + + parent_transform = mapping->Transform; + mapping->Transform = Transform; + mapping->GetTranForward = GetTranForward; + mapping->GetTranInverse = GetTranInverse; + +/* Store replacement pointers for methods which will be over-ridden by + new member functions implemented here. */ + object->Equal = Equal; + mapping->MapMerge = MapMerge; + +/* Declare the destructor and copy constructor. */ + astSetDelete( (AstObjectVtab *) vtab, Delete ); + astSetCopy( (AstObjectVtab *) vtab, Copy ); + +/* Declare the class dump function. */ + astSetDump( vtab, Dump, "PolyMap", "Polynomial 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 void IterInverse( AstPolyMap *this, AstPointSet *out, AstPointSet *result, + int *status ){ +/* +* Name: +* IterInverse + +* Purpose: +* Use an iterative method to evaluate the original inverse transformation +* of a PolyMap at a set of (original) output positions. + +* Type: +* Private function. + +* Synopsis: +* void IterInverse( AstPolyMap *this, AstPointSet *out, AstPointSet *result, +* int *status ) + +* Description: +* This function transforms a set of PolyMap positions using the original +* inverse transformation of the PolyMap (i.e. the Negated attribute +* is assumed to be zero). An iterative Newton-Raphson method is used +* which only required the original forward transformation of the PolyMap +* to be defined. + +* Parameters: +* this +* The PolyMap. +* out +* A PointSet holding the positions that are to be transformed using +* the original inverse transformation. These corresponds to +* outputs of the original (i.e. uninverted) PolyMap +* result +* A PointSet into which the transformed positions are to be stored. +* These corresponds to inputs of the original (i.e. uninverted) PolyMap + +* status +* Pointer to the inherited status variable. + +*/ + +/* Local Variables: */ + AstMapping *lintrunc; + AstPointSet *work; + AstPointSet **ps_jac; + AstPolyMap **jacob; + double *vec; + double *pb; + double **ptr_work; + double ***ptr_jac; + double *mat; + double **ptr_out; + double **ptr_in; + double *pa; + double det; + double maxerr; + double vlensq; + double xlensq; + double xx; + int *flags; + int *iw; + int fwd; + int icol; + int icoord; + int ipoint; + int irow; + int iter; + int maxiter; + int nconv; + int ncoord; + int npoint; + int sing; + +/* Check inherited status */ + if( !astOK ) return; + +/* Check the PolyMap has equal numbers of inputs and outputs. */ + ncoord = astGetNin( this ); + if( ncoord != astGetNout( this ) ) { + astError( AST__INTER, "astTransform(%s): Supplied %s has unequal numbers" + " of inputs and outputs and therefore an iterative inverse " + "cannot be used (internal AST Programming errpr).", status, + astGetClass(this), astGetClass(this) ); + } + +/* Get information about the Jacobian matrix for the forward polynomial + transformation. This matrix is a ncoord X ncoord matrix, in which + element (row=I,col=J) is the rate of change of output coord I with + respect to input coord J, of the supplied PolyMap's forward transformation. + The numerical values of the matrix vary depending on where it is + evaluated within the input space of the PolyMap. For this reason, the + "jacob" variable holds a vector of "ncoord" PolyMaps. The outputs of + each of these PolyMaps corresponds to a single column in the Jacobian + matrix. */ + jacob = GetJacobian( this, status ); + +/* Get the number of points to be transformed. */ + npoint = astGetNpoint( out ); + +/* Get another PointSet to hold intermediate results. */ + work = astPointSet( npoint, ncoord, " ", status ); + +/* See if the PolyMap has been inverted.*/ + fwd = !astGetInvert( this ); + +/* Get pointers to the data arrays for all PointSets. Note, here "in" and + "out" refer to inputs and outputs of the PolyMap (i.e. the forward + transformation). These are respectively *outputs* and *inputs* of the + inverse transformation. */ + ptr_in = astGetPoints( result ); /* Returned input positions */ + ptr_out = astGetPoints( out ); /* Supplied output positions */ + ptr_work = astGetPoints( work ); /* Work space */ + +/* Allocate an array of PointSets to hold the elements of the Jacobian + matrix. */ + ptr_jac = astMalloc( sizeof( double ** )*ncoord ); + ps_jac = astCalloc( ncoord, sizeof( AstPointSet * ) ); + if( astOK ) { + for( icoord = 0; icoord < ncoord; icoord++ ) { + ps_jac[ icoord ] = astPointSet( npoint, ncoord, " ", status ); + ptr_jac[ icoord ] = astGetPoints( ps_jac[ icoord ] ); + } + } + +/* Allocate an array to hold flags indicating if each position has + converged. Initialise it to hold zero at every element. */ + flags = astCalloc( npoint, sizeof( int ) ); + +/* Allocate memory to hold the Jacobian matrix at a single point. */ + mat = astMalloc( sizeof( double )*ncoord*ncoord ); + +/* Allocate memory to hold the offset vector. */ + vec = astMalloc( sizeof( double )*ncoord ); + +/* Allocate memory to hold work space for palDmat. */ + iw = astMalloc( sizeof( int )*ncoord ); + +/* Check pointers can be used safely. */ + if( astOK ) { + +/* Store the initial guess at the required input positions. These are + determined by transforming the supplied output positions using the + inverse of a linear truncation of the PolyMap's forward + transformation. */ + lintrunc = LinearGuess( this, status ); + (void) astTransform( lintrunc, out, 0, result ); + lintrunc = astAnnul( lintrunc ); + +/* Get the maximum number of iterations to perform. */ + maxiter = astGetNiterInverse( this ); + +/* Get the target relative error for the returned input axis values, and + square it. */ + maxerr = astGetTolInverse( this ); + maxerr *= maxerr; + +/* Initialise the number of positions which have reached the required + accuracy. */ + nconv = 0; + +/* Loop round doing iterations of a Newton-Raphson algorithm, until + all points have achieved the required relative error, or the + maximum number of iterations have been performed. */ + for( iter = 0; iter < maxiter && nconv < npoint && astOK; iter++ ) { + +/* Use the original forward transformation of the supplied PolyMap to + transform the current guesses at the required input positions into + the corresponding output positions. Store the results in the "work" + PointSet. */ + (void) astTransform( this, result, fwd, work ); + +/* Modify the work PointSet so that it holds the offsets from the output + positions produced by the current input position guesses, and the + required output positions. */ + for( icoord = 0; icoord < ncoord; icoord++ ) { + pa = ptr_out[ icoord ]; + pb = ptr_work[ icoord ]; + for( ipoint = 0; ipoint< npoint; ipoint++,pa++,pb++ ) { + if( *pa != AST__BAD && *pb != AST__BAD ){ + *pb = *pa - *pb; + } else { + *pb = AST__BAD; + } + } + } + +/* Evaluate the elements of the Jacobian matrix at the current input + position guesses. */ + for( icoord = 0; icoord < ncoord; icoord++ ) { + (void) astTransform( jacob[ icoord ], result, 1, ps_jac[ icoord ] ); + } + +/* For each position, we now invert the matrix equation + + Dy = Jacobian.Dx + + to find a guess at the vector (dx) holding the offsets from the + current input positions guesses to their required values. Loop over all + points. */ + for( ipoint = 0; ipoint < npoint; ipoint++ ) { + +/* Do not change positions that have already converged. */ + if( !flags[ ipoint ] ) { + +/* Get the numerical values for the elements of the Jacobian matrix at + the current point. */ + pa = mat; + for( irow = 0; irow < ncoord; irow++ ) { + for( icol = 0; icol < ncoord; icol++ ) { + *(pa++) = ptr_jac[ icol ][ irow ][ ipoint ]; + } + +/* Store the offset from the current output position to the required + output position. */ + vec[ irow ] = ptr_work[ irow ][ ipoint ]; + } + +/* Find the corresponding offset from the current input position to the required + input position. */ + palDmat( ncoord, mat, vec, &det, &sing, iw ); + +/* If the matrix was singular, the input position cannot be evaluated so + store a bad value for it and indicate it has converged. */ + if( sing ) { + for( icoord = 0; icoord < ncoord; icoord++ ) { + ptr_in[ icoord ][ ipoint ] = AST__BAD; + } + flags[ ipoint ] = 1; + nconv++; + +/* Otherwise, update the input position guess. */ + } else { + vlensq = 0.0; + xlensq = 0.0; + pa = vec; + for( icoord = 0; icoord < ncoord; icoord++,pa++ ) { + xx = ptr_in[ icoord ][ ipoint ] + (*pa); + ptr_in[ icoord ][ ipoint ] = xx; + xlensq += xx*xx; + vlensq += (*pa)*(*pa); + } + +/* Check for convergence. */ + if( vlensq <= maxerr*xlensq ) { + flags[ ipoint ] = 1; + nconv++; + } + } + } + } + } + } + +/* Free resources. */ + vec = astFree( vec ); + iw = astFree( iw ); + mat = astFree( mat ); + flags = astFree( flags ); + work = astAnnul( work ); + + if( ps_jac ) { + for( icoord = 0; icoord < ncoord; icoord++ ) { + ps_jac[ icoord ] = astAnnul( ps_jac[ icoord ] ); + } + ps_jac = astFree( ps_jac ); + } + + ptr_jac = astFree( ptr_jac ); + +} + +static AstMapping *LinearGuess( AstPolyMap *this, int *status ){ +/* +* Name: +* LinearTruncation + +* Purpose: +* Get a Mapping representing a linear approximation of a PolyMap + +* Type: +* Private function. + +* Synopsis: +* AstMapping *LinearGuess( AstPolyMap *this, int *status ) + +* Description: +* This function returns a linear Mapping approximating the original +* forward transformation of the supplied PolyMap (i.e. the Invert +* flag is assume to be zero). The linear and constant terms in the +* supplied PolyMap are used for all outputs that have such terms. Any +* other outputs are assumed to be equal to the corresponding inputs. +* If the forward transformation of the PolyMap is defined, then it is +* used to determine the returned Mapping. Otherwise, the inverse +* transformation of the PolyMap is used. The forward transformation of +* the returned Mapping always represents the forward transformation of +* the original (i.e. uninverted) PolyMap. + +* Parameters: +* this +* Pointer to the PolyMap. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* The returned linear Mapping, or NULL if an error occurs. + +*/ + +/* Local Variables: */ + AstMapping *result; /* Pointer to returned Mapping */ + AstMapping *mm; /* MatrixMap representing linear terms */ + AstMapping *sm; /* ShiftMap representing offset terms */ + double **coeff; /* Pointer to coefficient value arrays */ + double *matrix; /* Linear scaling matrix */ + double *vector; /* Offset vector */ + int ***power; /* Pointer to coefficient power arrays */ + int *ncoeff; /* Pointer to no. of coefficients */ + int flag; /* Indicates nature of current coeff */ + int ico; /* Coefficient index */ + int iin; /* Index of transformation input */ + int inverse; /* PolyMap inverse transformation selected? */ + int iout; /* Index of transformation output */ + int jin; /* Transformation input for current coeff */ + int nin; /* No. of i/ps for selected transformation */ + int nout; /* No. of o/ps for selected transformation */ + int ok; /* Current output is defined? */ + int pow; /* Power for current coeff and input */ + +/* Initialise returned value */ + result = NULL; + +/* Check inherited status */ + if( !astOK ) return result; + +/* If the required Mapping has already been determined, return a pointer + to it. */ + if( this->lintrunc ) return astClone( this->lintrunc ); + +/* Get pointers to the arrays holding the required coefficient + values and powers. Use the forward transformation of the original + (uninverted) PolyMap if it is defined and the inverse transformation + otherwise. */ + if ( this->ncoeff_f ){ + ncoeff = this->ncoeff_f; + coeff = this->coeff_f; + power = this->power_f; + nin = astGetNin( this ); + nout = astGetNout( this ); + inverse = 0; + + } else { + ncoeff = this->ncoeff_i; + coeff = this->coeff_i; + power = this->power_i; + nout = astGetNin( this ); + nin = astGetNout( this ); + inverse = 1; + } + +/* Allocate memory to hold a matrix holding the linear terms of the + PolyMap transformation selected above. Initialise it to hold zeros. */ + matrix = astCalloc( nin*nout, sizeof( double * ) ); + +/* Allocate memory to hold a vector holding the offset terms of the + PolyMap transformation selected above. Initialise it to hold zeros. */ + vector = astCalloc( nout, sizeof( double * ) ); + if( astOK ){ + +/* Loop round all the outputs of the transformation selected above. */ + + for( iout = 0; iout < nout; iout++ ){ + +/* Loop round all the coefficients for the current transformation output. */ + for( ico = 0; ico < ncoeff[iout]; ico++ ){ + +/* Initialise a tristate flag that indicates if the current coefficient is + constant (0), linear (1) or something else (2). */ + flag = 0; + +/* Look for the first input axis for which the current coeff uses a + non-zero power. If the power is not one, the current coeff is neither + linear nor constant, and so we can pass on to the next coeff. If any + subsequent non-zero power is found, the coeff is not linear. */ + jin = 0; + for( iin = 0; iin < nin; iin++ ){ + pow = power[iout][ico][iin]; + if( pow ) { + if( pow != 1 ) { + flag = 2; + break; + } else if( flag == 1 ) { + flag = 2; + break; + } else { + flag = 1; + jin = iin; + } + } + } + +/* If the coeff is a constant term, insert it into the offset vector. */ + if( flag == 0 ) { + vector[ iout ] = coeff[ iout ][ ico ]; + +/* If the coeff is a linear term, insert it into the matrix. */ + } else if( flag == 1 ) { + matrix[ jin + iout*nin ] = coeff[ iout ][ ico ]; + } + } + } + +/* If any PolyMap output has no linear terms, the matrix will not be + invertable. So insert a 1.0 value on the diagonal of such rows. */ + for( iout = 0; iout < nout; iout++ ){ + ok = 0; + for( iin = 0; iin < nin; iin++ ) { + if( matrix[ iin + iout*nin ] != 0.0 ) { + ok = 1; + break; + } + } + if( ! ok ) matrix[ iout + iout*nin ] = 1.0; + } + +/* Create a MatrixMap to represent the matrix. */ + mm = (AstMapping *) astMatrixMap( nin, nout, 0, matrix, " ", status ); + +/* If the MatrixMap cannot be inverted, use a unit map instead. */ + if( !astGetTranInverse( mm ) ) { + mm = astAnnul( mm ); + mm = (AstMapping *) astUnitMap( nin, " ", status ); + } + +/* Create a ShiftMap to represent the offset. */ + sm = (AstMapping *) astShiftMap( nout, vector, " ", status ); + +/* Combine them in series into a single compound Mapping. */ + result = (AstMapping *) astCmpMap( mm, sm, 1, " ", status ); + +/* If this CmpMap represents the inverse transformation of the original + (i.e. uninverted) PolyMap, invert it so that it represents the forward + transformation. */ + if( inverse ) astInvert( result ); + +/* Free resources */ + mm = astAnnul( mm ); + sm = astAnnul( sm ); + +/* Store the Mapping in the PolyMap structure so we do not need to + recalculate it every time it is needed. */ + if( astOK ) this->lintrunc = astClone( result ); + } + +/* Free resources */ + matrix = astFree( matrix ); + vector = astFree( vector ); + +/* Return the required Mapping. */ + return result; + +} + +static void LMFunc1D( const double *p, double *hx, int m, int n, void *adata ){ +/* +* Name: +* LMFunc1D + +* Purpose: +* Evaluate a test 1D polynomal. + +* Type: +* Private function. + +* Synopsis: +* void LMFunc1D( const double *p, double *hx, int m, int n, void *adata ) + +* Description: +* This function finds the residuals implied by a supplied set of +* candidate polynomial coefficients. Each residual is a candidate +* polynomial evaluated at one of the sample points, minus the +* supplied target value for the polynomial at that test point. +* +* The minimisation process minimises the sum of the squared residuals. + +* Parameters: +* p +* An array of "m" coefficients for the candidate polynomial. The +* coefficients are ordered C0, C1, C2, etc. +* hx +* An array in which to return the "n" residuals. The residual at +* sample "k" is returned in element (k). +* m +* The length of the "p" array. This should be equal to order. +* n +* The length of the "hx" array. This should be equal to nsamp. +* adata +* Pointer to a structure holding the sample positions and values, +* and other information. + +*/ + +/* Local Variables: */ + AstMinPackData *data; + double *px1; + double *py; + const double *vp; + double *vr; + double res; + int k; + int w1; + +/* Get a pointer to the data structure. */ + data = (AstMinPackData *) adata; + +/* Initialise a pointer to the current returned residual value. */ + vr = hx; + +/* Initialise a pointer to the sampled Y values for the polynomial output. */ + py = data->y[ 0 ]; + +/* Initialise a pointer to the powers of the input X values at the curent + (i.e. first) sample. */ + px1 = data->xp1; + +/* Loop over the index of the sample to which this residual refers. */ + for( k = 0; k < data->nsamp; k++ ) { + +/* Initialise a pointer to the first coefficient (the constant term) for the + polynomial output coordinate. */ + vp = p; + +/* Initialise this residual to hold the sampled Y value. Increment the + pointer to the next sampled value for the current polynomial output. */ + res = -( *(py++) ); + +/* Loop over the coefficients. */ + for( w1 = 0; w1 < data->order; w1++ ) { + +/* Increment the residual by the value of the current term Cw1*(x1^w1). + Increment the pointer to the next coefficient (C). Also increment the + pointer to the next higher power of X1. */ + res += ( *(vp++) )*( *(px1++) ); + } + +/* Store the complete residual in the returned array, and increment the + pointer to the next residual. */ + *(vr++) = res; + } +} + +static void LMFunc2D( const double *p, double *hx, int m, int n, void *adata ){ +/* +* Name: +* LMFunc2D + +* Purpose: +* Evaluate a test 2D polynomal. + +* Type: +* Private function. + +* Synopsis: +* void LMFunc2D( const double *p, double *hx, int m, int n, void *adata ) + +* Description: +* This function finds the residuals implied by a supplied set of +* candidate polynomial coefficients. Each residual is a candidate +* polynomial (either P1 or P2) evaluated at one of the sample points +* (x1,x2), minus the supplied target value for the polynomial at that +* test point. +* +* The minimisation process minimises the sum of the squared residuals. + +* Parameters: +* p +* An array of "m" coefficients for the candidate polynomial. All the +* coefficients for polynomial P1 come first, followed by those for P2. +* Within each each polynomial, the coefficients are order C00, C10, +* C01, C20, C11, C02, C30, C21, C12, C03, etc. So the coefficient +* of (x1^j*x2^k) (=Cjk) for polynomial Pi is stored in element +* [k + (j + k)*(j + k + 1)/2 + i*order*(order+1)/2] of the "p" array. +* hx +* An array in which to return the "n" residuals. The residual at +* sample "k" for polynomial "i" is returned in element (k + nsamp*i). +* m +* The length of the "p" array. This should be equal to order*(order+1). +* n +* The length of the "hx" array. This should be equal to 2*nsamp. +* adata +* Pointer to a structure holding the sample positions and values, +* and other information. + +*/ + +/* Local Variables: */ + AstMinPackData *data; + const double *vp0; + const double *vp; + double *py; + double *px1; + double *px10; + double *px20; + double *vr; + double res; + double *px2; + int iout; + int k; + int w12; + int w2; + +/* Get a pointer to the data structure. */ + data = (AstMinPackData *) adata; + +/* Initialise a pointer to the current returned residual value. */ + vr = hx; + +/* Initialise a pointer to the first coefficient (the constant term) for the + current (i.e. first) polynomial output coordinate. */ + vp0 = p; + +/* Loop over each polynomial output coordinate. */ + for( iout = 0; iout < 2; iout++ ) { + +/* Initialise a pointer to the sampled Y values for the first polynomial + output. */ + py = data->y[ iout ]; + +/* Initialise pointers to the powers of the input X values at the curent + (i.e. first) sample. */ + px10 = data->xp1; + px20 = data->xp2; + +/* Loop over the index of the sample to which this residual refers. */ + for( k = 0; k < data->nsamp; k++ ) { + +/* Reset the pointer to the first coefficient (the constant term) + for the current polynomial output coordinate. */ + vp = vp0; + +/* Initialise this residual to hold the sampled Y value. Increment the + pointer to the next sampled value for the current polynomial output. */ + res = -( *(py++) ); + +/* The w12 value is the sum of the powers of X1 and X2. So w12=0 + corresponds to the constant term in the polynomial, and (e.g.) w12=6 + corresponds to all the terms for which the sum of the powerss (w1+w2) + is 6. Loop over all possible w12 values. */ + for( w12 = 0; w12 < data->order; w12++ ) { + +/* The next coeff refers to (x1^w12)*(x2^0). Get pointers to the values + holding x1^w12 and x2^0. */ + px1 = px10++; + px2 = px20; + +/* Loop over powers of x2. The corresponding power of x1 is "w12-w2", but + is not explicitly needed here. So x1 moves down from w12 to zero, as + w2 moves up from zero to w12. */ + for( w2 = 0; w2 <= w12; w2++ ) { + +/* Increment the residual by the value of the current term Cw1w2*(x1^w1)*(x2^w2). + Increment the pointer tio the next coefficient (C). Also decrement the + pointer to the next lower power of X1, and increment the pointer to the next + higher power of X2. */ + res += ( *(vp++) )*( *(px1--) )*( *(px2++) ); + } + } + +/* Move on to the x2 powers for the next sample. Don't need to do this + for x1 since px10 is incremented within the above loop. */ + px20 += data->order; + +/* Store the complete residual in the returned array, and increment the + pointer to the next residual. */ + *(vr++) = res; + } + +/* Get a pointer to the first coefficient (the constant term) for the + next polynomial output coordinate. */ + vp0 += data->order*( 1 + data->order )/2; + } +} + +static void LMJacob1D( const double *p, double *jac, int m, int n, void *adata ){ +/* +* Name: +* LMJacob1D + +* Purpose: +* Evaluate the Jacobian matrix of a test 1D polynomal. + +* Type: +* Private function. + +* Synopsis: +* void LMJacob1D( const double *p, double *jac, int m, int n, void *adata ) + +* Description: +* This function finds the Jacobian matrix that describes the rate of +* change of every residual with respect to every polynomial coefficient. +* Each residual is a candidate polynomial evaluated at one of the sample +* points minus the supplied target value for the polynomial at that test +* point. +* +* For a polynomial the Jacobian matrix is constant (i.e. does not +* depend on the values of the polynomial coefficients). So we only +* evaluate it on the first call. + +* Parameters: +* p +* An array of "m" coefficients for the candidate polynomial. +* jac +* An array in which to return the "m*n" elements of the Jacobian +* matrix. The rate of change of residual "r" with respect to +* coefficient "c" is returned in element "r + c*n". The residual +* at sample "k" of polynomial Pi has an "r" index of (k + nsamp*i). +* The coefficient of (x1^j) for polynomial Pi has a "c" index +* of j. +* m +* The length of the "p" array. This should be equal to order. +* n +* The number of residuals. This should be equal to nsamp. +* adata +* Pointer to a structure holding the sample positions and values, +* and other information. + +*/ + +/* Local Variables: */ + AstMinPackData *data; + int k; + int w1; + +/* Get a pointer to the data structure. */ + data = (AstMinPackData *) adata; + +/* The Jacobian of the residuals with respect to the polynomial + coefficients is constant (i.e. does not depend on the values of the + polynomial coefficients). So we only need to calculate it once. If + this is the first call, calculate the Jacobian and return it in "jac". + otherwise, just return immediately retaining the supplied "jac" values + (which will be the values returned by the previous call to this + function). */ + if( data->init_jac ) { + data->init_jac = 0; + +/* Loop over all residuals. */ + for( k = 0; k < n; k++ ) { + +/* Loop over all parameters (i.e. polynomial coefficients). */ + for( w1 = 0; w1 < m; w1++ ) { + +/* Store the Jacobian. */ + jac[ k + w1*n ] = (data->xp1)[ w1 + k*data->order ]; + } + } + } +} + +static void LMJacob2D( const double *p, double *jac, int m, int n, void *adata ){ +/* +* Name: +* LMJacob2D + +* Purpose: +* Evaluate the Jacobian matrix of a test 2D polynomal. + +* Type: +* Private function. + +* Synopsis: +* void LMJacob2D( const double *p, double *jac, int m, int n, void *adata ) + +* Description: +* This function finds the Jacobian matrix that describes the rate of +* change of every residual with respect to every polynomial coefficient. +* Each residual is a candidate polynomial (either P1 or P2) evaluated +* at one of the sample points (x1,x2), minus the supplied target value +* for the polynomial at that test point. +* +* For a polynomial the Jacobian matrix is constant (i.e. does not +* depend on the values of the polynomial coefficients). So we only +* evaluate it on the first call. + +* Parameters: +* p +* An array of "m" coefficients for the candidate polynomial. All the +* coefficients for polynomial P1 come first, followed by those for P2. +* Within each each polynomial, the coefficients are order C00, C10, +* C01, C20, C11, C02, C30, C21, C12, C03, etc. So the coefficient +* of (x1^j*x2^k) (=Cjk) for polynomial Pi is stored in element +* [k + (j + k)*(j + k + 1)/2 + i*order*(order+1)/2] of the "p" array. +* jac +* An array in which to return the "m*n" elements of the Jacobian +* matrix. The rate of change of residual "r" with respect to +* coefficient "c" is returned in element "r + c*n". The residual +* at sample "k" of polynomial Pi has an "r" index of (k + nsamp*i). +* The coefficient of (x1^j*x2^k) for polynomial Pi has a "c" index +* of [k + (j + k)*(j + k + 1)/2 + i*order*(order+1)/2]. +* m +* The length of the "p" array. This should be equal to order*(order+1). +* n +* The number of residuals. This should be equal to 2*nsamp. +* adata +* Pointer to a structure holdin gthe sample positions and values, +* and other information. + +*/ + +/* Local Variables: */ + AstMinPackData *data; + int iout; + int k; + int ncof; + int vp; + int vr; + int w1; + int w12; + int w2; + +/* Get a pointer to the data structure. */ + data = (AstMinPackData *) adata; + +/* The Jacobian of the residuals with respect to the polynomial + coefficients is constant (i.e. does not depend on the values of the + polynomial coefficients). So we only need to calculate it once. If + this is the first call, calculate the Jacobian and return it in "jac". + otherwise, just return immediately retaining the supplied "jac" values + (which will be the values returned by the previous call to this + function). */ + if( data->init_jac ) { + data->init_jac = 0; + +/* Store the number of coefficients in one polynomial. */ + ncof = data->order*( 1 + data->order )/2; + +/* Loop over all residuals. */ + for( vr = 0; vr < n; vr++ ) { + +/* Calculate the polynomial output index, and sample index, that creates + the current residual. */ + iout = vr/data->nsamp; + k = vr - iout*data->nsamp; + +/* Loop over all parameters (i.e. polynomial coefficients). */ + for( vp = 0; vp < m; vp++ ) { + +/* If this coefficient is not used in the creation of the current + polynomial output value, then the Jacobian value is zero. */ + if( vp/ncof != iout ) { + jac[ vr + vp*n ] = 0.0; + +/* Otherwise, get the powers of the two polynomial inputs, to which + the current coefficient relates. */ + } else { + w12 = ( -1.0 + sqrt( 1.0 + 8.0*(vp - iout*ncof) ) )/2.0; + w2 = vp - iout*ncof - w12*( w12 + 1 )/2; + w1 = w12 - w2; + +/* Store the Jacobian. */ + jac[ vr + vp*n ] = (data->xp1)[ w1 + k*data->order ]* + (data->xp2)[ w2 + k*data->order ]; + } + } + } + } +} + +#if defined(THREAD_SAFE) +static int ManageLock( AstObject *this_object, int mode, int extra, + AstObject **fail, int *status ) { +/* +* Name: +* ManageLock + +* Purpose: +* Manage the thread lock on an Object. + +* Type: +* Private function. + +* Synopsis: +* #include "object.h" +* AstObject *ManageLock( AstObject *this, int mode, int extra, +* AstObject **fail, int *status ) + +* Class Membership: +* PolyMap member function (over-rides the astManageLock protected +* method inherited from the parent class). + +* Description: +* This function manages the thread lock on the supplied Object. The +* lock can be locked, unlocked or checked by this function as +* deteremined by parameter "mode". See astLock for details of the way +* these locks are used. + +* Parameters: +* this +* Pointer to the Object. +* mode +* An integer flag indicating what the function should do: +* +* AST__LOCK: Lock the Object for exclusive use by the calling +* thread. The "extra" value indicates what should be done if the +* Object is already locked (wait or report an error - see astLock). +* +* AST__UNLOCK: Unlock the Object for use by other threads. +* +* AST__CHECKLOCK: Check that the object is locked for use by the +* calling thread (report an error if not). +* extra +* Extra mode-specific information. +* fail +* If a non-zero function value is returned, a pointer to the +* Object that caused the failure is returned at "*fail". This may +* be "this" or it may be an Object contained within "this". Note, +* the Object's reference count is not incremented, and so the +* returned pointer should not be annulled. A NULL pointer is +* returned if this function returns a value of zero. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* A local status value: +* 0 - Success +* 1 - Could not lock or unlock the object because it was already +* locked by another thread. +* 2 - Failed to lock a POSIX mutex +* 3 - Failed to unlock a POSIX mutex +* 4 - Bad "mode" value supplied. + +* Notes: +* - This function attempts to execute even if an error has already +* occurred. +*/ + +/* Local Variables: */ + AstPolyMap *this; + int ic; + int result; + int nc; + +/* Initialise */ + result = 0; + +/* Check the supplied pointer is not NULL. */ + if( !this_object ) return result; + +/* Obtain a pointer to the PolyMap structure. */ + this = (AstPolyMap *) this_object; + +/* Invoke the ManageLock method inherited from the parent class. */ + if( !result ) result = (*parent_managelock)( this_object, mode, extra, + fail, status ); + +/* Invoke the astManageLock method on any Objects contained within + the supplied Object. */ + if( this->jacobian ) { + nc = astGetNin( this ); + for( ic = 0; ic < nc && result; ic++ ) { + result = astManageLock( (this->jacobian)[ ic ], mode, + extra, fail ); + } + } + + return result; + +} +#endif + +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 PolyMap. + +* 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: +* PolyMap 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 PolyMap 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 PolyMap 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 PolyMap which is to be merged with +* its neighbours. This should be a cloned copy of the PolyMap +* pointer contained in the array element "(*map_list)[where]" +* (see below). This pointer will not be annulled, and the +* PolyMap it identifies will not be modified by this function. +* where +* Index in the "*map_list" array (below) at which the pointer +* to the nominated PolyMap 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: */ + AstPolyMap *pmap0; /* Nominated PolyMap */ + AstPolyMap *pmap1; /* Neighbouring PolyMap */ + int i; /* Index of neighbour */ + int nin; /* Number of input coordinates for nominated PolyMap */ + int nout; /* Number of output coordinates for nominated PolyMap */ + int ok; /* Are PolyMaps equivalent? */ + int oldinv0; /* Original Invert value in pmap0 */ + int oldinv1; /* Original Invert value in pmap1 */ + int result; /* Result value to return */ + +/* Initialise. */ + result = -1; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Save a pointer to the nominated PolyMap. */ + pmap0 = (AstPolyMap *) ( *map_list )[ where ]; + +/* The only simplification which can currently be performed is to merge a PolyMap + with its own inverse. This can only be done in series. Obviously, + there are potentially other simplications which could be performed, but + time does not currently allow these to be coded. */ + if( series ) { + +/* Temporarily set the Invert flag of the nominated PolyMap to the + required value, first saving the original value so that it can be + re-instated later. */ + oldinv0 = astGetInvert( pmap0 ); + astSetInvert( pmap0, ( *invert_list )[ where ] ); + +/* Get the number of inputs and outputs to the nominated PolyMap. */ + nin = astGetNin( pmap0 ); + nout = astGetNout( pmap0 ); + +/* Check each neighbour. */ + for( i = where - 1; i <= where + 1; i += 2 ) { + +/* Continue with the next pass if the neighbour does not exist. */ + if( i < 0 || i >= *nmap ) continue; + +/* Continue with the next pass if this neighbour is not a PermMap. */ + if( ! astIsAPolyMap( ( *map_list )[ i ] ) ) continue; + +/* Get a pointer to it. */ + pmap1 = (AstPolyMap *) ( *map_list )[ i ]; + +/* Check it is used in the opposite direction to the nominated PolyMap. */ + if( ( *invert_list )[ i ] == ( *invert_list )[ where ] ) continue; + +/* We use the astEqual method to check that the two PolyMaps are equal. + But at the moment they may not be equal because they may have + different Invert flags. Therefore, temporarily set the invert flag + of the neighbour so that it is the same as the nominated PolyMap, + first saving the original value so that it can be re-instated later. + Note, we have already checked that the two PolyMaps are used in opposite + directions within the CmpMap. */ + oldinv1 = astGetInvert( pmap1 ); + astSetInvert( pmap1, ( *invert_list )[ where ] ); + +/* Use astEqual to check that the coefficients etc are equal in the two + PolyMaps. */ + ok = astEqual( pmap0, pmap1 ); + +/* Re-instate the original value of the Invert flag in the neighbour. */ + astSetInvert( pmap1, oldinv1 ); + +/* Pass on to the next neighbour if the current neighbour differs from + the nominated PolyMap. */ + if( !ok ) continue; + +/* If we get this far, then the nominated PolyMap and the current + neighbour cancel each other out, so replace each by a UnitMap. */ + pmap0 = astAnnul( pmap0 ); + pmap1 = astAnnul( pmap1 ); + if( i < where ) { + ( *map_list )[ where ] = (AstMapping *) astUnitMap( nout, "", status ); + ( *map_list )[ i ] = (AstMapping *) astUnitMap( nout, "", status ); + ( *invert_list )[ where ] = 0; + ( *invert_list )[ i ] = 0; + result = i; + } else { + ( *map_list )[ where ] = (AstMapping *) astUnitMap( nin, "", status ); + ( *map_list )[ i ] = (AstMapping *) astUnitMap( nin, "", status ); + ( *invert_list )[ where ] = 0; + ( *invert_list )[ i ] = 0; + result = where; + } + +/* Leave the loop. */ + break; + } + +/* If the nominated PolyMap was not replaced by a UnitMap, then + re-instate its original value for the Invert flag. */ + if( pmap0 ) astSetInvert( pmap0, oldinv0 ); + } + +/* Return the result. */ + return result; +} + +static int MPFunc1D( void *p, int m, int n, const double *x, double *fvec, + double *fjac, int ldfjac, int iflag ) { +/* +* Name: +* MPFunc1D + +* Purpose: +* Evaluate a test 1D polynomal or its Jacobian. + +* Type: +* Private function. + +* Synopsis: +* int MPFunc1D( void *p, int m, int n, const double *x, double *fvec, +* double *fjac, int ldfjac, int iflag ) + +* Description: +* This function finds the residuals or Jacobian implied by a supplied +* set of candidate polynomial coefficients. Each residual is a candidate +* polynomial evaluated at one of the sample points, minus the +* supplied target value for the polynomial at that test point. +* +* The minimisation process minimises the sum of the squared residuals. + +* Parameters: +* p +* A pointer to data needed to evaluate the required residuals or +* Jacobians. +* m +* The number of residuals. +* n +* The number of coefficients. +* x +* An array of "n" coefficients for the candidate polynomial. +* fvec +* An array in which to return the "m" residuals. +* fjac +* An array in which to return the "mxn" Jacobian values. +* ldflac +* Unused (should always be equal to "m"). +* iflag +* If 1 calculate the residuals, otherwise calculate the Jacobians. + +*/ + +/* If required, calculate the function values at X, and return this + vector in fvec. Do not alter fjac. */ + if( iflag == 1 ) { + LMFunc1D( x, fvec, n, m, p ); + +/* Otherwise, calculate the Jacobian values at X, and return this + matrix in fjac. Do not alter fvec. */ + } else { + LMJacob1D( x, fjac, n, m, p ); + } + + return 0; +} + +static int MPFunc2D( void *p, int m, int n, const double *x, double *fvec, + double *fjac, int ldfjac, int iflag ) { +/* +* Name: +* MPFunc1D + +* Purpose: +* Evaluate a test 2D polynomal or its Jacobian. + +* Type: +* Private function. + +* Synopsis: +* int MPFunc2D( void *p, int m, int n, const double *x, double *fvec, +* double *fjac, int ldfjac, int iflag ) + +* Description: +* This function finds the residuals or Jacobian implied by a supplied +* set of candidate polynomial coefficients. Each residual is a candidate +* polynomial evaluated at one of the sample points, minus the +* supplied target value for the polynomial at that test point. +* +* The minimisation process minimises the sum of the squared residuals. + +* Parameters: +* p +* A pointer to data needed to evaluate the required residuals or +* Jacobians. +* m +* The number of residuals. +* n +* The number of coefficients. +* x +* An array of "n" coefficients for the candidate polynomial. +* fvec +* An array in which to return the "m" residuals. +* fjac +* An array in which to return the "mxn" Jacobian values. +* ldflac +* Unused (should always be equal to "m"). +* iflag +* If 1 calculate the residuals, otherwise calculate the Jacobians. + +*/ + +/* If required, calculate the function values at X, and return this + vector in fvec. Do not alter fjac. */ + if( iflag == 1 ) { + LMFunc2D( x, fvec, n, m, p ); + +/* Otherwise, calculate the Jacobian values at X, and return this + matrix in fjac. Do not alter fvec. */ + } else { + LMJacob2D( x, fjac, n, m, p ); + + } + + return 0; +} + +static void PolyCoeffs( AstPolyMap *this, int forward, int nel, double *coeffs, + int *ncoeff, int *status ){ +/* +*++ +* Name: +c astPolyCoeffs +f AST_POLYCOEFFS + +* Purpose: +* Retrieve the coefficient values used by a PolyMap. + +* Type: +* Public function. + +* Synopsis: +c #include "polymap.h" +c void astPolyCoeffs( AstPolyMap *this, int forward, int nel, double *coeffs, +c int *ncoeff ) +f CALL AST_POLYCOEFFS( THIS, FORWARD, NEL, COEFFS, NCOEFF, STATUS ) + +* Class Membership: +* PolyMap method. + +* Description: +* This function returns the coefficient values used by either the +* forward or inverse transformation of a PolyMap, in the same form +* that they are supplied to the PolyMap constructor. +* +* Usually, you should call this method first with +c "nel" +f NEL +* set to zero to determine the number of coefficients used by the +* PolyMap. This allows you to allocate an array of the correct size to +* hold all coefficient data. You should then call this method a +* second time to get the coefficient data. + +* Parameters: +c this +f THIS = INTEGER (Given) +* Pointer to the original Mapping. +c forward +f FORWARD = LOGICAL (Given) +c If non-zero, +f If .TRUE., +* the coefficients of the forward PolyMap transformation are +* returned. Otherwise the inverse transformation coefficients are +* returned. +c nel +f NEL = INTEGER (Given) +* The length of the supplied +c "coeffs" +f COEFFS +* array. It should be at least "ncoeff*( nin + 2 )" if +c "foward" is non-zero, +f FORWARD is .TRUE., +* and "ncoeff*( nout + 2 )" otherwise, where "ncoeff" is the +* number of coefficients to be returned. If a value of zero +* is supplied, no coefficient values are returned, but the +* number of coefficients used by the transformation is still +* returned in +c "ncoeff". +f NCOEFF. +c coeffs +f COEFFS( NEL ) = DOUBLE PRECISION (Returned) +* An array in which to return the coefficients used by the +* requested transformation of the PolyMap. Ignored if +c "nel" is zero. +f NEL is zero. +* The coefficient data is returned in the form in which it is +* supplied to the PolyMap constructor. That is, each group of +* "2 + nin" or "2 + nout" adjacent elements describe a single +* coefficient of the forward or inverse transformation. See the +* PolyMap constructor documentation for further details. +* +* If the supplied array is too short to hold all the coefficients, +* trailing coefficients are excluded. If the supplied array is +* longer than needed to hold all the coefficients, trailing +* elements are filled with zeros. +c ncoeff +f NCOEFF = INTEGER (Returned) +* The number of coefficients used by the requested transformation. +* A value of zero is returned if the transformation does not +* have any defining polynomials. A value is returned for this argument +* even if +c "nel" is zero. +f NEL is zero. +f STATUS = INTEGER (Given and Returned) +f The global status. + +*-- +*/ + +/* Local Variables: */ + double **coeff; + int ***power; + int *nco; + int *pp; + int iax; + int icoeff; + int iel; + int ipoly; + int nax; + int npoly; + +/* Initialise */ + *ncoeff = 0; + +/* Check the inherited status. */ + if ( !astOK ) return; + +/* Fill any supplied array with zeros. */ + if( nel ) memset( coeffs, 0, nel*sizeof( *coeffs ) ); + +/* Get the values to use, taking account of whether the PolyMap has been + inverted or not. */ + if( forward != astGetInvert( this ) ){ + nco = this->ncoeff_f; + power = this->power_f; + coeff = this->coeff_f; + npoly = astGetNout( this ); + nax = astGetNin( this ); + } else { + nco = this->ncoeff_i; + power = this->power_i; + coeff = this->coeff_i; + npoly = astGetNin( this ); + nax = astGetNout( this ); + } + +/* Notheg to do if there are no coeffs. */ + if( nco && power && coeff ) { + +/* Initialise index of next value to store in returned array. */ + iel = 0; + +/* Loop round each 1D polynomial. */ + for( ipoly = 0; ipoly < npoly; ipoly++ ){ + +/* Loop round coefficients. */ + for( icoeff = 0; icoeff < nco[ ipoly ]; icoeff++ ) { + +/* Store the coefficient value in the next element of the returned array, + if the array is not already full. Increment the pointer to the next + element. */ + if( iel < nel ) coeffs[ iel++ ] = coeff[ipoly][icoeff]; + +/* Next store the integer index of the PolyMap input or output which uses + the coefficient within its defining polynomial (the first axis has + index 1). */ + if( iel < nel ) coeffs[ iel++ ] = ipoly + 1; + +/* The remaining elements of the group give the integer powers to use + with each input or output coordinate value. */ + pp = power[ipoly][icoeff]; + for( iax = 0; iax < nax; iax++,pp++ ) { + if( iel < nel ) coeffs[ iel++ ] = *pp; + } + } + +/* Increment the total number of coefficients used by the transformation. */ + *ncoeff += nco[ ipoly ]; + } + } +} + +static void PolyPowers( AstPolyMap *this, double **work, int ncoord, + const int *mxpow, double **ptr, int point, + int fwd, int *status ){ +/* +*+ +* Name: +* astPolyPowers + +* Purpose: +* Find the required powers of the input axis values. + +* Type: +* Protected function. + +* Synopsis: +* #include "polymap.h" +* void astPolyPowers( AstPolyMap *this, double **work, int ncoord, +* const int *mxpow, double **ptr, int point, +* int fwd ) + +* Class Membership: +* PolyMap virtual function. + +* Description: +* This function is used by astTransform to calculate the powers of +* the axis values for a single input position. In the case of +* sub-classes, the powers may not be simply powers of the supplied +* axis values but may be more complex quantities such as a Chebyshev +* polynomial of the required degree evaluated at the input axis values. + +* Parameters: +* this +* Pointer to the PolyMap. +* work +* An array of "ncoord" pointers, each pointing to an array of +* length "max(2,mxpow)". The required values are placed in this +* array on exit. +* ncoord +* The number of axes. +* mxpow +* Pointer to an array holding the maximum power required of each +* axis value. Should have "ncoord" elements. +* ptr +* An array of "ncoord" pointers, each pointing to an array holding +* the axis values. Each of these arrays of axis values must have +* at least "point+1" elements. +* point +* The zero based index of the point within "ptr" that holds the +* axis values to be exponentiated. +* fwd +* Do the supplied coefficients define the foward transformation of +* the PolyMap? +*- +*/ + +/* Local Variables; */ + double *pwork; + double x; + int coord; + int ip; + +/* Check the local error status. */ + if ( !astOK ) return; + +/* For the base PolyMap class, this method simply raises each input axis + value to the required power. Loop over all input axes. */ + for( coord = 0; coord < ncoord; coord++ ) { + +/* Get a pointer to the array in which the powers of the current axis + value are to be returned. */ + pwork = work[ coord ]; + +/* Anything to the power zero is 1.0. */ + pwork[ 0 ] = 1.0; + +/* Get the input axis value. If it is bad, store bad values for all + remaining powers. */ + x = ptr[ coord ][ point ]; + if( x == AST__BAD ) { + for( ip = 1; ip <= mxpow[ coord ]; ip++ ) pwork[ ip ] = AST__BAD; + +/* Otherwise, form and store the required powers of the input axis value. */ + } else { + for( ip = 1; ip <= mxpow[ coord ]; ip++ ) { + pwork[ ip ] = pwork[ ip - 1 ]*x; + } + } + } +} + +static AstPolyMap *PolyTran( AstPolyMap *this, int forward, double acc, + double maxacc, int maxorder, const double *lbnd, + const double *ubnd, int *status ){ +/* +*++ +* Name: +c astPolyTran +f AST_POLYTRAN + +* Purpose: +* Fit a PolyMap inverse or forward transformation. + +* Type: +* Public function. + +* Synopsis: +c #include "polymap.h" +c AstPolyMap *astPolyTran( AstPolyMap *this, int forward, double acc, +c double maxacc, int maxorder, const double *lbnd, +c const double *ubnd ) +f RESULT = AST_POLYTRAN( THIS, FORWARD, ACC, MAXACC, MAXORDER, LBND, +f UBND, STATUS ) + +* Class Membership: +* PolyMap method. + +* Description: +* This function creates a new PolyMap which is a copy of the supplied +* PolyMap, in which a specified transformation (forward or inverse) +* has been replaced by a new polynomial transformation. The +* coefficients of the new transformation are estimated by sampling +* the other transformation and performing a least squares polynomial +* fit in the opposite direction to the sampled positions and values. +* +* This method can only be used on (1-input,1-output) or (2-input,2-output) +* PolyMaps. +* +* The transformation to create is specified by the +c "forward" parameter. +f FORWARD parameter. +* In what follows "X" refers to the inputs of the PolyMap, and "Y" to +* the outputs of the PolyMap. The forward transformation transforms +* input values (X) into output values (Y), and the inverse transformation +* transforms output values (Y) into input values (X). Within a PolyMap, +* each transformation is represented by an independent set of +* polynomials, P_f or P_i: Y=P_f(X) for the forward transformation and +* X=P_i(Y) for the inverse transformation. +* +c The "forward" +f The FORWARD +* parameter specifies the transformation to be replaced. If it is +c non-zero, +f is .TRUE., +* a new forward transformation is created +* by first finding the input values (X) using the inverse transformation +* (which must be available) at a regular grid of points (Y) covering a +* rectangular region of the PolyMap's output space. The coefficients of +* the required forward polynomial, Y=P_f(X), are chosen in order to +* minimise the sum of the squared residuals between the sampled values +* of Y and P_f(X). +* +c If "forward" is zero (probably the most likely case), +f If FORWARD is .FALSE. (probably the most likely case), +* a new inverse transformation is created by +* first finding the output values (Y) using the forward transformation +* (which must be available) at a regular grid of points (X) covering a +* rectangular region of the PolyMap's input space. The coefficients of +* the required inverse polynomial, X=P_i(Y), are chosen in order to +* minimise the sum of the squared residuals between the sampled values +* of X and P_i(Y). +* +* This fitting process is performed repeatedly with increasing +* polynomial orders (starting with linear) until the target +* accuracy is achieved, or a specified maximum order is reached. If +* the target accuracy cannot be achieved even with this maximum-order +* polynomial, the best fitting maximum-order polynomial is returned so +* long as its accuracy is better than +c "maxacc". +f MAXACC. +* If it is not, a NULL pointer is returned but no error is reported. + +* Parameters: +c this +f THIS = INTEGER (Given) +* Pointer to the original Mapping. +c forward +f FORWARD = LOGICAL (Given) +c If non-zero, +f If .TRUE., +* the forward PolyMap transformation is replaced. Otherwise the +* inverse transformation is replaced. +c acc +f ACC = DOUBLE (Given) +* The target accuracy, expressed as a geodesic distance within +* the PolyMap's input space (if +c "forward" is zero) or output space (if "forward" is non-zero). +f FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.). +c maxacc +f MAXACC = DOUBLE (Given) +* The maximum allowed accuracy for an acceptable polynomial, +* expressed as a geodesic distance within the PolyMap's input +* space (if +c "forward" is zero) or output space (if "forward" is non-zero). +f FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.). +c maxorder +f MAXORDER = INTEGER (Given) +* The maximum allowed polynomial order. This is one more than the +* maximum power of either input axis. So for instance, a value of +* 3 refers to a quadratic polynomial. Note, cross terms with total +* powers greater than or equal to +c maxorder +f MAXORDER +* are not inlcuded in the fit. So the maximum number of terms in +* each of the fitted polynomials is +c maxorder*(maxorder+1)/2. +f MAXORDER*(MAXORDER+1)/2. +c lbnd +f LBND( * ) = DOUBLE PRECISION (Given) +c Pointer to an +f An +* array holding the lower bounds of a rectangular region within +* the PolyMap's input space (if +c "forward" is zero) or output space (if "forward" is non-zero). +f FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.). +* The new polynomial will be evaluated over this rectangle. The +* length of this array should equal the value of the PolyMap's Nin +* or Nout attribute, depending on +c "forward". +f FORWARD. +c ubnd +f UBND( * ) = DOUBLE PRECISION (Given) +c Pointer to an +f An +* array holding the upper bounds of a rectangular region within +* the PolyMap's input space (if +c "forward" is zero) or output space (if "forward" is non-zero). +f FORWARD is .FALSE.) or output space (if FORWARD is .TRUE.). +* The new polynomial will be evaluated over this rectangle. The +* length of this array should equal the value of the PolyMap's Nin +* or Nout attribute, depending on +c "forward". +f FORWARD. +f STATUS = INTEGER (Given and Returned) +f The global status. + +* Returned Value: +c astPolyTran() +f AST_POLYTRAN = INTEGER +* A pointer to the new PolyMap. +c A NULL pointer +f AST__NULL +* will be returned if the fit fails to achieve the accuracy +* specified by +c "maxacc", +f MAXACC, +* but no error will be reported. + +* Applicability: +* PolyMap +* All PolyMaps have this method. +* ChebyMap +c The ChebyMap implementation of this method allows +c NULL pointers to be supplied for "lbnd" and/or "ubnd", +c in which case the corresponding bounds supplied when the ChebyMap +c was created are used. +* The returned PolyMap will be a ChebyMap, and the new transformation +* will be defined as a weighted sum of Chebyshev functions of the +* first kind. + +* Notes: +* - The IterInverse attribute is always cleared in the returned PolyMap. +* This means that the returned PolyMap will always use the new fit by +* default, rather than the iterative inverse, regardless of the setting +* of IterInverse in the supplied PolyMap. +* - This function can only be used on 1D or 2D PolyMaps which have +* the same number of inputs and outputs. +* - 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. +*-- +*/ + +/* Local Variables: */ + AstPolyMap *result; + int ok; + +/* Initialise. */ + result = NULL; + +/* Check the inherited status. */ + if ( !astOK ) return result; + +/* Take a copy of the supplied PolyMap. */ + result = astCopy( this ); + +/* Replace the required transformation. */ + ok = ReplaceTransformation( result, forward, acc, maxacc, maxorder, lbnd, + ubnd, status ); + +/* If an error occurred, or the fit was not good enough, annul the returned + PolyMap. */ + if ( !ok || !astOK ) { + result = astAnnul( result ); + +/* Otherwise, ensure the new fit is used in preference to the iterative + inverse by clearing the IterInverse attribute. */ + } else { + astClearIterInverse( result ); + } + +/* Return the result. */ + return result; +} + +static int ReplaceTransformation( AstPolyMap *this, int forward, double acc, + double maxacc, int maxorder, const double *lbnd, + const double *ubnd, int *status ){ +/* +* Name: +* ReplaceTransformation + +* Purpose: +* Create a new inverse or forward transformation for a PolyMap. + +* Type: +* Private function. + +* Synopsis: +* int ReplaceTransformation( AstPolyMap *this, int forward, double acc, +* double maxacc, int maxorder, const double *lbnd, +* const double *ubnd, int *status ) + +* Description: +* This function creates a new forward or inverse transformation for +* the supplied PolyMap (replacing any existing transformation), by +* sampling the other transformation and performing a least squares +* polynomial fit to the sample positions and values. +* +* The transformation to create is specified by the "forward" parameter. +* In what follows "X" refers to the inputs of the PolyMap, and "Y" to +* the outputs of the PolyMap. The forward transformation transforms +* input values (X) into output values (Y), and the inverse transformation +* transforms output values (Y) into input values (X). Within a PolyMap, +* each transformation is represented by an independent set of +* polynomials: Y=P_f(X) for the forward transformation and X=P_i(Y) +* for the inverse transformation. +* +* If "forward" is zero then a new inverse transformation is created by +* first finding the output values (Y) using the forward transformation +* (which must be available) at a regular grid of points (X) covering a +* rectangular region of the PolyMap's input space. The coefficients of +* the required inverse polynomial, X=P_i(Y), are chosen in order to +* minimise the sum of the squared residuals between the sampled values +* of X and P_i(Y). +* +* If "forward" is non-zero then a new forward transformation is created +* by first finding the input values (X) using the inverse transformation +* (which must be available) at a regular grid of points (Y) covering a +* rectangular region of the PolyMap's output space. The coefficients of +* the required forward polynomial, Y=P_f(X), are chosen in order to +* minimise the sum of the squared residuals between the sampled values +* of Y and P_f(X). +* +* This fitting process is performed repeatedly with increasing +* polynomial orders (starting with linear) until the target +* accuracy is achieved, or a specified maximum order is reached. If +* the target accuracy cannot be achieved even with this maximum-order +* polynomial, the best fitting maximum-order polynomial is returned so +* long as its accuracy is better than "maxacc". + +* Parameters: +* this +* The PolyMap. +* forward +* If non-zero, then the forward PolyMap transformation is +* replaced. Otherwise the inverse transformation is replaced. +* acc +* The target accuracy, expressed as a geodesic distance within +* the PolyMap's input space (if "forward" is zero) or output +* space (if "forward" is non-zero). +* maxacc +* The maximum allowed accuracy for an acceptable polynomial, +* expressed as a geodesic distance within the PolyMap's input +* space (if "forward" is zero) or output space (if "forward" is +* non-zero). +* maxorder +* The maximum allowed polynomial order. This is one more than the +* maximum power of either input axis. So for instance, a value of +* 3 refers to a quadratic polynomial. Note, cross terms with total +* powers greater than or equal to maxorder are not inlcuded in the +* fit. So the maximum number of terms in each of the fitted polynomials +* is maxorder*(maxorder+1)/2. +* lbnd +* An array holding the lower bounds of a rectangular region within +* the PolyMap's input space (if "forward" is zero) or output +* space (if "forward" is non-zero). The new polynomial will be +* evaluated over this rectangle. +* ubnd +* An array holding the upper bounds of a rectangular region within +* the PolyMap's input space (if "forward" is zero) or output +* space (if "forward" is non-zero). The new polynomial will be +* evaluated over this rectangle. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* - Non-zero if a fit was performed succesfully, to at least the +* maximum allowed accuracy. Zero if the fit failed, or an error +* occurred. + +* Notes: +* - No error is reported if the fit fails to achieve the required +* maximum accuracy. +* - An error is reported if the transformation that is not being +* replaced is not defined. +* - An error is reported if the PolyMap does not have equal numbers +* of inputs and outputs. +* - An error is reported if the PolyMap has more than 2 inputs or outputs. + +*/ + +/* Local Variables: */ + double **table; + double *cofs; + double racc; + double scales[ 4 ]; + int ndim; + int ncof; + int nsamp; + int order; + int result; + +/* Check inherited status */ + if( !astOK ) return 0; + +/* Check the PolyMap can be used. */ + ndim = astGetNin( this ); + if( astGetNout( this ) != ndim && astOK ) { + astError( AST__BADNI, "astPolyTran(%s): Supplied %s has " + "different number of inputs (%d) and outputs (%d).", + status, astGetClass( this ), astGetClass( this ), ndim, + astGetNout( this ) ); + + } else if( ndim > 2 && astOK ) { + astError( AST__BADNI, "astPolyTran(%s): Supplied %s has " + "too many inputs and outputs (%d) - must be 1 or 2.", + status, astGetClass( this ), astGetClass( this ), ndim ); + } + + if( forward != astGetInvert( this ) ){ + if( ! this->ncoeff_i && astOK ) { + astError( AST__NODEF, "astPolyTran(%s): Supplied %s has " + "no inverse transformation.", status, astGetClass( this ), + astGetClass( this ) ); + } + } else { + if( ! this->ncoeff_f && astOK ) { + astError( AST__NODEF, "astPolyTran(%s): Supplied %s has " + "no forward transformation.", status, astGetClass( this ), + astGetClass( this ) ); + } + } + +/* Check the bounds can be used. */ + if( !lbnd || !ubnd ) { + astError( AST__NODEF, "astPolyTran(%s): No upper and/or lower bounds " + "supplied.", status, astGetClass( this ) ); + } else if( lbnd[ 0 ] >= ubnd[ 0 ] && astOK ) { + astError( AST__NODEF, "astPolyTran(%s): Supplied upper " + "bound for the first axis (%g) is less than or equal to the " + "supplied lower bound (%g).", status, astGetClass( this ), + lbnd[ 0 ], ubnd[ 0 ] ); + } else if( ndim == 2 && lbnd[ 1 ] >= ubnd[ 1 ] && astOK ) { + astError( AST__NODEF, "astPolyTran(%s): Supplied upper " + "bound for the second axis (%g) is less than or equal to the " + "supplied lower bound (%g).", status, astGetClass( this ), + lbnd[ 1 ], ubnd[ 1 ] ); + } + +/* Initialise pointer to work space. */ + table = NULL; + +/* Loop over increasing polynomial orders until the required accuracy is + achieved, up to the specified maximum order. The "order" value is one more + than the maximum power in the polynomial (so a quadratic has "order" 3). */ + if( maxorder < 2 ) maxorder = 2; + for( order = 2; order <= maxorder; order++ ) { + +/* First do 2D PolyMaps. */ + if( ndim == 2 ) { + +/* Sample the requested polynomial transformation at a grid of points. This + grid covers the user-supplied region, using 2*order points on each + axis. */ + table = SamplePoly2D( this, !forward, table, lbnd, ubnd, 2*order, + &nsamp, scales, status ); + +/* Fit the polynomial. Always fit a linear polynomial ("order" 2) to any + dummy second axis. If successfull, replace the PolyMap transformation + and break out of the order loop. */ + cofs = FitPoly2D( this, forward, nsamp, acc, order, table, scales, + &ncof, &racc, status ); + +/* Now do 1D PolyMaps. */ + } else { + table = SamplePoly1D( this, !forward, table, lbnd[ 0 ], ubnd[ 0 ], + 2*order, &nsamp, scales, status ); + cofs = FitPoly1D( this, forward, nsamp, acc, order, table, scales, + &ncof, &racc, status ); + } + +/* If the fit was succesful, replace the PolyMap transformation and break + out of the order loop. */ + if( cofs && ( racc < acc || ( racc < maxacc && order == maxorder ) ) ) { + StoreArrays( this, forward, ncof, cofs, status ); + break; + } else { + cofs = astFree( cofs ); + } + } + +/* If no fit was produced, return zero. */ + result = cofs ? 1 : 0; + +/* Free resources. */ + cofs = astFree( cofs ); + table = astFreeDouble( table ); + +/* Return the result. */ + return result; +} + +static double **SamplePoly1D( AstPolyMap *this, int forward, double **table, + double lbnd, double ubnd, int npoint, int *nsamp, + double scales[2], int *status ){ +/* +* Name: +* SamplePoly1D + +* Purpose: +* Create a table of input and output positions for a 1D PolyMap. + +* Type: +* Private function. + +* Synopsis: +* double **SamplePoly1D( AstPolyMap *this, int forward, double **table, +* double lbnd, double ubnd, int npoint, int *nsamp, +* double scales[2], int *status ) + +* Description: +* This function creates a table containing samples of the requested +* polynomial transformation at a grid of input points. This grid covers +* the user-supplied region, using "npoint" points. + +* Parameters: +* this +* The PolyMap. +* forward +* If non-zero, then the forward PolyMap transformation is sampled. +* Otherwise the inverse transformation is sampled. +* table +* Pointer to a previous table created by this function, which is +* to be re-used, or NULL. +* lbnd +* The lower bounds of the region within the PolyMap's 1D input space +* (if "forward" is non-zero) or output space (if "forward" is zero). +* The new polynomial will be evaluated over this region. +* ubnd +* The upper bounds of the region within the PolyMap's 1D input space +* (if "forward" is non-zero) or output space (if "forward" is zero). +* The new polynomial will be evaluated over this region. +* npoint +* The number of points to use. +* nsamp +* Address of an int in which to return the total number of samples +* in the returned table. +* scales +* Array in which to return the scaling factors for the two columns +* of the returned table. Multiplying the returned table values by +* the scale factor produces PolyMap input or output axis values. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Pointer to an array of 2 pointers. Each of these pointers points +* to an array of "nsamp" doubles, being the sampled values for y1 +* and x1 in that order. Here x1 is the input value for the sampled +* transformation (these are spaced on the regular grid specified +* by lbnd, ubnd and npoint), and y1 is the output position produced +* by the sampled transformation. The returned values are scaled so +* that each column has an RMS value of 1.0. The scaling factors that +* convert scaled values into original values are returned in "scales". +* The returned pointer should be freed using astFreeDouble when no +* longer needed. + +*/ + +/* Local Variables: */ + AstPointSet *ps1; + AstPointSet *ps2; + double **result; + double *p0; + double *p1; + double *ptr1[ 1 ]; + double *ptr2[ 1 ]; + double delta0; + double rms; + double sum; + double val0; + int i; + int icol; + +/* Initialise returned value */ + result = table; + *nsamp = 0; + +/* Check inherited status */ + if( !astOK ) return result; + +/* Ensure we have a table of the correct size. */ + *nsamp = npoint; + if( !result ) result = astCalloc( 2, sizeof( double * ) ); + if( result ) { + for( i = 0; i < 2; i++ ) { + result[ i ] = astRealloc( result[ i ] , (*nsamp)*sizeof( double ) ); + } + } + +/* Work out the step sizes for the grid. */ + delta0 = ( ubnd - lbnd )/( npoint - 1 ); + +/* Create a PointSet to hold the grid of input positions. Use column 1 + of the table to hold the PointSet values. */ + ps1 = astPointSet( *nsamp, 1, " ", status ); + ptr1[ 0 ] = result[ 1 ]; + astSetPoints( ps1, ptr1 ); + +/* Create a PointSet to hold the grid of output positions. Use column 0 + of the table to hold the PointSet values. */ + ps2 = astPointSet( *nsamp, 1, " ", status ); + ptr2[ 0 ] = result[ 0 ]; + astSetPoints( ps2, ptr2 ); + if( astOK ) { + +/* Calculate the grid of input positions and store in the PointSet and + therefore also in the returned table. Rounding error may cause the + final point to be just over the upper bound, so we leave the final + point out of the loop and set it explicitly to the upper bound + afterwards. */ + val0 = lbnd; + p0 = ptr1[ 0 ]; + for( i = 0; i < npoint-1; i++ ) { + *(p0++) = val0; + val0 += delta0; + } + *p0 = ubnd; + +/* Transform the input grid to get the output grid. */ + (void) astTransform( this, ps1, forward, ps2 ); + +/* Scale each column in turn. */ + for( icol = 0; icol < 2; icol++ ) { + +/* Find the RMS of the values in the column. */ + sum = 0.0; + p0 = result[ icol ]; + p1 = p0 + (*nsamp); + for( ; p0 < p1; p0++ ) sum += ( *p0 )*( *p0 ); + rms = sqrt( sum/(*nsamp) ); + +/* Divide the table values by the RMS. */ + p0 = result[ icol ]; + p1 = p0 + (*nsamp); + for( ; p0 < p1; p0++ ) *p0 /= rms; + +/* Return the RMS as the scale factor. */ + scales[ icol ] = rms; + } + } + +/* Free resources */ + ps1 = astAnnul( ps1 ); + ps2 = astAnnul( ps2 ); + +/* If an error occurred, free the returned array. */ + if( !astOK ) result = astFreeDouble( result ); + +/* Return a pointer to the table. */ + return result; +} + +static double **SamplePoly2D( AstPolyMap *this, int forward, double **table, + const double *lbnd, const double *ubnd, int npoint, + int *nsamp, double scales[4], int *status ){ +/* +* Name: +* SamplePoly2D + +* Purpose: +* Create a table of input and output positions for a 2D PolyMap. + +* Type: +* Private function. + +* Synopsis: +* double **SamplePoly2D( AstPolyMap *this, int forward, double **table, +* const double *lbnd, const double *ubnd, int npoint, +* int *nsamp, double scales[4], int *status ) + +* Description: +* This function creates a table containing samples of the requested +* polynomial transformation at a grid of input points. This grid covers +* the user-supplied region, using "npoint" points on each axis. + +* Parameters: +* this +* The PolyMap. +* forward +* If non-zero, then the forward PolyMap transformation is sampled. +* Otherwise the inverse transformation is sampled. +* table +* Pointer to a previous table created by this function, which is +* to be re-used, or NULL. +* lbnd +* An array holding the lower bounds of a rectangular region within +* the PolyMap's input space (if "forward" is non-zero) or output +* space (if "forward" is zero). The new polynomial will be +* evaluated over this rectangle. +* ubnd +* An array holding the upper bounds of a rectangular region within +* the PolyMap's input space (if "forward" is non-zero) or output +* space (if "forward" is zero). The new polynomial will be +* evaluated over this rectangle. +* npoint +* The number of points along each edge of the grid. +* nsamp +* Address of an int in which to return the total number of samples +* in the returned table. +* scales +* Array in which to return the scaling factors for the four +* columns of the returned table. Multiplying the returned table +* values by the scale factor produces PolyMap input or output axis +* values. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* Pointer to an array of 4 pointers. Each of these pointers points +* to an array of "nsamp" doubles, being the sampled values for y1, +* y2, x1 or x2 in that order. Here (x1,x2) are the input values +* for the sampled transformation (these are spaced on the regular +* grid specified by lbnd, ubnd and npoint), and (y1,y2) are the +* output positions produced by the sampled transformation. The +* returned values are scaled so that each column has an RMS value +* of 1.0. The scaling factors that convert scaled values into +* original values are returned in "scales". The returned pointer +* should be freed using astFreeDouble when no longer needed. + +*/ + +/* Local Variables: */ + AstPointSet *ps1; + AstPointSet *ps2; + double **result; + double *p0; + double *p1; + double *ptr1[ 2 ]; + double *ptr2[ 2 ]; + double delta1; + double delta0; + double rms; + double sum; + double val0; + double val1; + int i; + int icol; + int j; + +/* Initialise returned value */ + result = table; + *nsamp = 0; + +/* Check inherited status */ + if( !astOK ) return result; + +/* Ensure we have a table of the correct size. */ + *nsamp = npoint*npoint; + if( !result ) result = astCalloc( 4, sizeof( double * ) ); + if( result ) { + for( i = 0; i < 4; i++ ) { + result[ i ] = astRealloc( result[ i ] , (*nsamp)*sizeof( double ) ); + } + } + +/* Work out the step sizes for the grid. */ + delta0 = ( ubnd[ 0 ] - lbnd[ 0 ] )/( npoint - 1 ); + delta1 = ( ubnd[ 1 ] - lbnd[ 1 ] )/( npoint - 1 ); + +/* Create a PointSet to hold the grid of input positions. Use columns 2 + and 3 of the table to hold the PointSet values. */ + ps1 = astPointSet( *nsamp, 2, " ", status ); + ptr1[ 0 ] = result[ 2 ]; + ptr1[ 1 ] = result[ 3 ]; + astSetPoints( ps1, ptr1 ); + +/* Create a PointSet to hold the grid of output positions. Use columns 0 + and 1 of the table to hold the PointSet values. */ + ps2 = astPointSet( *nsamp, 2, " ", status ); + ptr2[ 0 ] = result[ 0 ]; + ptr2[ 1 ] = result[ 1 ]; + astSetPoints( ps2, ptr2 ); + if( astOK ) { + +/* Calculate the grid of input positions and store in the PointSet and + therefore also in the returned table. Rounding error may cause the + final point on each axis to be just over the upper bound, so we leave + the final point out of the loop and set it explicitly to the upper bound + afterwards. */ + val0 = lbnd[ 0 ]; + p0 = ptr1[ 0 ]; + p1 = ptr1[ 1 ]; + for( i = 0; i < npoint - 1; i++ ) { + val1 = lbnd[ 1 ]; + for( j = 0; j < npoint-1; j++ ) { + *(p0++) = val0; + *(p1++) = val1; + val1 += delta1; + } + *(p0++) = val0; + *(p1++) = ubnd[ 1 ]; + val0 += delta0; + } + + + val1 = lbnd[ 1 ]; + for( j = 0; j < npoint-1; j++ ) { + *(p0++) = ubnd[ 0 ]; + *(p1++) = val1; + val1 += delta1; + } + *(p0++) = ubnd[ 0 ]; + *(p1++) = ubnd[ 1 ]; + + +/* Transform the input grid to get the output grid. */ + (void) astTransform( this, ps1, forward, ps2 ); + +/* Scale each pair of columns in turn. Use the same scale factor for + each axis in order to ensure an isotropic metric. */ + for( icol = 0; icol < 4; icol += 2 ) { + +/* Find the RMS of the values in the two columns. */ + sum = 0.0; + p0 = result[ icol ]; + p1 = p0 + (*nsamp); + for( ; p0 < p1; p0++ ) sum += ( *p0 )*( *p0 ); + + p0 = result[ icol + 1 ]; + p1 = p0 + (*nsamp); + for( ; p0 < p1; p0++ ) sum += ( *p0 )*( *p0 ); + + rms = sqrt( sum/(2*(*nsamp)) ); + +/* Divide the table values by the RMS. */ + p0 = result[ icol ]; + p1 = p0 + (*nsamp); + for( ; p0 < p1; p0++ ) *p0 /= rms; + + p0 = result[ icol + 1 ]; + p1 = p0 + (*nsamp); + for( ; p0 < p1; p0++ ) *p0 /= rms; + +/* Return the RMS as the scale factor. */ + scales[ icol ] = rms; + scales[ icol + 1 ] = rms; + } + } + +/* Free resources */ + ps1 = astAnnul( ps1 ); + ps2 = astAnnul( ps2 ); + +/* If an error occurred, free the returned array. */ + if( !astOK ) result = astFreeDouble( result ); + +/* Return a pointer to the table. */ + return result; +} + +static void SetAttrib( AstObject *this_object, const char *setting, int *status ) { +/* +* Name: +* SetAttrib + +* Purpose: +* Set an attribute value for a PolyMap. + +* Type: +* Private function. + +* Synopsis: +* #include "polymap.h" +* void SetAttrib( AstObject *this, const char *setting ) + +* Class Membership: +* PolyMap member function (over-rides the astSetAttrib protected +* method inherited from the Mapping class). + +* Description: +* This function assigns an attribute value for a PolyMap, the +* attribute and its value being specified by means of a string of +* the form: +* +* "attribute= value " +* +* Here, "attribute" specifies the attribute name and should be in +* lower case with no white space present. The value to the right +* of the "=" should be a suitable textual representation of the +* value to be assigned and this will be interpreted according to +* the attribute's data type. White space surrounding the value is +* only significant for string attributes. + +* Parameters: +* this +* Pointer to the PolyMap. +* setting +* Pointer to a null-terminated string specifying the new attribute +* value. +*/ + +/* Local Variables: */ + AstPolyMap *this; /* Pointer to the PolyMap structure */ + double dval; /* Floating point attribute value */ + int ival; /* Integer attribute value */ + int len; /* Length of setting string */ + int nc; /* Number of characters read by astSscanf */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain a pointer to the PolyMap structure. */ + this = (AstPolyMap *) this_object; + +/* Obtain the length of the setting string. */ + len = (int) strlen( setting ); + +/* Test for each recognised attribute in turn, using "astSscanf" to parse + the setting string and extract the attribute value (or an offset to + it in the case of string values). In each case, use the value set + in "nc" to check that the entire string was matched. Once a value + has been obtained, use the appropriate method to set it. */ + +/* IterInverse. */ +/* ------------ */ + if ( nc = 0, + ( 1 == astSscanf( setting, "iterinverse= %d %n", &ival, &nc ) ) + && ( nc >= len ) ) { + astSetIterInverse( this, ival ); + +/* NiterInverse. */ +/* ------------- */ + } else if ( nc = 0, + ( 1 == astSscanf( setting, "niterinverse= %d %n", &ival, &nc ) ) + && ( nc >= len ) ) { + astSetNiterInverse( this, ival ); + +/* TolInverse. */ +/* ----------- */ + } else if ( nc = 0, + ( 1 == astSscanf( setting, "tolinverse= %lg %n", &dval, &nc ) ) + && ( nc >= len ) ) { + astSetTolInverse( this, dval ); + +/* If the attribute is still not recognised, pass it on to the parent + method for further interpretation. */ + } else { + (*parent_setattrib)( this_object, setting, status ); + } +} + +static void StoreArrays( AstPolyMap *this, int forward, int ncoeff, + const double *coeff, int *status ){ +/* +* Name: +* StoreArrays + +* Purpose: +* Store the dynamic arrays for a single transformation within a PolyMap + +* Type: +* Private function. + +* Synopsis: +* #include "polymap.h" +* void StoreArrays( AstPolyMap *this, int forward, int ncoeff, +* const double *coeff, int *status ) + +* Class Membership: +* PolyMap initialiser. + +* Description: +* This function sets up the arrays within a PolyMap structure that +* describes either the forward or inverse transformation. + +* Parameters: +* this +* The PolyMap. +* forward +* If non-zero, replace the forward transformation. Otherwise, +* replace the inverse transformation. +* ncoeff +* The number of non-zero coefficients necessary to define the +* specified transformation of the PolyMap. If zero is supplied, the +* transformation will be undefined. +* coeff +* An array containing "ncof*( 2 + nin )" elements. Each group +* of "2 + nin" adjacent elements describe a single coefficient of +* the transformation. Within each such group, the first +* element is the coefficient value; the next element is the +* integer index of the PolyMap output which uses the coefficient +* within its defining polynomial (the first output has index 1); +* the remaining elements of the group give the integer powers to +* use with each input coordinate value (powers must not be +* negative). +* status +* Pointer to inherited status. +*/ + +/* Local Variables: */ + const double *group; /* Pointer to start of next coeff. description */ + int *pows; /* Pointer to powers for current coeff. */ + int gsize; /* Length of each coeff. description */ + int i; /* Loop count */ + int ico; /* Index of next coeff. for current input or output */ + int iin; /* Input index extracted from coeff. description */ + int iout; /* Output index extracted from coeff. description */ + int j; /* Loop count */ + int nin; /* Number of orignal inputs */ + int nout; /* Number of original outputs */ + int pow; /* Power extracted from coeff. description */ + +/* Check the global status. */ + if ( !astOK ) return; + +/* First Free any existing arrays. */ + FreeArrays( this, forward, status ); + +/* Get the number of inputs and outputs of the original uninverted PolyMap. + Swap the transformation to be set if the PolyMap has been inverted. */ + if( astGetInvert( this ) ) { + nout = astGetNin( this ); + nin = astGetNout( this ); + forward = !forward; + } else { + nin = astGetNin( this ); + nout = astGetNout( this ); + } + +/* Now initialise the original forward transformation, if required. */ + if( forward && ncoeff > 0 ) { + +/* Create the arrays decribing the forward transformation. */ + this->ncoeff_f = astMalloc( sizeof( int )*(size_t) nout ); + this->mxpow_f = astMalloc( sizeof( int )*(size_t) nin ); + this->power_f = astMalloc( sizeof( int ** )*(size_t) nout ); + this->coeff_f = astMalloc( sizeof( double * )*(size_t) nout ); + if( astOK ) { + +/* Initialise the count of coefficients for each output coordinate to zero. */ + for( i = 0; i < nout; i++ ) this->ncoeff_f[ i ] = 0; + +/* Initialise max power for each input coordinate to zero. */ + for( j = 0; j < nin; j++ ) this->mxpow_f[ j ] = 0; + +/* Scan through the supplied forward coefficient array, counting the + number of coefficients which relate to each output. Also find the + highest power used for each input axis. Report errors if any unusable + values are found in the supplied array. */ + group = coeff; + gsize = 2 + nin; + for( i = 0; i < ncoeff && astOK; i++, group += gsize ) { + + iout = floor( group[ 1 ] + 0.5 ); + if( iout < 1 || iout > nout ) { + astError( AST__BADCI, "astInitPolyMap(%s): Forward " + "coefficient %d referred to an illegal output " + "coordinate %d.", status, astGetClass( this ), i + 1, + iout ); + astError( AST__BADCI, "This number should be in the " + "range 1 to %d.", status, nout ); + break; + } + + this->ncoeff_f[ iout - 1 ]++; + + for( j = 0; j < nin; j++ ) { + pow = floor( group[ 2 + j ] + 0.5 ); + if( pow < 0 ) { + astError( AST__BADPW, "astInitPolyMap(%s): Forward " + "coefficient %d has a negative power (%d) " + "for input coordinate %d.", status, + astGetClass( this ), i + 1, pow, j + 1 ); + astError( AST__BADPW, "All powers should be zero or " + "positive." , status); + break; + } + if( pow > this->mxpow_f[ j ] ) this->mxpow_f[ j ] = pow; + } + } + +/* Allocate the arrays to store the input powers associated with each + coefficient, and the coefficient values. Reset the coefficient count + for each axis to zero afterwards so that we can use the array as an index + to the next vacant slot within the following loop. */ + for( i = 0; i < nout; i++ ) { + this->power_f[ i ] = astMalloc( sizeof( int * )* + (size_t) this->ncoeff_f[ i ] ); + this->coeff_f[ i ] = astMalloc( sizeof( double )* + (size_t) this->ncoeff_f[ i ] ); + this->ncoeff_f[ i ] = 0; + } + + if( astOK ) { + +/* Extract the coefficient values and powers form the supplied array and + store them in the arrays created above. */ + group = coeff; + for( i = 0; i < ncoeff && astOK; i++, group += gsize ) { + iout = floor( group[ 1 ] + 0.5 ) - 1; + ico = ( this->ncoeff_f[ iout ] )++; + this->coeff_f[ iout ][ ico ] = group[ 0 ]; + + pows = astMalloc( sizeof( int )*(size_t) nin ); + this->power_f[ iout ][ ico ] = pows; + if( astOK ) { + for( j = 0; j < nin; j++ ) { + pows[ j ] = floor( group[ 2 + j ] + 0.5 ); + } + } + } + } + } + } + +/* Now initialise the original inverse transformation, if required. */ + if( !forward && ncoeff > 0 ) { + +/* Create the arrays decribing the inverse transformation. */ + this->ncoeff_i = astMalloc( sizeof( int )*(size_t) nin ); + this->mxpow_i = astMalloc( sizeof( int )*(size_t) nout ); + this->power_i = astMalloc( sizeof( int ** )*(size_t) nin ); + this->coeff_i = astMalloc( sizeof( double * )*(size_t) nin ); + if( astOK ) { + +/* Initialise the count of coefficients for each input coordinate to zero. */ + for( i = 0; i < nin; i++ ) this->ncoeff_i[ i ] = 0; + +/* Initialise max power for each output coordinate to zero. */ + for( j = 0; j < nout; j++ ) this->mxpow_i[ j ] = 0; + +/* Scan through the supplied inverse coefficient array, counting the + number of coefficients which relate to each input. Also find the + highest power used for each output axis. Report errors if any unusable + values are found in the supplied array. */ + group = coeff; + + gsize = 2 + nout; + for( i = 0; i < ncoeff && astOK; i++, group += gsize ) { + + iin = floor( group[ 1 ] + 0.5 ); + if( iin < 1 || iin > nin ) { + astError( AST__BADCI, "astInitPolyMap(%s): Inverse " + "coefficient %d referred to an illegal input " + "coordinate %d.", status, astGetClass( this ), + i + 1, iin ); + astError( AST__BADCI, "This number should be in the " + "range 1 to %d.", status, nin ); + break; + } + + this->ncoeff_i[ iin - 1 ]++; + + for( j = 0; j < nout; j++ ) { + pow = floor( group[ 2 + j ] + 0.5 ); + if( pow < 0 ) { + astError( AST__BADPW, "astInitPolyMap(%s): Inverse " + "coefficient %d has a negative power (%d) " + "for output coordinate %d.", status, + astGetClass( this ), i + 1, pow, j + 1 ); + astError( AST__BADPW, "All powers should be zero or " + "positive." , status); + break; + } + if( pow > this->mxpow_i[ j ] ) this->mxpow_i[ j ] = pow; + } + } + +/* Allocate the arrays to store the output powers associated with each + coefficient, and the coefficient values. Reset the coefficient count + for each axis to zero afterwards so that we can use the array as an index + to the next vacant slot within the following loop. */ + for( i = 0; i < nin; i++ ) { + this->power_i[ i ] = astMalloc( sizeof( int * )* + (size_t) this->ncoeff_i[ i ] ); + this->coeff_i[ i ] = astMalloc( sizeof( double )* + (size_t) this->ncoeff_i[ i ] ); + this->ncoeff_i[ i ] = 0; + } + + if( astOK ) { + +/* Extract the coefficient values and powers form the supplied array and + store them in the arrays created above. */ + group = coeff; + for( i = 0; i < ncoeff && astOK; i++, group += gsize ) { + iin = floor( group[ 1 ] + 0.5 ) - 1; + ico = ( this->ncoeff_i[ iin ] )++; + this->coeff_i[ iin ][ ico ] = group[ 0 ]; + + pows = astMalloc( sizeof( int )*(size_t) nout ); + this->power_i[ iin ][ ico ] = pows; + if( astOK ) { + for( j = 0; j < nout; j++ ) { + pows[ j ] = floor( group[ 2 + j ] + 0.5 ); + } + } + } + } + } + } +} + +static int TestAttrib( AstObject *this_object, const char *attrib, int *status ) { +/* +* Name: +* TestAttrib + +* Purpose: +* Test if a specified attribute value is set for a PolyMap. + +* Type: +* Private function. + +* Synopsis: +* #include "polymap.h" +* int TestAttrib( AstObject *this, const char *attrib, int *status ) + +* Class Membership: +* PolyMap member function (over-rides the astTestAttrib protected +* method inherited from the Mapping class). + +* Description: +* This function returns a boolean result (0 or 1) to indicate whether +* a value has been set for one of a PolyMap's attributes. + +* Parameters: +* this +* Pointer to the PolyMap. +* attrib +* Pointer to a null-terminated string specifying the attribute +* name. This should be in lower case with no surrounding white +* space. +* status +* Pointer to the inherited status variable. + +* Returned Value: +* One if a value has been set, otherwise zero. + +* 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: */ + AstPolyMap *this; /* Pointer to the PolyMap structure */ + int result; /* Result value to return */ + +/* Initialise. */ + result = 0; + +/* Check the global error status. */ + if ( !astOK ) return result; + +/* Obtain a pointer to the PolyMap structure. */ + this = (AstPolyMap *) this_object; + +/* Check the attribute name and test the appropriate attribute. */ + +/* IterInverse. */ +/* ------------ */ + if ( !strcmp( attrib, "iterinverse" ) ) { + result = astTestIterInverse( this ); + +/* NiterInverse. */ +/* ------------- */ + } else if ( !strcmp( attrib, "niterinverse" ) ) { + result = astTestNiterInverse( this ); + +/* TolInverse. */ +/* ----------- */ + } else if ( !strcmp( attrib, "tolinverse" ) ) { + result = astTestTolInverse( this ); + +/* If the attribute is still not recognised, pass it on to the parent + method for further interpretation. */ + } else { + result = (*parent_testattrib)( this_object, attrib, status ); + } + +/* Return the result, */ + return result; +} + +static AstPointSet *Transform( AstMapping *this, AstPointSet *in, + int forward, AstPointSet *out, int *status ) { +/* +* Name: +* Transform + +* Purpose: +* Apply a PolyMap to transform a set of points. + +* Type: +* Private function. + +* Synopsis: +* #include "polymap.h" +* AstPointSet *Transform( AstMapping *this, AstPointSet *in, +* int forward, AstPointSet *out, int *status ) + +* Class Membership: +* PolyMap member function (over-rides the astTransform protected +* method inherited from the Mapping class). + +* Description: +* This function takes a PolyMap and a set of points encapsulated in a +* PointSet and transforms the points. + +* Parameters: +* this +* Pointer to the PolyMap. +* 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 PolyMap being applied. +* - The number of coordinate values per point in the output PointSet will +* equal the number of rows in the PolyMap 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 */ + AstPolyMap *map; /* Pointer to PolyMap to be applied */ + double **coeff; /* Pointer to coefficient value arrays */ + double **ptr_in; /* Pointer to input coordinate data */ + double **ptr_out; /* Pointer to output coordinate data */ + double **work; /* Pointer to exponentiated axis values */ + double *outcof; /* Pointer to next coefficient value */ + double outval; /* Output axis value */ + double term; /* Term to be added to output value */ + double xp; /* Exponentiated input axis value */ + int ***power; /* Pointer to coefficient power arrays */ + int **outpow; /* Pointer to next set of axis powers */ + int *mxpow; /* Pointer to max used power for each input */ + int *ncoeff; /* Pointer to no. of coefficients */ + int in_coord; /* Index of output coordinate */ + int ico; /* Coefficient index */ + int nc; /* No. of coefficients in polynomial */ + 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 */ + int pow; /* Next axis power */ + +/* Check the global error status. */ + if ( !astOK ) return NULL; + +/* Obtain a pointer to the PolyMap. */ + map = (AstPolyMap *) 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 ); + +/* Determine whether to apply the original forward or inverse mapping, + according to the direction specified and whether the mapping has been + inverted. */ + if ( astGetInvert( map ) ) forward = !forward; + +/* We will now extend the parent astTransform method by performing the + calculations needed to generate the output coordinate values. */ + +/* If we are using the original inverse transformatiom, and the IterInverse + attribute is non-zero, use an iterative inverse algorithm rather than any + inverse transformation defined within the PolyMap. */ + if( !forward && astGetIterInverse(map) ) { + IterInverse( map, in, result, status ); + +/* Otherwise, 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. */ + } else { + ncoord_in = astGetNcoord( in ); + ncoord_out = astGetNcoord( result ); + npoint = astGetNpoint( in ); + ptr_in = astGetPoints( in ); + ptr_out = astGetPoints( result ); + +/* Get a pointer to the arrays holding the required coefficient + values and powers, according to the direction of mapping required. */ + if ( forward ) { + ncoeff = map->ncoeff_f; + coeff = map->coeff_f; + power = map->power_f; + mxpow = map->mxpow_f; + } else { + ncoeff = map->ncoeff_i; + coeff = map->coeff_i; + power = map->power_i; + mxpow = map->mxpow_i; + } + +/* Allocate memory to hold the required powers of the input axis values. */ + work = astMalloc( sizeof( double * )*(size_t) ncoord_in ); + for( in_coord = 0; in_coord < ncoord_in; in_coord++ ) { + work[ in_coord ] = astMalloc( sizeof( double )* + (size_t) ( astMAX( 2, mxpow[in_coord]+1 ) ) ); + } + +/* Perform coordinate arithmetic. */ +/* ------------------------------ */ + if ( astOK ) { + +/* Loop to apply the polynomial to each point in turn.*/ + for ( point = 0; point < npoint; point++ ) { + +/* Find the required powers of the input axis values and store them + in the work array. Note, using a virtual method here slows the PolyMap + Transform function down by about 5%, compared to doing the equivalent + calculations in-line. But we need some way to allow the ChebyMap class + to over-ride the calculation of the powers, so we must do something + like this. If the 5% slow-down is too much, it can be reduced down to + about 2% by replacing the invocation of the astPolyPowers_ interface + function with a direct call to the implementation function itself. + This involves replacing the astPolyPowers call below with this: + + (**astMEMBER(this,PolyMap,PolyPowers))( (AstPolyMap *) this, work, ncoord_in, + mxpow, ptr_in, point, forward, status ); + + The above could be wrapped up in an alternative implementation of the + astPolyPowers macro, so that it looks the same as the existing code. + In fact, this scheme could be more widely used to speed up invocation + of virtual functions within AST. The disadvantage is that the interface + functions for some virtual methods includes some extra processing, + over and above simply invoking the implementation function. + + Of course the other way to get rid of the 5% slow down, is to + revert to using in-line code below, and then replicate this entire + function in the ChebyMap class, making suitable changes to use + Chebyshev functions in place of simple powers. But that is bad + structuring... */ + astPolyPowers( this, work, ncoord_in, mxpow, ptr_in, point, + forward ); + +/* Loop round each output. */ + for( out_coord = 0; out_coord < ncoord_out; out_coord++ ) { + +/* Initialise the output value. */ + outval = 0.0; + +/* Get pointers to the coefficients and powers for this output. */ + outcof = coeff[ out_coord ]; + outpow = power[ out_coord ]; + +/* Loop round all polynomial coefficients.*/ + nc = ncoeff[ out_coord ]; + for ( ico = 0; ico < nc && outval != AST__BAD; + ico++, outcof++, outpow++ ) { + +/* Initialise the current term to be equal to the value of the coefficient. + If it is bad, store a bad output value. */ + term = *outcof; + if( term == AST__BAD ) { + outval = AST__BAD; + +/* Otherwise, loop round all inputs */ + } else { + for( in_coord = 0; in_coord < ncoord_in; in_coord++ ) { + +/* Get the power of the current input axis value used by the current + coefficient. If it is zero, pass on. */ + pow = (*outpow)[ in_coord ]; + if( pow > 0 ) { + +/* Get the axis value raised to the appropriate power. */ + xp = work[ in_coord ][ pow ]; + +/* If bad, set the output value bad and break. */ + if( xp == AST__BAD ) { + outval = AST__BAD; + break; + +/* Otherwise multiply the current term by the exponentiated axis value. */ + } else { + term *= xp; + } + } + } + } + +/* Increment the output value by the current term of the polynomial. */ + if( outval != AST__BAD ) outval += term; + + } + +/* Store the output value. */ + ptr_out[ out_coord ][ point ] = outval; + + } + } + } + +/* Free resources. */ + for( in_coord = 0; in_coord < ncoord_in; in_coord++ ) { + work[ in_coord ] = astFree( work[ in_coord ] ); + } + work = astFree( work ); + } + +/* Return a pointer to the output PointSet. */ + 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). */ + +/* +*att++ +* Name: +* IterInverse + +* Purpose: +* Provide an iterative inverse transformation? + +* Type: +* Public attribute. + +* Synopsis: +* Integer (boolean). + +* Description: +* This attribute indicates whether the original inverse transformation +* of the PolyMap should be implemented via an iterative Newton-Raphson +* approximation that uses the forward transformation to transform +* candidate input positions until an output position is found which +* is close to the required output position. By default, an iterative +* inverse is provided if, and only if, no inverse polynomial was supplied +* when the PolyMap was constructed. +* +* Note, the term "inverse transformation" here refers to the inverse +* transformation of the original PolyMap, ignoring any subsequent +* inversions. Also, "input" and "output" refer to the inputs and +* outputs of the original PolyMap. + +* The NiterInverse and TolInverse attributes provide parameters that +* control the behaviour of the inverse approximation method. +* +* The iterative inverse returns AST__BAD axis values at positions +* for which the inverse transformation is undefined. For instance, +* if the forward transformation is y = x*x, the iterative inverse +* will return x = AST__BAD at y = -1. If the inverse transformation +* is multiply defined, the position returned by the iterative inverse +* will be the position of the solution that is closest to the +* supplied position. For instance, using the above example, y = x*x, +* the iterative inverse will return x = +2 at y = 4, because x = +2 +* is the closest solution to 4 (the other solution is x = -2). + +* Applicability: +* PolyMap +* All PolyMaps have this attribute. +* ChebyMap +* The ChebyMap class does not currently provide an option for an +* iterative inverse, and so the IterInverse value is always zero. +* Setting or clearing the IterInverse attribute of a ChebyMap has +* no effect. + +* Notes: +* - The transformation replaced by the iterative algorithm is the +* transformation from the original PolyMap output space to the +* original PolyMap input space (i.e. the input and output spaces as +* defined by the arguments of the PolyMap constructor). This is still +* the case even if the PolyMap has subsequently been inverted. In +* other words if a PolyMap is created and then inverted, setting +* the IterInverse to a non-zero value will replace the forward +* transformation of the inverted PolyMap (i.e. the inverse +* transformation of the original PolyMap). It is not possible to +* replace the other transformation (i.e. from the original PolyMap +* input space to the original PolyMap output space) with an iterative +* algorithm. +* - If a PolyMap that has an iterative inverse transformation is +* subsequently inverted, the inverted PolyMap will have an iterative +* forward transformation. +* - An iterative inverse can only be used if the PolyMap has equal +* numbers of inputs and outputs, as given by the Nin and Nout +* attributes. An error will be reported if IterInverse is set non-zero +* for a PolyMap that does not meet this requirement. + +*att-- +*/ +astMAKE_CLEAR(PolyMap,IterInverse,iterinverse,-INT_MAX) +astMAKE_GET(PolyMap,IterInverse,int,0,( ( this->iterinverse == -INT_MAX ) ? + (this->ncoeff_i == 0) : this->iterinverse )) +astMAKE_SET(PolyMap,IterInverse,int,iterinverse, + (((astGetNin(this)==astGetNout(this))||!value)?((value?1:0)):(astError(AST__ATTIN,"astSetIterInverse(%s):" + "Cannot use an iterative inverse because the %s has unequal numbers of " + "inputs and outputs.", status, astGetClass(this),astGetClass(this)),this->iterinverse))) +astMAKE_TEST(PolyMap,IterInverse,( this->iterinverse != -INT_MAX )) + +/* NiterInverse. */ +/* --------- */ +/* +*att++ +* Name: +* NiterInverse + +* Purpose: +* Maximum number of iterations for the iterative inverse transformation. + +* Type: +* Public attribute. + +* Synopsis: +* Integer. + +* Description: +* This attribute controls the iterative inverse transformation +* used if the IterInverse attribute is non-zero. +* +* Its value gives the maximum number of iterations of the +* Newton-Raphson algorithm to be used for each transformed position. +* The default value is 4. See also attribute TolInverse. + +* Applicability: +* PolyMap +* All PolyMaps have this attribute. + +*att-- +*/ +astMAKE_CLEAR(PolyMap,NiterInverse,niterinverse,-INT_MAX) +astMAKE_GET(PolyMap,NiterInverse,int,0,( this->niterinverse == -INT_MAX ? 4 : this->niterinverse)) +astMAKE_SET(PolyMap,NiterInverse,int,niterinverse,value) +astMAKE_TEST(PolyMap,NiterInverse,( this->niterinverse != -INT_MAX )) + +/* TolInverse. */ +/* ----------- */ +/* +*att++ +* Name: +* TolInverse + +* Purpose: +* Target relative error for the iterative inverse transformation. + +* Type: +* Public attribute. + +* Synopsis: +* Floating point. + +* Description: +* This attribute controls the iterative inverse transformation +* used if the IterInverse attribute is non-zero. +* +* Its value gives the target relative error in the axis values of +* each transformed position. Further iterations will be performed +* until the target relative error is reached, or the maximum number +* of iterations given by attribute NiterInverse is reached. + +* The default value is 1.0E-6. + +* Applicability: +* PolyMap +* All PolyMaps have this attribute. +*att-- +*/ +astMAKE_CLEAR(PolyMap,TolInverse,tolinverse,AST__BAD) +astMAKE_GET(PolyMap,TolInverse,double,0.0,( this->tolinverse == AST__BAD ? 1.0E-6 : this->tolinverse)) +astMAKE_SET(PolyMap,TolInverse,double,tolinverse,value) +astMAKE_TEST(PolyMap,TolInverse,( this->tolinverse != AST__BAD )) + +/* Copy constructor. */ +/* ----------------- */ +static void Copy( const AstObject *objin, AstObject *objout, int *status ) { +/* +* Name: +* Copy + +* Purpose: +* Copy constructor for PolyMap objects. + +* Type: +* Private function. + +* Synopsis: +* void Copy( const AstObject *objin, AstObject *objout, int *status ) + +* Description: +* This function implements the copy constructor for PolyMap 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 +* coefficients associated with the input PolyMap. +*/ + + +/* Local Variables: */ + AstPolyMap *in; /* Pointer to input PolyMap */ + AstPolyMap *out; /* Pointer to output PolyMap */ + int nin; /* No. of input coordinates */ + int nout; /* No. of output coordinates */ + int i; /* Loop count */ + int j; /* Loop count */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain pointers to the input and output PolyMaps. */ + in = (AstPolyMap *) objin; + out = (AstPolyMap *) 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->ncoeff_f = NULL; + out->power_f = NULL; + out->coeff_f = NULL; + out->mxpow_f = NULL; + + out->ncoeff_i = NULL; + out->power_i = NULL; + out->coeff_i = NULL; + out->mxpow_i = NULL; + + out->jacobian = NULL; + out->lintrunc = NULL; + +/* Get the number of inputs and outputs of the uninverted Mapping. */ + nin = ( (AstMapping *) in )->nin; + nout = ( (AstMapping *) in )->nout; + +/* Copy the number of coefficients associated with each output of the + forward transformation of the uninverted Mapping. */ + if( in->ncoeff_f ) { + out->ncoeff_f = (int *) astStore( NULL, (void *) in->ncoeff_f, + sizeof( int )*(size_t) nout ); + +/* Copy the maximum power of each input axis value used by the forward + transformation. */ + out->mxpow_f = (int *) astStore( NULL, (void *) in->mxpow_f, + sizeof( int )*(size_t) nin ); + +/* Copy the coefficient values used by the forward transformation. */ + if( in->coeff_f ) { + out->coeff_f = astMalloc( sizeof( double * )*(size_t) nout ); + if( astOK ) { + for( i = 0; i < nout; i++ ) { + out->coeff_f[ i ] = (double *) astStore( NULL, (void *) in->coeff_f[ i ], + sizeof( double )*(size_t) in->ncoeff_f[ i ] ); + } + } + } + +/* Copy the input axis powers associated with each coefficient of the forward + transformation. */ + if( in->power_f ) { + out->power_f = astMalloc( sizeof( int ** )*(size_t) nout ); + if( astOK ) { + for( i = 0; i < nout; i++ ) { + out->power_f[ i ] = astMalloc( sizeof( int * )*(size_t) in->ncoeff_f[ i ] ); + if( astOK ) { + for( j = 0; j < in->ncoeff_f[ i ]; j++ ) { + out->power_f[ i ][ j ] = (int *) astStore( NULL, (void *) in->power_f[ i ][ j ], + sizeof( int )*(size_t) nin ); + } + } + } + } + } + } + +/* Do the same for the inverse transformation. */ + if( in->ncoeff_i ) { + out->ncoeff_i = (int *) astStore( NULL, (void *) in->ncoeff_i, + sizeof( int )*(size_t) nin ); + + out->mxpow_i = (int *) astStore( NULL, (void *) in->mxpow_i, + sizeof( int )*(size_t) nout ); + + if( in->coeff_i ) { + out->coeff_i = astMalloc( sizeof( double * )*(size_t) nin ); + if( astOK ) { + for( i = 0; i < nin; i++ ) { + out->coeff_i[ i ] = (double *) astStore( NULL, (void *) in->coeff_i[ i ], + sizeof( double )*(size_t) in->ncoeff_i[ i ] ); + } + } + } + + if( in->power_i ) { + out->power_i = astMalloc( sizeof( int ** )*(size_t) nin ); + if( astOK ) { + for( i = 0; i < nin; i++ ) { + out->power_i[ i ] = astMalloc( sizeof( int * )*(size_t) in->ncoeff_i[ i ] ); + if( astOK ) { + for( j = 0; j < in->ncoeff_i[ i ]; j++ ) { + out->power_i[ i ][ j ] = (int *) astStore( NULL, (void *) in->power_i[ i ][ j ], + sizeof( int )*(size_t) nout ); + } + } + } + } + } + } + +/* Copy the linear truncation of the PolyMap - if it has been found. */ + if( in->lintrunc ) out->lintrunc = astCopy( in->lintrunc ); + +/* If an error has occurred, free all the resources allocated above. */ + if( !astOK ) { + FreeArrays( out, 1, status ); + FreeArrays( out, 0, status ); + } + + return; + +} + +/* Destructor. */ +/* ----------- */ +static void Delete( AstObject *obj, int *status ) { +/* +* Name: +* Delete + +* Purpose: +* Destructor for PolyMap objects. + +* Type: +* Private function. + +* Synopsis: +* void Delete( AstObject *obj, int *status ) + +* Description: +* This function implements the destructor for PolyMap 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: */ + AstPolyMap *this; + int nc; + int ic; + int lstat; + int error; + +/* Obtain a pointer to the PolyMap structure. */ + this = (AstPolyMap *) obj; + +/* Free the arrays. */ + FreeArrays( this, 1, status ); + FreeArrays( this, 0, status ); + +/* Free the resources used to store the Jacobian of the forward + transformation. */ + if( this->jacobian ) { + +/* Get the number of PolyMap inputs. We need to clear any error status + first since astGetNin returns zero if an error has occurred. The + Jacobian will only be non-NULL if the number of inputs and outputs + are equal. */ + error = !astOK; + if( error ) { + lstat = astStatus; + astClearStatus; + } + nc = astGetNin( this ); + if( error ) astSetStatus( lstat ); + + for( ic = 0; ic < nc; ic++ ) { + (this->jacobian)[ ic ] = astAnnul( (this->jacobian)[ ic ] ); + } + this->jacobian = astFree( this->jacobian ); + } + +/* Free the linear truncation. */ + if( this->lintrunc ) this->lintrunc = astAnnul( this->lintrunc ); +} + +/* Dump function. */ +/* -------------- */ +static void Dump( AstObject *this_object, AstChannel *channel, int *status ) { +/* +* Name: +* Dump + +* Purpose: +* Dump function for PolyMap 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 PolyMap class to an output Channel. + +* Parameters: +* this +* Pointer to the PolyMap 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: */ + AstPolyMap *this; /* Pointer to the PolyMap structure */ + char buff[ KEY_LEN + 1 ]; /* Buffer for keyword string */ + char comm[ 100 ]; /* Buffer for comment string */ + double dval; /* Floating point attribute value */ + int i; /* Loop index */ + int iv; /* Vectorised keyword index */ + int ival; /* Integer value */ + int j; /* Loop index */ + int k; /* Loop index */ + int nin; /* No. of input coords */ + int nout; /* No. of output coords */ + int set; /* Attribute value set? */ + +/* Check the global error status. */ + if ( !astOK ) return; + +/* Obtain a pointer to the PolyMap structure. */ + this = (AstPolyMap *) this_object; + +/* Find the number of inputs and outputs of the uninverted Mapping. */ + nin = ( (AstMapping *) this )->nin; + nout = ( (AstMapping *) this )->nout; + +/* Write out values representing the instance variables for the + PolyMap class. */ + +/* First do the forward transformation arrays. Check they are used. */ + if( this->ncoeff_f ) { + +/* Store the maximum power of each input axis value used by the forward + transformation. */ + for( i = 0; i < nin; i++ ){ + (void) sprintf( buff, "MPF%d", i + 1 ); + (void) sprintf( comm, "Max. power of input %d in any forward polynomial", i + 1 ); + astWriteInt( channel, buff, 1, 1, (this->mxpow_f)[ i ], comm ); + } + +/* Store the number of coefficients associated with each output of the forward + transformation. */ + for( i = 0; i < nout; i++ ){ + (void) sprintf( buff, "NCF%d", i + 1 ); + (void) sprintf( comm, "No. of coeff.s for forward polynomial %d", i + 1 ); + astWriteInt( channel, buff, 1, 1, (this->ncoeff_f)[ i ], comm ); + } + +/* Store the coefficient values used by the forward transformation. */ + iv = 1; + for( i = 0; i < nout; i++ ){ + for( j = 0; j < this->ncoeff_f[ i ]; j++, iv++ ){ + if( (this->coeff_f)[ i ][ j ] != AST__BAD ) { + (void) sprintf( buff, "CF%d", iv ); + (void) sprintf( comm, "Coeff %d of forward polynomial %d", j + 1, i + 1 ); + astWriteDouble( channel, buff, 1, 1, (this->coeff_f)[ i ][ j ], comm ); + } + } + } + +/* Store the input axis powers associated with each coefficient of the forward + transformation. */ + iv = 1; + for( i = 0; i < nout; i++ ){ + for( j = 0; j < this->ncoeff_f[ i ]; j++ ){ + for( k = 0; k < nin; k++, iv++ ){ + if( (this->power_f)[ i ][ j ][ k ] > 0 ) { + (void) sprintf( buff, "PF%d", iv ); + (void) sprintf( comm, "Power of i/p %d for coeff %d of fwd poly %d", k + 1, j + 1, i + 1 ); + astWriteDouble( channel, buff, 1, 1, (this->power_f)[ i ][ j ][ k ], comm ); + } + } + } + } + } + +/* Now do the inverse transformation arrays. Check they are used. */ + if( this->ncoeff_i ) { + +/* Store the maximum power of each output axis value used by the inverse + transformation. */ + for( i = 0; i < nout; i++ ){ + (void) sprintf( buff, "MPI%d", i + 1 ); + (void) sprintf( comm, "Max. power of output %d in any inverse polynomial", i + 1 ); + astWriteInt( channel, buff, 1, 1, (this->mxpow_i)[ i ], comm ); + } + +/* Store the number of coefficients associated with each input of the inverse + transformation. */ + for( i = 0; i < nin; i++ ){ + (void) sprintf( buff, "NCI%d", i + 1 ); + (void) sprintf( comm, "No. of coeff.s for inverse polynomial %d", i + 1 ); + astWriteInt( channel, buff, 1, 1, (this->ncoeff_i)[ i ], comm ); + } + +/* Store the coefficient values used by the inverse transformation. */ + iv = 1; + for( i = 0; i < nin; i++ ){ + for( j = 0; j < this->ncoeff_i[ i ]; j++, iv++ ){ + if( (this->coeff_i)[ i ][ j ] != AST__BAD ) { + (void) sprintf( buff, "CI%d", iv ); + (void) sprintf( comm, "Coeff %d of inverse polynomial %d", j + 1, i + 1 ); + astWriteDouble( channel, buff, 1, 1, (this->coeff_i)[ i ][ j ], comm ); + } + } + } + +/* Store the output axis powers associated with each coefficient of the inverse + transformation. */ + iv = 1; + for( i = 0; i < nin; i++ ){ + for( j = 0; j < this->ncoeff_i[ i ]; j++ ){ + for( k = 0; k < nout; k++, iv++ ){ + if( (this->power_i)[ i ][ j ][ k ] > 0 ) { + (void) sprintf( buff, "PI%d", iv ); + (void) sprintf( comm, "Power of o/p %d for coeff %d of inv poly %d", k + 1, j + 1, i + 1 ); + astWriteDouble( channel, buff, 1, 1, (this->power_i)[ i ][ j ][ k ], comm ); + } + } + } + } + } + +/* Use an iterative inverse? */ + set = TestIterInverse( this, status ); + ival = set ? GetIterInverse( this, status ) : astGetIterInverse( this ); + astWriteInt( channel, "IterInv", set, 0, ival, ival ? "Use an iterative inverse" : "Do not use an iterative inverse" ); + +/* Max number of iterations for iterative inverse. */ + set = TestNiterInverse( this, status ); + ival = set ? GetNiterInverse( this, status ) : astGetNiterInverse( this ); + astWriteInt( channel, "NiterInv", set, 0, ival, "Max number of iterations for iterative inverse" ); + +/* Target relative error for iterative inverse. */ + set = TestTolInverse( this, status ); + dval = set ? GetTolInverse( this, status ) : astGetTolInverse( this ); + astWriteDouble( channel, "TolInv", set, 0, dval, "Target relative error for iterative inverse" ); + +/* Undefine macros local to this function. */ +#undef KEY_LEN +} + +/* Standard class functions. */ +/* ========================= */ +/* Implement the astIsAPolyMap and astCheckPolyMap functions using the macros + defined for this purpose in the "object.h" header file. */ +astMAKE_ISA(PolyMap,Mapping) +astMAKE_CHECK(PolyMap) + +AstPolyMap *astPolyMap_( int nin, int nout, int ncoeff_f, const double coeff_f[], + int ncoeff_i, const double coeff_i[], const char *options, int *status, ...){ +/* +*++ +* Name: +c astPolyMap +f AST_POLYMAP + +* Purpose: +* Create a PolyMap. + +* Type: +* Public function. + +* Synopsis: +c #include "polymap.h" +c AstPolyMap *astPolyMap( int nin, int nout, int ncoeff_f, const double coeff_f[], +c int ncoeff_i, const double coeff_i[], +c const char *options, ... ) +f RESULT = AST_POLYMAP( NIN, NOUT, NCOEFF_F, COEFF_F, NCOEFF_I, COEFF_I, +f OPTIONS, STATUS ) + +* Class Membership: +* PolyMap constructor. + +* Description: +* This function creates a new PolyMap and optionally initialises +* its attributes. +* +* A PolyMap is a form of Mapping which performs a general polynomial +* transformation. Each output coordinate is a polynomial function of +* all the input coordinates. The coefficients are specified separately +* for each output coordinate. The forward and inverse transformations +* are defined independantly by separate sets of coefficients. If no +* inverse transformation is supplied, the default behaviour is to use +* an iterative method to evaluate the inverse based only on the forward +* transformation (see attribute IterInverse). + +* Parameters: +c nin +f NIN = INTEGER (Given) +* The number of input coordinates. +c nout +f NOUT = INTEGER (Given) +* The number of output coordinates. +c ncoeff_f +f NCOEFF_F = INTEGER (Given) +* The number of non-zero coefficients necessary to define the +* forward transformation of the PolyMap. If zero is supplied, the +* forward transformation will be undefined. +c coeff_f +f COEFF_F( * ) = DOUBLE PRECISION (Given) +* An array containing +c "ncoeff_f*( 2 + nin )" elements. Each group of "2 + nin" +f "NCOEFF_F*( 2 + NIN )" elements. Each group of "2 + NIN" +* adjacent elements describe a single coefficient of the forward +* transformation. Within each such group, the first element is the +* coefficient value; the next element is the integer index of the +* PolyMap output which uses the coefficient within its defining +* polynomial (the first output has index 1); the remaining elements +* of the group give the integer powers to use with each input +* coordinate value (powers must not be negative, and floating +* point values are rounded to the nearest integer). +c If "ncoeff_f" is zero, a NULL pointer may be supplied for "coeff_f". +* +* For instance, if the PolyMap has 3 inputs and 2 outputs, each group +* consisting of 5 elements, A groups such as "(1.2, 2.0, 1.0, 3.0, 0.0)" +* describes a coefficient with value 1.2 which is used within the +* definition of output 2. The output value is incremented by the +* product of the coefficient value, the value of input coordinate +* 1 raised to the power 1, and the value of input coordinate 2 raised +* to the power 3. Input coordinate 3 is not used since its power is +* specified as zero. As another example, the group "(-1.0, 1.0, +* 0.0, 0.0, 0.0 )" describes adds a constant value -1.0 onto +* output 1 (it is a constant value since the power for every input +* axis is given as zero). +* +c Each final output coordinate value is the sum of the "ncoeff_f" terms +c described by the "ncoeff_f" groups within the supplied array. +f Each final output coordinate value is the sum of the "NCOEFF_F" terms +f described by the "NCOEFF_F" groups within the supplied array. +c ncoeff_i +f NCOEFF_I = INTEGER (Given) +* The number of non-zero coefficients necessary to define the +* inverse transformation of the PolyMap. If zero is supplied, the +* default behaviour is to use an iterative method to evaluate the +* inverse based only on the forward transformation (see attribute +* IterInverse). +c coeff_i +f COEFF_I( * ) = DOUBLE PRECISION (Given) +* An array containing +c "ncoeff_i*( 2 + nout )" elements. Each group of "2 + nout" +f "NCOEFF_I*( 2 + NOUT )" elements. Each group of "2 + NOUT" +* adjacent elements describe a single coefficient of the inverse +c transformation, using the same schame as "coeff_f", +f transformation, using the same schame as "COEFF_F", +* except that "inputs" and "outputs" are transposed. +c If "ncoeff_i" is zero, a NULL pointer may be supplied for "coeff_i". +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 PolyMap. 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 PolyMap. 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 astPolyMap() +f AST_POLYMAP = INTEGER +* A pointer to the new PolyMap. + +* Notes: +* - 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. +*-- +*/ + +/* Local Variables: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstPolyMap *new; /* Pointer to new PolyMap */ + 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 PolyMap, allocating memory and initialising the + virtual function table as well if necessary. */ + new = astInitPolyMap( NULL, sizeof( AstPolyMap ), !class_init, + &class_vtab, "PolyMap", nin, nout, + ncoeff_f, coeff_f, ncoeff_i, coeff_i ); + +/* 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 PolyMap'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 PolyMap. */ + return new; +} + +AstPolyMap *astPolyMapId_( int nin, int nout, int ncoeff_f, const double coeff_f[], + int ncoeff_i, const double coeff_i[], const char *options, ... ){ +/* +* Name: +* astPolyMapId_ + +* Purpose: +* Create a PolyMap. + +* Type: +* Private function. + +* Synopsis: +* #include "polymap.h" +* AstPolyMap *astPolyMap( int nin, int nout, int ncoeff_f, const double coeff_f[], +* int ncoeff_i, const double coeff_i[], +* const char *options, ... ) + +* Class Membership: +* PolyMap constructor. + +* Description: +* This function implements the external (public) interface to the +* astPolyMap constructor function. It returns an ID value (instead +* of a true C pointer) to external users, and must be provided +* because astPolyMap_ 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 astPolyMap_ 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 astPolyMap_. + +* Returned Value: +* The ID value associated with the new PolyMap. +*/ + +/* Local Variables: */ + astDECLARE_GLOBALS /* Pointer to thread-specific global data */ + AstPolyMap *new; /* Pointer to new PolyMap */ + 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 PolyMap, allocating memory and initialising the + virtual function table as well if necessary. */ + new = astInitPolyMap( NULL, sizeof( AstPolyMap ), !class_init, + &class_vtab, "PolyMap", nin, nout, + ncoeff_f, coeff_f, ncoeff_i, coeff_i ); + +/* 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 PolyMap'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 PolyMap. */ + return astMakeId( new ); +} + +AstPolyMap *astInitPolyMap_( void *mem, size_t size, int init, + AstPolyMapVtab *vtab, const char *name, + int nin, int nout, int ncoeff_f, const double coeff_f[], + int ncoeff_i, const double coeff_i[], int *status ){ +/* +*+ +* Name: +* astInitPolyMap + +* Purpose: +* Initialise a PolyMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "polymap.h" +* AstPolyMap *astInitPolyMap( void *mem, size_t size, int init, +* AstPolyMapVtab *vtab, const char *name, +* int nin, int nout, int ncoeff_f, +* const double coeff_f[], int ncoeff_i, +* const double coeff_i[] ) + +* Class Membership: +* PolyMap initialiser. + +* Description: +* This function is provided for use by class implementations to initialise +* a new PolyMap object. It allocates memory (if necessary) to accommodate +* the PolyMap plus any additional data associated with the derived class. +* It then initialises a PolyMap 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 PolyMap at the start of the memory passed via the +* "vtab" parameter. + +* Parameters: +* mem +* A pointer to the memory in which the PolyMap is to be initialised. +* This must be of sufficient size to accommodate the PolyMap data +* (sizeof(PolyMap)) 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 PolyMap (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 PolyMap +* structure, so a valid value must be supplied even if not required for +* allocating memory. +* init +* A logical flag indicating if the PolyMap'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 PolyMap. +* 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. +* ncoeff_f +* The number of non-zero coefficients necessary to define the +* forward transformation of the PolyMap. If zero is supplied, the +* forward transformation will be undefined. +* coeff_f +* An array containing "ncoeff_f*( 2 + nin )" elements. Each group +* of "2 + nin" adjacent elements describe a single coefficient of +* the forward transformation. Within each such group, the first +* element is the coefficient value; the next element is the +* integer index of the PolyMap output which uses the coefficient +* within its defining polynomial (the first output has index 1); +* the remaining elements of the group give the integer powers to +* use with each input coordinate value (powers must not be +* negative) +* +* For instance, if the PolyMap has 3 inputs and 2 outputs, each group +* consisting of 5 elements, A groups such as "(1.2, 2.0, 1.0, 3.0, 0.0)" +* describes a coefficient with value 1.2 which is used within the +* definition of output 2. The output value is incremented by the +* product of the coefficient value, the value of input coordinate +* 1 raised to the power 1, and the value of input coordinate 2 raised +* to the power 3. Input coordinate 3 is not used since its power is +* specified as zero. As another example, the group "(-1.0, 1.0, +* 0.0, 0.0, 0.0 )" describes adds a constant value -1.0 onto +* output 1 (it is a constant value since the power for every input +* axis is given as zero). +* +* Each final output coordinate value is the sum of the "ncoeff_f" terms +* described by the "ncoeff_f" groups within the supplied array. +* ncoeff_i +* The number of non-zero coefficients necessary to define the +* inverse transformation of the PolyMap. If zero is supplied, the +* inverse transformation will be undefined. +* coeff_i +* An array containing +* "ncoeff_i*( 2 + nout )" elements. Each group of "2 + nout" +* adjacent elements describe a single coefficient of the inverse +* transformation, using the same schame as "coeff_f", except that +* "inputs" and "outputs" are transposed. + +* Returned Value: +* A pointer to the new PolyMap. + +* 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: */ + AstPolyMap *new; /* Pointer to new PolyMap */ + +/* Check the global status. */ + if ( !astOK ) return NULL; + +/* If necessary, initialise the virtual function table. */ + if ( init ) astInitPolyMapVtab( vtab, name ); + +/* Initialise a Mapping structure (the parent class) as the first component + within the PolyMap structure, allocating memory if necessary. Specify that + the Mapping should be defined in both the forward and inverse directions. */ + new = (AstPolyMap *) astInitMapping( mem, size, 0, + (AstMappingVtab *) vtab, name, + nin, nout, 1, 1 ); + if ( astOK ) { + +/* Initialise the PolyMap data. */ +/* ---------------------------- */ + +/* First initialise the pointers in case of errors. */ + new->ncoeff_f = NULL; + new->power_f = NULL; + new->coeff_f = NULL; + new->mxpow_f = NULL; + + new->ncoeff_i = NULL; + new->power_i = NULL; + new->coeff_i = NULL; + new->mxpow_i = NULL; + +/* Store the forward transformation. */ + StoreArrays( new, 1, ncoeff_f, coeff_f, status ); + +/* Store the inverse transformation. */ + StoreArrays( new, 0, ncoeff_i, coeff_i, status ); + +/* Other class attributes. */ + new->iterinverse = -INT_MAX; + new->niterinverse = -INT_MAX; + new->tolinverse = AST__BAD; + new->jacobian = NULL; + new->lintrunc = NULL; + +/* If an error occurred, clean up by deleting the new PolyMap. */ + if ( !astOK ) new = astDelete( new ); + } + +/* Return a pointer to the new PolyMap. */ + return new; +} + +AstPolyMap *astLoadPolyMap_( void *mem, size_t size, + AstPolyMapVtab *vtab, const char *name, + AstChannel *channel, int *status ) { +/* +*+ +* Name: +* astLoadPolyMap + +* Purpose: +* Load a PolyMap. + +* Type: +* Protected function. + +* Synopsis: +* #include "polymap.h" +* AstPolyMap *astLoadPolyMap( void *mem, size_t size, +* AstPolyMapVtab *vtab, const char *name, +* AstChannel *channel ) + +* Class Membership: +* PolyMap loader. + +* Description: +* This function is provided to load a new PolyMap 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 +* PolyMap 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 PolyMap at the start of the memory +* passed via the "vtab" parameter. + + +* Parameters: +* mem +* A pointer to the memory into which the PolyMap is to be +* loaded. This must be of sufficient size to accommodate the +* PolyMap data (sizeof(PolyMap)) 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 PolyMap (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 PolyMap 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(AstPolyMap) is used instead. +* vtab +* Pointer to the start of the virtual function table to be +* associated with the new PolyMap. If this is NULL, a pointer +* to the (static) virtual function table for the PolyMap 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 "PolyMap" is used instead. + +* Returned Value: +* A pointer to the new PolyMap. + +* 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: */ + AstPolyMap *new; /* Pointer to the new PolyMap */ + char buff[ KEY_LEN + 1 ]; /* Buffer for keyword string */ + int i; /* Loop index */ + int iv; /* Vectorised keyword index */ + int j; /* Loop index */ + int k; /* Loop index */ + int nin; /* No. of input coords */ + int nout; /* No. of output coords */ + int undef; /* Is the transformation undefined? */ + +/* 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 PolyMap. In this case the + PolyMap belongs to this class, so supply appropriate values to be + passed to the parent class loader (and its parent, etc.). */ + if ( !vtab ) { + size = sizeof( AstPolyMap ); + vtab = &class_vtab; + name = "PolyMap"; + +/* If required, initialise the virtual function table for this class. */ + if ( !class_init ) { + astInitPolyMapVtab( 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 PolyMap. */ + new = astLoadMapping( mem, size, (AstMappingVtab *) vtab, name, + channel ); + + if ( astOK ) { + +/* Get the number of inputs and outputs for the uninverted Mapping. */ + nin = ( (AstMapping *) new )->nin; + nout = ( (AstMapping *) new )->nout; + +/* Read input data. */ +/* ================ */ +/* Request the input Channel to read all the input data appropriate to + this class into the internal "values list". */ + astReadClassData( channel, "PolyMap" ); + +/* Allocate memory to hold the forward arrays. */ + new->ncoeff_f = astMalloc( sizeof( int )*(size_t) nout ); + new->mxpow_f = astMalloc( sizeof( int )*(size_t) nin ); + new->power_f = astMalloc( sizeof( int ** )*(size_t) nout ); + new->coeff_f = astMalloc( sizeof( double * )*(size_t) nout ); + if( astOK ) { + +/* Assume the forward transformation is defined. */ + undef = 0; + +/* Get the maximum power of each input axis value used by the forward + transformation. Set a flag "undef" if no values relating to the + forward transformation are found (this indicates that the forward + transformation is not defined). */ + for( i = 0; i < nin && !undef; i++ ){ + (void) sprintf( buff, "mpf%d", i + 1 ); + (new->mxpow_f)[ i ] = astReadInt( channel, buff, INT_MAX ); + if( (new->mxpow_f)[ i ] == INT_MAX ) undef = 1; + } + +/* Get the number of coefficients associated with each output of the forward + transformation. */ + for( i = 0; i < nout && !undef; i++ ){ + (void) sprintf( buff, "ncf%d", i + 1 ); + (new->ncoeff_f)[ i ] = astReadInt( channel, buff, INT_MAX ); + if( (new->ncoeff_f)[ i ] == INT_MAX ) undef = 1; + } + +/* Get the coefficient values used by the forward transformation. This + uses new style vectorised key names if available. Otherwise it uses + old style indexed names (which were superceded by vectorised names + because they are shorter and so work better with FitsChans). */ + iv = 0; + for( i = 0; i < nout && !undef; i++ ){ + + (new->coeff_f)[ i ] = astMalloc( sizeof( double )* + (size_t) new->ncoeff_f[ i ] ); + if( astOK ) { + for( j = 0; j < new->ncoeff_f[ i ]; j++ ){ + (void) sprintf( buff, "cf%d", ++iv ); + (new->coeff_f)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD ); + if( (new->coeff_f)[ i ][ j ] == AST__BAD ) { + (void) sprintf( buff, "cf%d_%d", i + 1, j + 1 ); + (new->coeff_f)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD ); + } + } + } + } + +/* Get the input axis powers associated with each coefficient of the forward + transformation. */ + iv = 0; + for( i = 0; i < nout && !undef; i++ ){ + (new->power_f)[ i ] = astMalloc( sizeof( int * )* + (size_t) new->ncoeff_f[ i ] ); + if( astOK ) { + for( j = 0; j < new->ncoeff_f[ i ]; j++ ){ + (new->power_f)[ i ][ j ] = astMalloc( sizeof( int )* (size_t) nin ); + if( astOK ) { + for( k = 0; k < nin; k++ ){ + (void) sprintf( buff, "pf%d", ++iv ); + (new->power_f)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 ); + if( (new->power_f)[ i ][ j ][ k ] == 0 ) { + (void) sprintf( buff, "pf%d_%d_%d", i + 1, j + 1, k + 1 ); + (new->power_f)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 ); + } + } + } + } + } + } + +/* Free the arrays if the forward transformation is undefined. */ + if( undef ) { + new->ncoeff_f = astFree( new->ncoeff_f ); + new->mxpow_f = astFree( new->mxpow_f ); + new->power_f = astFree( new->power_f ); + new->coeff_f = astFree( new->coeff_f ); + } + } + +/* Allocate memory to hold the inverse arrays. */ + new->ncoeff_i = astMalloc( sizeof( int )*(size_t) nin ); + new->mxpow_i = astMalloc( sizeof( int )*(size_t) nout ); + new->power_i = astMalloc( sizeof( int ** )*(size_t) nin ); + new->coeff_i = astMalloc( sizeof( double * )*(size_t) nin ); + if( astOK ) { + +/* Assume the inverse transformation is defined. */ + undef = 0; + +/* Get the maximum power of each output axis value used by the inverse + transformation. Set a flag "undef" if no values relating to the + inverse transformation are found (this indicates that the inverse + transformation is not defined). */ + for( i = 0; i < nout && !undef; i++ ){ + (void) sprintf( buff, "mpi%d", i + 1 ); + (new->mxpow_i)[ i ] = astReadInt( channel, buff, INT_MAX ); + if( (new->mxpow_i)[ i ] == INT_MAX ) undef = 1; + } + +/* Get the number of coefficients associated with each input of the inverse + transformation. */ + for( i = 0; i < nin && !undef; i++ ){ + (void) sprintf( buff, "nci%d", i + 1 ); + (new->ncoeff_i)[ i ] = astReadInt( channel, buff, INT_MAX ); + if( (new->ncoeff_i)[ i ] == INT_MAX ) undef = 1; + } + +/* Get the coefficient values used by the inverse transformation. */ + iv = 0; + for( i = 0; i < nin && !undef; i++ ){ + + (new->coeff_i)[ i ] = astMalloc( sizeof( double )* + (size_t) new->ncoeff_i[ i ] ); + if( astOK ) { + for( j = 0; j < new->ncoeff_i[ i ]; j++ ){ + (void) sprintf( buff, "ci%d", ++iv ); + (new->coeff_i)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD ); + if( (new->coeff_i)[ i ][ j ] == AST__BAD ) { + (void) sprintf( buff, "ci%d_%d", i + 1, j + 1 ); + (new->coeff_i)[ i ][ j ] = astReadDouble( channel, buff, AST__BAD ); + } + } + } + } + +/* Get the output axis powers associated with each coefficient of the inverse + transformation. */ + iv = 0; + for( i = 0; i < nin && !undef; i++ ){ + (new->power_i)[ i ] = astMalloc( sizeof( int * )* + (size_t) new->ncoeff_i[ i ] ); + if( astOK ) { + for( j = 0; j < new->ncoeff_i[ i ]; j++ ){ + (new->power_i)[ i ][ j ] = astMalloc( sizeof( int )* (size_t) nout ); + if( astOK ) { + for( k = 0; k < nout; k++ ){ + (void) sprintf( buff, "pi%d", ++iv ); + (new->power_i)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 ); + if( (new->power_i)[ i ][ j ][ k ] == 0 ) { + (void) sprintf( buff, "pi%d_%d_%d", i + 1, j + 1, k + 1 ); + (new->power_i)[ i ][ j ][ k ] = astReadInt( channel, buff, 0 ); + } + } + } + } + } + } + +/* Free the arrays if the inverse transformation is undefined. */ + if( undef ) { + new->ncoeff_i = astFree( new->ncoeff_i ); + new->mxpow_i = astFree( new->mxpow_i ); + new->power_i = astFree( new->power_i ); + new->coeff_i = astFree( new->coeff_i ); + } + } + +/* Whether to use an iterative inverse transformation. */ + new->iterinverse = astReadInt( channel, "iterinv", -INT_MAX ); + if ( TestIterInverse( new, status ) ) SetIterInverse( new, new->iterinverse, status ); + +/* Max number of iterations for iterative inverse transformation. */ + new->niterinverse = astReadInt( channel, "niterinv", -INT_MAX ); + if ( TestNiterInverse( new, status ) ) SetNiterInverse( new, new->niterinverse, status ); + +/* Target relative error for iterative inverse transformation. */ + new->tolinverse = astReadDouble( channel, "tolinv", AST__BAD ); + if ( TestTolInverse( new, status ) ) SetTolInverse( new, new->tolinverse, status ); + +/* The Jacobian of the PolyMap's forward transformation has not yet been + found. */ + new->jacobian = NULL; + +/* The linear truncation of the PolyMap has not yet been found. */ + new->lintrunc = NULL; + +/* If an error occurred, clean up by deleting the new PolyMap. */ + if ( !astOK ) new = astDelete( new ); + } + +/* Return the new PolyMap 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. */ + +void astPolyPowers_( AstPolyMap *this, double **work, int ncoord, + const int *mxpow, double **ptr, int point, int fwd, + int *status ){ + if ( !astOK ) return; + (**astMEMBER(this,PolyMap,PolyPowers))( this, work, ncoord, mxpow, ptr, + point, fwd, status ); +} + +AstPolyMap *astPolyTran_( AstPolyMap *this, int forward, double acc, + double maxacc, int maxorder, const double *lbnd, + const double *ubnd, int *status ){ + if ( !astOK ) return NULL; + return (**astMEMBER(this,PolyMap,PolyTran))( this, forward, acc, + maxacc, maxorder, lbnd, + ubnd, status ); +} + +void astPolyCoeffs_( AstPolyMap *this, int forward, int nel, double *array, + int *ncoeff, int *status ){ + if ( !astOK ) return; + (**astMEMBER(this,PolyMap,PolyCoeffs))( this, forward, nel, + array, ncoeff, status ); +} + +void astFitPoly1DInit_( AstPolyMap *this, int forward, double **table, + AstMinPackData *data, double *scales, + int *status ){ + if ( !astOK ) return; + (**astMEMBER(this,PolyMap,FitPoly1DInit))( this, forward, table, data, scales, + status ); +} + + +void astFitPoly2DInit_( AstPolyMap *this, int forward, double **table, + AstMinPackData *data, double *scales, + int *status ){ + if ( !astOK ) return; + (**astMEMBER(this,PolyMap,FitPoly2DInit))( this, forward, table, data, scales, + status ); +} + + + + |