Index: /trunk/psLib/src/astronomy/psAstrometry.c
===================================================================
--- /trunk/psLib/src/astronomy/psAstrometry.c	(revision 3496)
+++ /trunk/psLib/src/astronomy/psAstrometry.c	(revision 3497)
@@ -8,6 +8,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-02-17 19:26:23 $
+ *  @version $Revision: 1.60 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-03-24 22:36:16 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -26,4 +26,5 @@
 #include "psConstants.h"
 #include "psAstronomyErrors.h"
+#include "psMatrix.h"
 #include "psTrace.h"
 #include "psLogMsg.h"
@@ -106,4 +107,7 @@
              (X2 * -((E/B) / (F - ((C*E)/B)))) +
              (Y2 * (1.0 / (F - ((C*E)/B))));
+ 
+XXX: Since thre is now a general psPlaneTransformInvertTmp() function, we
+should rename this.
  *****************************************************************************/
 static psPlaneTransform *invertPlaneTransform(psPlaneTransform *transform)
@@ -908,2 +912,552 @@
 }
 
+
+
+
+
+/*****************************************************************************
+psPlaneTransformInvertHEY(out, in, region, nSamples): this is an earlier
+version of this function which doesnot separate the code for building and
+solving the matrixequations.
+ 
+XXX: Delete this codeonce you know the other version works.
+ *****************************************************************************/
+psPlaneTransform *psPlaneTransformInvertHEY(psPlaneTransform *out,
+        const psPlaneTransform *in,
+        psRegion *region,
+        int nSamples)
+{
+    PS_PTR_CHECK_NULL(in, NULL);
+    if (isProjectionLinear((psPlaneTransform *) in)) {
+        return(invertPlaneTransform((psPlaneTransform *) in));
+    }
+    psPlaneTransform *myPT = NULL;
+    psPlane *inCoord = psPlaneAlloc();
+    psPlane *outCoord = psPlaneAlloc();
+
+    PS_PTR_CHECK_NULL(region, NULL);
+    PS_INT_COMPARE(0, nSamples, NULL);
+
+    // XXX: Is this correct?
+    psS32 order = PS_MAX(in->x->nX, in->x->nY);
+
+    //
+    // Allocate a new psPlaneTransform if "out" is NULL, or has the wrong size.
+    //
+    if (out == NULL) {
+        myPT = psPlaneTransformAlloc(order, order);
+    } else {
+        if (!((out->x->nX == order) &&
+                (out->x->nY == order) &&
+                (out->y->nX == order) &&
+                (out->y->nY == order))) {
+            psFree(out);
+            myPT = psPlaneTransformAlloc(order, order);
+        } else {
+            myPT = out;
+        }
+    }
+    // XXX: Initialize myPT?
+
+    //
+    // Create fake polynomial to use in evaluation
+    //
+    psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
+    for (int i = 0; i < order; i++) {
+        for (int j = 0; j < order; j++) {
+            fakePoly->coeff[i][j] = 1.0; // Set all coeffecients to 1
+            fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate
+        }
+    }
+
+    //
+    // Create a grid of xin,yin --> xout,yout
+    //
+    psVector *xIn = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
+    psVector *yIn = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
+    psVector *xOut = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
+    psVector *yOut = psVectorAlloc(nSamples * nSamples, PS_TYPE_F32);
+
+    //
+    // Initialize the grid of points
+    //
+    for (int yint = 0; yint < nSamples; yint++) {
+        inCoord->y = region->y0 + ((psF32) yint) * ((region->y1 - region->y0) / ((psF32) nSamples));
+        for (int xint = 0; xint < nSamples; xint++) {
+            inCoord->x = region->x0 + ((psF32) xint) * ((region->x1 - region->x0) / ((psF32) nSamples));
+
+            (void)psPlaneTransformApply(outCoord, in, inCoord);
+            xOut->data.F32[yint*nSamples + xint] = inCoord->x;
+            yOut->data.F32[yint*nSamples + xint] = inCoord->y;
+            xIn->data.F32[yint*nSamples + xint] = outCoord->x;
+            yIn->data.F32[yint*nSamples + xint] = outCoord->y;
+        }
+    }
+
+    //
+    // Initialise the matrix and vectors
+    //
+    psS32 nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
+    psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
+    psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
+    psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
+    for (psS32 i = 0; i < nCoeff; i++) {
+        for (psS32 j = 0; j < nCoeff; j++) {
+            matrix->data.F64[i][j] = 0.0;
+        }
+        xVector->data.F64[i] = 0.0;
+        yVector->data.F64[i] = 0.0;
+    }
+
+    //
+    // Iterate over the grid points
+    //
+    for (psS32 g = 0; g < nSamples*nSamples; g++) {
+        // Iterate over the polynomial coefficients, accumulating the matrix and vectors
+
+        for (psS32 i = 0, ijIndex = 0; i < order; i++) {
+            for (psS32 j = 0; j < order - i; j++, ijIndex++) {
+
+                fakePoly->mask[i][j] = 0;
+                psF64 ijPoly = psDPolynomial2DEval(fakePoly, (psF64) xIn->data.F32[g], (psF64) yIn->data.F32[g]);
+                fakePoly->mask[i][j] = 1;
+
+                for (psS32 m = 0, mnIndex = 0; m < order; m++) {
+                    for (psS32 n = 0; n < order - m; n++, mnIndex++) {
+                        fakePoly->mask[m][n] = 0;
+                        psF64 mnPoly = psDPolynomial2DEval(fakePoly, (psF64) xIn->data.F32[g], (psF64) yIn->data.F32[g]);
+                        fakePoly->mask[m][n] = 1;
+
+                        matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
+                    }
+                }
+
+                xVector->data.F64[ijIndex] += ijPoly * (psF64)xOut->data.F32[g];
+                yVector->data.F64[ijIndex] += ijPoly * (psF64)yOut->data.F32[g];
+            }
+        }
+    }
+
+    //
+    // Solution via LU Decomposition
+    //
+    psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
+    psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
+    psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
+    psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
+
+    //
+    // Stuff coefficients into transformation
+    //
+    for (psS32 i = 0, ijIndex = 0; i < order; i++) {
+        for (psS32 j = 0; j < order - i; j++, ijIndex++) {
+            myPT->x->coeff[i][j] = xSolution->data.F64[ijIndex];
+            myPT->y->coeff[i][j] = ySolution->data.F64[ijIndex];
+        }
+    }
+
+    return(myPT);
+}
+
+
+/*****************************************************************************
+multiplyCoeffs(trans1, trans2): Takes two 2-D polynomials as input and
+multiplies them.  Basically, for each non-zero coeff in the trans1 coeff[][]
+array, you must multiply by all non-zero coeffs in trans2.
+ 
+XXX: Inefficient in that the out polynomial is allocated every time.
+ *****************************************************************************/
+psDPolynomial2D *multiplyDPoly2D(psDPolynomial2D *trans1,
+                                 psDPolynomial2D *trans2)
+{
+    psS32 orderX = (trans1->nX + trans2->nX) - 1;
+    psS32 orderY = (trans1->nX + trans2->nX) - 1;
+
+    psDPolynomial2D *out = psDPolynomial2DAlloc(orderX, orderY, PS_POLYNOMIAL_ORD);
+    for (psS32 i = 0 ; i < out->nX; i++) {
+        for (psS32 j = 0 ; j < out->nY; j++) {
+            out->coeff[i][j] = 0.0;
+            out->mask[i][j] = 0;
+        }
+    }
+
+    for (psS32 t1x = 0 ; t1x < trans1->nX ; t1x++) {
+        for (psS32 t1y = 0 ; t1y < trans1->nY ; t1y++) {
+            if (0.0 != trans1->coeff[t1x][t1y]) {
+                for (psS32 t2x = 0 ; t2x < trans2->nX ; t2x++) {
+                    for (psS32 t2y = 0 ; t2y < trans2->nY ; t2y++) {
+                        out->coeff[t1x+t2x][t1y+t2y]+= (trans1->coeff[t1x][t1y] * trans2->coeff[t2x][t2y]);
+                    }
+                }
+            }
+        }
+    }
+    return(out);
+}
+
+
+
+
+/*****************************************************************************
+psPlaneTransformCombineTmp(out, trans1, trans2)
+ 
+XXX: Much room for optimization.  Currently, we call the polyMultiply
+routine far too many times.
+ *****************************************************************************/
+psPlaneTransform *psPlaneTransformCombineTmp(psPlaneTransform *out,
+        const psPlaneTransform *trans1,
+        const psPlaneTransform *trans2)
+{
+    PS_PTR_CHECK_NULL(trans1, NULL);
+    PS_PTR_CHECK_NULL(trans2, NULL);
+
+    //
+    // Determine the size of the new psPlaneTransform.
+    //
+    // PS_MAX(  Number of x terms in T2->x * number of x terms in T1->x,
+    //          Number of y terms in T2->x * number of x terms in T1->y,
+    psS32 orderXnX = PS_MAX((trans2->x->nX * trans1->x->nX),
+                            (trans2->x->nY * trans1->y->nX));
+    psS32 orderXnY = PS_MAX((trans2->x->nX * trans1->x->nY),
+                            (trans2->x->nY * trans1->y->nY));
+
+    psS32 orderYnX = PS_MAX((trans2->y->nX * trans1->x->nX),
+                            (trans2->y->nY * trans1->y->nX));
+    psS32 orderYnY = PS_MAX((trans2->y->nX * trans1->x->nY),
+                            (trans2->y->nY * trans1->y->nY));
+    psS32 orderX = PS_MAX(orderXnX, orderYnX);
+    psS32 orderY = PS_MAX(orderXnY, orderYnY);
+
+    //
+    // Allocate the new psPlaneTransform, if necessary.
+    //
+    psPlaneTransform *myPT = NULL;
+    if (out == NULL) {
+        myPT = psPlaneTransformAlloc(orderX, orderY);
+    } else {
+        if ((out->x->nX == orderX) && (out->x->nY == orderY) &&
+                (out->y->nX == orderX) && (out->y->nY == orderY)) {
+            myPT = out;
+        } else {
+            psFree(out);
+            myPT = psPlaneTransformAlloc(orderX, orderY);
+        }
+    }
+
+    //
+    // Initialize the new psPlaneTransform, if necessary.
+    //
+    for (psS32 i = 0 ; i < orderX ; i++) {
+        for (psS32 j = 0 ; j < orderY ; j++) {
+            myPT->x->coeff[i][j] = 0.0;
+            myPT->x->mask[i][j] = 0;
+            myPT->y->coeff[i][j] = 0.0;
+            myPT->y->mask[i][j] = 0;
+        }
+    }
+
+    //
+    // For each term (a * x^i * y^j) in trans2, we substitute the appropriate
+    // equation from trans1, and raise it to the appropriate power.  This is
+    // done via the multiplyDPoly2D().  The result is a polynomial (currPoly)
+    // and its coefficients are added into the myPT coeff matrix.
+    //
+    // XXX: This is horribly inefficient in that the trans1 polys are repeatedly
+    // multiplied against themselves.  This can easily be improved.
+    //
+    for (psS32 t2x = 0 ; t2x < trans2->x->nX ; t2x++) {
+        for (psS32 t2y = 0 ; t2y < trans2->x->nY ; t2y++) {
+            psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
+            currPoly->coeff[0][0] = 1.0;
+            currPoly->mask[0][0] = 0;
+            psDPolynomial2D *newPoly = NULL;
+
+            if (trans2->x->mask[t2x][t2y] == 0) {
+
+                // Must raise trans1->y to the t2y-power.
+                for (psS32 c = 0 ; c < t2y; c++) {
+                    newPoly = multiplyDPoly2D(currPoly, trans1->y);
+                    psFree(currPoly);
+                    currPoly = newPoly;
+                }
+
+                // Must raise trans1->x to the t2x-power.
+                for (psS32 c = 0 ; c < t2x; c++) {
+                    newPoly = multiplyDPoly2D(currPoly, trans1->x);
+                    psFree(currPoly);
+                    currPoly = newPoly;
+                }
+
+                // Set the appropriate coeffs in myPT->x
+                for (psS32 i = 0 ; i < currPoly->nX ; i++) {
+                    for (psS32 j = 0 ; j < currPoly->nY ; j++) {
+                        myPT->x->coeff[i][j]+= currPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y];
+                    }
+                }
+            }
+            psFree(currPoly);
+        }
+    }
+
+
+
+    for (psS32 t2x = 0 ; t2x < trans2->y->nX ; t2x++) {
+        for (psS32 t2y = 0 ; t2y < trans2->y->nY ; t2y++) {
+            psDPolynomial2D *currPoly = psDPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
+            currPoly->coeff[0][0] = 1.0;
+            currPoly->mask[0][0] = 0;
+            psDPolynomial2D *newPoly = NULL;
+
+            if (trans2->y->mask[t2x][t2y] == 0) {
+
+                // Must raise trans1->y to the t2y-power.
+                for (psS32 c = 0 ; c < t2y; c++) {
+                    newPoly = multiplyDPoly2D(currPoly, trans1->y);
+                    psFree(currPoly);
+                    currPoly = newPoly;
+                }
+
+                // Must raise trans1->x to the t2x-power.
+                for (psS32 c = 0 ; c < t2x; c++) {
+                    newPoly = multiplyDPoly2D(currPoly, trans1->x);
+                    psFree(currPoly);
+                    currPoly = newPoly;
+                }
+
+                // Set the appropriate coeffs in myPT->x
+                for (psS32 i = 0 ; i < currPoly->nX ; i++) {
+                    for (psS32 j = 0 ; j < currPoly->nY ; j++) {
+                        myPT->y->coeff[i][j]+= currPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y];
+                    }
+                }
+            }
+            psFree(currPoly);
+        }
+    }
+
+    return(myPT);
+}
+
+/*****************************************************************************
+psPlaneTranformFitTmp(trans, source, dest, nRejIter, sigmaClip)
+ 
+XXX: What about nRejIter?  Iterations?
+XXX: Use static vectors for internal data.
+ *****************************************************************************/
+bool psPlaneTranformFitTmp(psPlaneTransform *trans,
+                           const psArray *source,
+                           const psArray *dest,
+                           int nRejIter,
+                           float sigmaClip)
+{
+    PS_PTR_CHECK_NULL(trans, NULL);
+    PS_PTR_CHECK_NULL(source, NULL);
+    PS_PTR_CHECK_NULL(dest, NULL);
+
+    // Ensure that the input transformation is symmetrical.
+    if ((trans->x->nX != trans->x->nY) ||
+            (trans->y->nX != trans->y->nY) ||
+            (trans->x->nX != trans->y->nX)) {
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Input transformation must have same nX==nY.");
+    }
+
+    psS32 numCoords = PS_MIN(source->n, dest->n);
+    // This is not really necessary because of above conditionals.
+    psS32 order = PS_MAX(trans->x->nX, trans->x->nY);
+
+    //
+    // Create fake polynomial to use in evaluation
+    //
+    psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
+    for (int i = 0; i < order; i++) {
+        for (int j = 0; j < order; j++) {
+            fakePoly->coeff[i][j] = 1.0;
+            fakePoly->mask[i][j] = 1;       // Mask all coefficients; unmask to evaluate
+        }
+    }
+
+    //
+    // Initialize the matrix and vectors
+    //
+    psS32 nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
+    psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
+    psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
+    psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
+    for (psS32 i = 0; i < nCoeff; i++) {
+        for (psS32 j = 0; j < nCoeff; j++) {
+            matrix->data.F64[i][j] = 0.0;
+        }
+        xVector->data.F64[i] = 0.0;
+        yVector->data.F64[i] = 0.0;
+    }
+
+    //
+    // Iterate over the grid points
+    //
+    for (psS32 g = 0; g < numCoords; g++) {
+        // Iterate over the polynomial coefficients, accumulating the matrix and vectors
+
+        for (psS32 i = 0, ijIndex = 0; i < order; i++) {
+            for (psS32 j = 0; j < order - i; j++, ijIndex++) {
+                fakePoly->mask[i][j] = 0;
+                psF64 xIn = ((psPlane *) source->data[g])->x;
+                psF64 yIn = ((psPlane *) source->data[g])->y;
+                psF64 xOut = ((psPlane *) dest->data[g])->x;
+                psF64 yOut = ((psPlane *) dest->data[g])->y;
+                psF64 ijPoly = psDPolynomial2DEval(fakePoly, xIn, yIn);
+                fakePoly->mask[i][j] = 1;
+
+                for (psS32 m = 0, mnIndex = 0; m < order; m++) {
+                    for (psS32 n = 0; n < order - m; n++, mnIndex++) {
+                        fakePoly->mask[m][n] = 0;
+                        psF64 mnPoly = psDPolynomial2DEval(fakePoly, xIn, yIn);
+                        fakePoly->mask[m][n] = 1;
+
+                        matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
+                    }
+                }
+
+                xVector->data.F64[ijIndex] += ijPoly * xOut;
+                yVector->data.F64[ijIndex] += ijPoly * yOut;
+            }
+        }
+    }
+
+    //
+    // Solution via LU Decomposition
+    //
+    psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
+    psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
+    psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
+    psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
+
+    //
+    // XXX: Should check the output of the matrix routines and return false if bad.
+    //
+
+    //
+    // Stuff coefficients into transformation
+    //
+    for (psS32 i = 0, ijIndex = 0; i < order; i++) {
+        for (psS32 j = 0; j < order - i; j++, ijIndex++) {
+            trans->x->coeff[i][j] = xSolution->data.F64[ijIndex];
+            trans->y->coeff[i][j] = ySolution->data.F64[ijIndex];
+        }
+    }
+
+    psFree(fakePoly);
+    psFree(permutation);
+    psFree(luMatrix);
+    psFree(xSolution);
+    psFree(ySolution);
+    psFree(matrix);
+    psFree(xVector);
+    psFree(yVector);
+
+    return(true);
+}
+
+
+/*****************************************************************************
+psPlaneTransformInvertTmp(out, in, region, nSamples)
+ 
+// XXX: Use static data structures.
+ *****************************************************************************/
+psPlaneTransform *psPlaneTransformInvertTmp(psPlaneTransform *out,
+        const psPlaneTransform *in,
+        psRegion *region,
+        int nSamples)
+{
+    PS_PTR_CHECK_NULL(in, NULL);
+    //
+    // If the transform is linear, then invert it exactly and return.
+    //
+    if (isProjectionLinear((psPlaneTransform *) in)) {
+        return(invertPlaneTransform((psPlaneTransform *) in));
+    }
+    PS_PTR_CHECK_NULL(region, NULL);
+    PS_INT_COMPARE(0, nSamples, NULL);
+
+    // Ensure that the input transformation is symmetrical.
+    if ((in->x->nX != in->x->nY) ||
+            (in->y->nX != in->y->nY) ||
+            (in->x->nX != in->y->nX)) {
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Input transformation must have same nX==nY.");
+    }
+    psS32 order = PS_MAX(in->x->nX, in->x->nY);
+
+    psPlaneTransform *myPT = NULL;
+    psPlane *inCoord = psPlaneAlloc();
+    psPlane *outCoord = psPlaneAlloc();
+
+    //
+    // Allocate a new psPlaneTransform if "out" is NULL, or has the wrong size.
+    //
+    if (out == NULL) {
+        myPT = psPlaneTransformAlloc(order, order);
+    } else {
+        if ((out->x->nX == order) && (out->x->nY == order) &&
+                (out->y->nX == order) && (out->y->nY == order)) {
+            myPT = out;
+        } else {
+            psFree(out);
+            myPT = psPlaneTransformAlloc(order, order);
+        }
+    }
+
+    //
+    // Copy the input transform to myPT.
+    //
+    for (psS32 i = 0 ; i < in->x->nX ; i++) {
+        for (psS32 j = 0 ; j < in->x->nY ; j++) {
+            myPT->x->coeff[i][j] = in->x->coeff[i][j];
+        }
+    }
+    for (psS32 i = 0 ; i < in->y->nX ; i++) {
+        for (psS32 j = 0 ; j < in->y->nY ; j++) {
+            myPT->y->coeff[i][j] = in->y->coeff[i][j];
+        }
+    }
+
+    //
+    // Create a grid of xin,yin --> xout,yout
+    //
+    psArray *inData = psArrayAlloc(nSamples * nSamples);
+    psArray *outData = psArrayAlloc(nSamples * nSamples);
+    for (psS32 i = 0 ; i < inData->n; i++) {
+        inData->data[i] = (psPtr *) psPlaneAlloc();
+        outData->data[i] = (psPtr *) psPlaneAlloc();
+    }
+
+    //
+    // Initialize the grid.
+    //
+    psS32 cnt = 0;
+    for (int yint = 0; yint < nSamples; yint++) {
+        inCoord->y = region->y0 + ((psF32) yint) * ((region->y1 - region->y0) / ((psF32) nSamples));
+        for (int xint = 0; xint < nSamples; xint++) {
+            inCoord->x = region->x0 + ((psF32) xint) * ((region->x1 - region->x0) / ((psF32) nSamples));
+            (void)psPlaneTransformApply(outCoord, in, inCoord);
+
+            ((psPlane *) outData->data[cnt])->x = inCoord->x;
+            ((psPlane *) outData->data[cnt])->y = inCoord->y;
+            ((psPlane *) inData->data[cnt])->x = outCoord->x;
+            ((psPlane *) inData->data[cnt])->y = outCoord->y;
+
+            cnt++;
+        }
+    }
+    bool rc = psPlaneTranformFitTmp(myPT, inData, outData, 10, 100.0);
+
+    psFree(inCoord);
+    psFree(outCoord);
+    psFree(inData);
+    psFree(outData);
+
+    if (rc == true) {
+        return(myPT);
+    }
+
+    // XXX: Generate an error message, or warning message.
+    return(NULL);
+}
Index: /trunk/psLib/src/astronomy/psAstrometry.h
===================================================================
--- /trunk/psLib/src/astronomy/psAstrometry.h	(revision 3496)
+++ /trunk/psLib/src/astronomy/psAstrometry.h	(revision 3497)
@@ -8,6 +8,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-02-17 19:26:23 $
+*  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-03-24 22:36:16 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -534,3 +534,32 @@
 );
 
+
+// XXX: These functions don't belong here.  Will migrate to psCoords.c later.
+// XXX: Doxygenate.
+psPlaneTransform *psPlaneTransformInvertTmp(
+    psPlaneTransform *out,
+    const psPlaneTransform *in,
+    psRegion *region,
+    int nSamples
+);
+
+// XXX: Doxygenate.
+psPlaneTransform *psPlaneTransformCombineTmp(
+    psPlaneTransform *out,
+    const psPlaneTransform *trans1,
+    const psPlaneTransform *trans2
+);
+
+// XXX: Doxygenate.
+bool psPlaneTranformFitTmp(
+    psPlaneTransform *trans,
+    const psArray *source,
+    const psArray *dest,
+    int nRejIter,
+    float sigmaClip
+);
+
+
+
+
 #endif