Index: trunk/psastro/src/psastroModelAdjust.c
===================================================================
--- trunk/psastro/src/psastroModelAdjust.c	(revision 15891)
+++ trunk/psastro/src/psastroModelAdjust.c	(revision 19515)
@@ -1,4 +1,7 @@
 # include "psastroStandAlone.h"
 # define NONLIN_TOL 0.001 /* tolerance in pixels */
+# define DEBUG 0
+
+bool psastroModelAdjustBoresite (pmFPAfile *output, pmChip *refChip);
 
 bool psastroModelAdjust (pmConfig *config) {
@@ -13,4 +16,18 @@
     }
 
+    // if we have not measured the boresite position, no adjustment is needed
+    bool fitBoresite = psMetadataLookupBool (&status, recipe, "PSASTRO.MODEL.FIT.BORESITE");
+    if (!status) psAbort ("Can't find recipe option PSASTRO.MODEL.FIT.BORESITE");
+
+    // as an alternative to fit the boresite from a rotation sequence, we can set the boresite
+    // relative to the reference chip coordinates
+    bool setBoresite = psMetadataLookupBool (&status, recipe, "PSASTRO.MODEL.SET.BORESITE");
+    if (!status) psAbort ("Can't find recipe option PSASTRO.MODEL.SET.BORESITE");
+
+    if (fitBoresite && setBoresite) {
+	psError(PS_ERR_IO, true, "invalid to choose both FIT.BORESITE and SET.BORESITE"); 
+	return false; 
+    } 
+	
     pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSASTRO.OUT.MODEL");
     if (!status) psAbort ("Can't find output pmFPAfile PSASTRO.OUT.MODEL");
@@ -29,7 +46,91 @@
     if (!refChip->toFPA) psAbort ("invalid astrometry for reference chip");
 
-    psPlane *boreCH = psPlaneAlloc();
-    psPlane *boreFP = psPlaneAlloc();    
-    psPlane *boreTP = psPlaneAlloc();    
+    // save the TPA region for tranformation inversions below
+    // psRegion *tpaRegion = pmAstromFPInTP (output->fpa);
+    psRegion *fpaRegion = pmAstromFPAExtent (output->fpa);
+
+    if (DEBUG) psastroDumpCorners ("corners.up.raw.dat", "corners.dn.raw.dat", output->fpa);
+
+    if (setBoresite) {
+	float boreXchip = psMetadataLookupF32 (&status, recipe, "PSASTRO.MODEL.BORESITE.X");
+	float boreYchip = psMetadataLookupF32 (&status, recipe, "PSASTRO.MODEL.BORESITE.Y");
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.X0", PS_META_REPLACE, "boresite parameter", boreXchip); 
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.Y0", PS_META_REPLACE, "boresite parameter", boreYchip); 
+    }
+
+    if (fitBoresite || setBoresite) {
+	psastroModelAdjustBoresite (output, refChip);
+    } else {
+	// FPA.BORE.X0,Y0 should be 0,0 in the focal plane, not the chip.  Ask for the
+	// coordinates which make refChip->toFPA(x,y) = (0,0)
+	psPlane *PT = psPlaneTransformGetCenter (refChip->toFPA, NONLIN_TOL);
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.X0", PS_META_REPLACE, "boresite parameter", PT->x); 
+	psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.Y0", PS_META_REPLACE, "boresite parameter", PT->y); 
+	psFree (PT);
+    }
+
+    // rotate the chip-to-FPA transforms to have 0.0 posangle for refChip; 
+    // compensate by rotating fpa to TPA transform
+
+    // get the current posangle of the ref chip
+    float chipAngle = atan2 (refChip->toFPA->y->coeff[1][0], refChip->toFPA->x->coeff[1][0]);
+    fprintf (stderr, "chipAngle: %f\n", chipAngle*PS_DEG_RAD);
+    // psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POSANGLE", PS_META_REPLACE, "boresite parameter", posangle);
+
+    // rotate the chip transforms
+    for (int i = 0; i < output->fpa->chips->n; i++) {
+	pmChip *chip = output->fpa->chips->data[i];
+	if (!chip->toFPA) continue;
+	// skip chips without astrometry
+
+	// save the region of this chip for the inversion below
+	psRegion *region = pmChipPixels (chip);
+
+	psPlaneTransform *toFPA = psPlaneTransformRotate (NULL, chip->toFPA, chipAngle);
+	psFree (chip->toFPA);
+	chip->toFPA = toFPA;
+
+	// invert the new fromFPA transform to get the new toFPA transform
+	psPlaneTransform *fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50);
+	psFree (chip->fromFPA);
+	chip->fromFPA = fromFPA;
+
+	psFree (region);
+
+	// save the transformation in the header
+	pmAstromWriteBilevelChip (chip->hdu->header, chip, NONLIN_TOL);
+    }
+
+    // get the current posangle of the fpa
+    float fpaAngle = atan2 (output->fpa->toTPA->y->coeff[1][0], output->fpa->toTPA->x->coeff[1][0]);
+    fprintf (stderr, "fpaAngle: %f\n", fpaAngle*PS_DEG_RAD);
+    // psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POSANGLE", PS_META_REPLACE, "boresite parameter", posangle);
+
+    // remove the fpa rotation to generate a rotation-free model
+    psPlaneTransform *toTPA = psPlaneTransformRotate (NULL, output->fpa->toTPA, fpaAngle);
+    psFree (output->fpa->toTPA);
+    output->fpa->toTPA = toTPA;
+
+    psFree (output->fpa->fromTPA);
+    output->fpa->fromTPA = psPlaneTransformInvert(NULL, output->fpa->toTPA, *fpaRegion, 50);
+
+    // the model now describes the unrotated focal-plane
+    if (DEBUG) psastroDumpCorners ("corners.up.rot.dat", "corners.dn.rot.dat", output->fpa);
+
+    psMetadata *header = output->fpa->hdu->header;
+    pmAstromWriteBilevelMosaic (header, output->fpa, NONLIN_TOL);
+
+    psFree (fpaRegion);
+
+    return true;
+}
+
+bool psastroModelAdjustBoresite (pmFPAfile *output, pmChip *refChip) {
+
+    bool status;
+
+    psPlane  *boreCH  = psPlaneAlloc();
+    psPlane  *boreFP  = psPlaneAlloc();    
+    psPlane  *boreTP  = psPlaneAlloc();    
     psSphere *boreSky = psSphereAlloc();    
 
@@ -46,91 +147,99 @@
 	// skip the chips without astrometry
 
-	psPlaneTransform *fromFPA = psPlaneTransformSetCenter (NULL, chip->fromFPA, boreFP->x, boreFP->y);
+	// save the FPA region of this chip for the inversion below
+	psRegion *region = pmChipPixels (chip);
+
+	// the current toFPA returns boreFP->x,y for the boresite; subtract this from the transformations
+	chip->toFPA->x->coeff[0][0] -= boreFP->x;
+	chip->toFPA->y->coeff[0][0] -= boreFP->y;
+
+	// psPlaneTransform *toFPA = psPlaneTransformSetCenter (NULL, chip->toFPA, -boreFP->x, -boreFP->y);
+	// psFree (chip->toFPA);
+	// chip->toFPA = toFPA;
+
+	// invert the new fromFPA transform to get the new toFPA transform
+	// the region used here is the region covered by the chip in the FPA
+	psPlaneTransform *fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50); 
 	psFree (chip->fromFPA);
 	chip->fromFPA = fromFPA;
 
-	// invert the new fromFPA transform to get the new toFPA transform
-	psRegion *region = pmChipPixels (chip);
-	psFree (chip->toFPA);
-	chip->toFPA = psPlaneTransformInvert(NULL, chip->fromFPA, *region, 50);
 	psFree (region);
     }
 
-    // XXX we have now shifted the (0,0) pixel of the focal plane to the true boresite from the
-    // reported boresite.  At this point, we need to shift the tangent plane to use the new
-    // center as well.  if the toTPA transform were not linear, we would need to modify fromFPA
-    // to yield 0,0 at the new boresite location (ie, find Po,Qo = toTPA(Lo,Mo), probably could modify the
-    // toTPA/fromTPA transforms to use the new ref pixel, but this would only give us a tranf
-
-    // save the old (L,M) = (0,0) TP coordinate
-    float Po = output->fpa->toTPA->x->coeff[0][0];
-    float Qo = output->fpa->toTPA->y->coeff[0][0];
-
-    // the new toTPA yields the same TP coordinates for FP coordinates offset by Lo,Mo
-    psPlaneTransform *toTPA = psPlaneTransformSetCenter (NULL, output->fpa->toTPA, boreFP->x, boreFP->y);
-    psFree (output->fpa->toTPA);
-    output->fpa->toTPA = toTPA;
-    
-    // we are going to shift P,Q to have toTPA(0,0) = Po, Qo.  
-    // find the sky coordinates of the 0,0 pixel for the new frame
-    boreTP->x = output->fpa->toTPA->x->coeff[0][0] - Po;
-    boreTP->y = output->fpa->toTPA->y->coeff[0][0] - Qo;
+    if (DEBUG) psastroDumpCorners ("corners.up.shf.dat", "corners.dn.shf.dat", output->fpa);
+
+    // we have now adjusted the chips to use the correct boresite position as the center of the focal-plane system.
+    // we now need to reconstruct the TP to FP transformation, starting from stars projected about this new boresite position.
+
+    // find the R,D of the new boresite (boreFP -> 0,0; 0,0 -> -boreFP)
+    boreFP->x = -boreFP->x;
+    boreFP->y = -boreFP->y;
+    psPlaneTransformApply (boreTP, output->fpa->toTPA, boreFP);
     psDeproject (boreSky, boreTP, output->fpa->toSky); // find the RA,DEC coord of the focal-plane coordinate
 
-    // adjust the new TP frame to yield the same old (L,M) = 0,0 TP coordinate:
-    output->fpa->toTPA->x->coeff[0][0] = Po;
-    output->fpa->toTPA->y->coeff[0][0] = Qo;
-
-    // set the projection to use the new (P,Q) = (0,0) position
-    output->fpa->toSky->R = boreSky->r;
-    output->fpa->toSky->D = boreSky->d;
-
+    psProjection *newSky = psProjectionAlloc (boreSky->r, boreSky->d, output->fpa->toSky->Xs, output->fpa->toSky->Ys, output->fpa->toSky->type);
+
+    // generate a collection of points on the sky using the old toTPA transformation and toSky projection, projected with the newSky projection
+    // this is the FPA coordinate range covered by the FP: 
     psRegion *fpaRegion = pmAstromFPAExtent (output->fpa);
-
-    psFree (output->fpa->fromTPA);
-    output->fpa->fromTPA = psPlaneTransformInvert(NULL, output->fpa->toTPA, *fpaRegion, 50);
-
-    // rotate the chip to FPA transforms to have 0.0 posangle; 
-    // compensate by rotating fpa to TPA transforms
-
-    // get the current posangle of the ref chip
-    float posangle = atan2 (refChip->toFPA->y->coeff[1][0], refChip->toFPA->x->coeff[1][0]);
-    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POSANGLE", PS_META_REPLACE, "boresite parameter", posangle);
-
-    // rotate the chip transforms
-    for (int i = 0; i < output->fpa->chips->n; i++) {
-	pmChip *chip = output->fpa->chips->data[i];
-	if (!chip->toFPA) continue;
-	// skip chips without astrometry
-
-	psPlaneTransform *toFPA = psPlaneTransformRotate (NULL, chip->toFPA, posangle);
-	psFree (chip->toFPA);
-	chip->toFPA = toFPA;
-
-	// invert the new fromFPA transform to get the new toFPA transform
-	psRegion *region = pmChipPixels (chip);
-	psFree (chip->fromFPA);
-	chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50);
-	psFree (region);
-
-	// XXX temporary output dump
-	psMetadata *header = chip->hdu->header;
-	// XXX to make the output single-level, this needs to be in a loop *after* the fromTPA rotation below
-	// pmAstromWriteWCS (header, output->fpa, chip, NONLIN_TOL);
-	pmAstromWriteBilevelChip (header, chip, NONLIN_TOL);
-    }
-
-    psPlaneTransform *fromTPA = psPlaneTransformRotate (NULL, output->fpa->fromTPA, posangle);
-    psFree (output->fpa->fromTPA);
-    output->fpa->fromTPA = fromTPA;
-
-    psFree (output->fpa->toTPA);
-    output->fpa->toTPA = psPlaneTransformInvert(NULL, output->fpa->fromTPA, *fpaRegion, 50);
-
-    psMetadata *header = output->fpa->hdu->header;
-    pmAstromWriteBilevelMosaic (header, output->fpa, NONLIN_TOL);
-
+    float dx = (fpaRegion->x1 - fpaRegion->x0) / 50.0;
+    float dy = (fpaRegion->y1 - fpaRegion->y0) / 50.0;
+
+    psPlane fp, tp;
+    psSphere sky;
+
+    psVector *FPx = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *FPy = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *TPx = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *TPy = psVectorAllocEmpty (100, PS_TYPE_F32);
+
+    // XXX a test: boreFP->x,y, should transform to tp.x,y = 0,0
+    fp.x = boreFP->x;
+    fp.y = boreFP->y;
+    psPlaneTransformApply (&tp, output->fpa->toTPA, &fp);
+    psDeproject (&sky, &tp, output->fpa->toSky); // find the RA,DEC coord of the focal-plane coordinate
+    psProject (&tp, &sky, newSky); // find the RA,DEC coord of the focal-plane coordinate
+
+    int Npts = 0;
+    for (fp.x = fpaRegion->x0; fp.x <= fpaRegion->x1; fp.x += dx) {
+	for (fp.y = fpaRegion->y0; fp.y <= fpaRegion->y1; fp.y += dy) {
+	    psPlaneTransformApply (&tp, output->fpa->toTPA, &fp);
+	    psDeproject (&sky, &tp, output->fpa->toSky); // find the RA,DEC coord of the focal-plane coordinate
+	    psProject (&tp, &sky, newSky); // find the RA,DEC coord of the focal-plane coordinate
+
+	    // we are fitting points in the NEW FP system to points in the NEW TP system
+	    FPx->data.F32[Npts] = fp.x - boreFP->x;
+	    FPy->data.F32[Npts] = fp.y - boreFP->y;
+	    TPx->data.F32[Npts] = tp.x;
+	    TPy->data.F32[Npts] = tp.y;
+	    psVectorExtend (FPx, 100, 1);
+	    psVectorExtend (FPy, 100, 1);
+	    psVectorExtend (TPx, 100, 1);
+	    psVectorExtend (TPy, 100, 1);
+	    Npts ++;
+	}
+    }
     psFree (fpaRegion);
 
+    // fit both up and down transformations to the same points
+    psVectorFitPolynomial2D (output->fpa->toTPA->x, NULL, 0, TPx, NULL, FPx, FPy);
+    psVectorFitPolynomial2D (output->fpa->toTPA->y, NULL, 0, TPy, NULL, FPx, FPy);
+    psVectorFitPolynomial2D (output->fpa->fromTPA->x, NULL, 0, FPx, NULL, TPx, TPy);
+    psVectorFitPolynomial2D (output->fpa->fromTPA->y, NULL, 0, FPy, NULL, TPx, TPy);
+
+    psFree (output->fpa->toSky);
+    output->fpa->toSky = newSky;
+
+    if (DEBUG) psastroDumpCorners ("corners.up.bore.dat", "corners.dn.bore.dat", output->fpa);
+
+    psFree (FPx);
+    psFree (FPy);
+    psFree (TPx);
+    psFree (TPy);
+
+    psFree (boreCH);
+    psFree (boreFP);
+    psFree (boreTP);
+    psFree (boreSky);
     return true;
 }