Index: trunk/psastro/src/psastroModelAnalysis.c
===================================================================
--- trunk/psastro/src/psastroModelAnalysis.c	(revision 15887)
+++ trunk/psastro/src/psastroModelAnalysis.c	(revision 15891)
@@ -12,4 +12,7 @@
     }
 
+    pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSASTRO.OUT.MODEL");
+    if (!status) psAbort ("Can't find output pmFPAfile PSASTRO.OUT.MODEL");
+
     // physical pixel scale in microns per pixel
     char *refChip = psMetadataLookupStr (&status, recipe, "PSASTRO.MODEL.REF.CHIP");
@@ -19,22 +22,18 @@
     } 
 
-    /* model analysis steps:
-
-    1) determine POS_ZERO:
-    * choose a reference chip
-    * what is posangle of chip
-    * posangle = FPA.POSANGLE - POS_ZERO
-    * POS_ZERO = FPA.POSANGLE - posangle
-
-    ** average over input images
-
-    2) determine boresite model
-
-    * find CRPIX1,2 for a series of images
-    * fit ellipse 
-
-    */
-
-    /***** find the POSANGLE offset *****/
+    /* model analysis:
+     *
+     * determine POS_ZERO via comparison of measured and reported posangles
+     * POS_ZERO = FPA.POSANGLE - posangle
+     *
+     * determine boresite model:
+     * X = Xo + R_X cos(FPA.POSANGLE - T_0) cos(P_0) + R_Y sin(FPA.POSANGLE - T_0) sin(P_0) 
+     * Y = Yo + R_Y sin(FPA.POSANGLE - T_0) cos(P_0) - R_X cos(FPA.POSANGLE - T_0) sin(P_0) 
+     * position of reported boresite in reference chip pixels
+     * Xo, Yo : true coordinate of boresite (rotator center) in reference chip pixels
+     * R_X, R_Y : amplitude of boresite offset
+     * T_0 : reference angle for rotator
+     * P_0 : orientation of boresite ellipse
+     */
 
     // select the input pmFPAfile pointers
@@ -44,25 +43,39 @@
     psArray *files = psListToArray (item->data.list);
 
-    // temp data vectors
+    // data storage vectors for measurements
     psVector *posZero  = psVectorAlloc (files->n, PS_TYPE_F32);
     psVector *Po       = psVectorAlloc (files->n, PS_TYPE_F32);
-    psVector *Lo       = psVectorAlloc (files->n, PS_TYPE_F32);
-    psVector *Mo       = psVectorAlloc (files->n, PS_TYPE_F32);
+    psVector *Xo       = psVectorAlloc (files->n, PS_TYPE_F32);
+    psVector *Yo       = psVectorAlloc (files->n, PS_TYPE_F32);
 
+    // counter for accepted measured values
     int n = 0;
 
-    // convert the headers for the input file into fpa astrometry terms
+    // Re-select the output fpa.  The output fpa needs to have valid astrometry for the
+    // refChip.  output->fpa is a copy of the pointer to one of the input->fpa, but the choice
+    // is arbitrary.  select a new one that has an existing ref chip
+    psFree (output->fpa);
+    output->fpa = NULL;
+
+    // extract the relevant measured and reported values from the reference chip
     for (int i = 0; i < files->n; i++) {
 	psMetadataItem *file = files->data[i];
 	pmFPAfile *input = file->data.V;
 
+	// reported rotator position angle (this should perhaps be ROT, not POS)
 	double POSANGLE = psMetadataLookupF64 (&status, input->fpa->concepts, "FPA.POSANGLE"); 
 	if (!status) psAbort ("missing FPA.POSANGLE");
 
-    	// get chip from name
+    	// get reference chip from name
 	pmChip *chip = pmConceptsChipFromName (input->fpa, refChip);
 	if (!chip) psAbort ("invalid chip name for reference");
 
+	fprintf (stderr, "input %d : %zx : %zx : %zx\n", i, (size_t) input, (size_t) chip, (size_t) chip->toFPA);
 	if (!chip->toFPA) continue;
+
+	if (!output->fpa) {
+	    // this one matches
+	    output->fpa = psMemIncrRefCounter(input->fpa);
+	}
 
 	// we have two measurements of the posangle (may be parity flipped to different quadrants)
@@ -72,13 +85,13 @@
 	while (posZero->data.F32[n] <   0.0) posZero->data.F32[n] += 360.0;
 
-	Po->data.F32[n] = POSANGLE * PM_RAD_DEG;
-	Lo->data.F32[n] = chip->toFPA->x->coeff[0][0];
-	Mo->data.F32[n] = chip->toFPA->y->coeff[0][0];
-	fprintf (stderr, "%d : %f %f : %f = %f - %f\n", i, Lo->data.F32[n], Mo->data.F32[n], posZero->data.F32[n], POSANGLE, posangle);
+	Po->data.F32[n] = POSANGLE * PM_RAD_DEG; // reported position angle
+	Xo->data.F32[n] = chip->fromFPA->x->coeff[0][0]; // reported boresite x position in ref chip coordinates
+	Yo->data.F32[n] = chip->fromFPA->y->coeff[0][0]; // reported boresite y position in ref chip coordinates
+	fprintf (stderr, "%d : %f %f : %f = %f - %f\n", i, Xo->data.F32[n], Yo->data.F32[n], posZero->data.F32[n], POSANGLE, posangle);
 	n ++;
     }
       
-    Lo->n = n;
-    Mo->n = n;
+    Xo->n = n;
+    Yo->n = n;
     Po->n = n;
     posZero->n = n;
@@ -88,6 +101,16 @@
 
     fprintf (stderr, "pos zero %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POS_ZERO", PS_META_REPLACE, "offset between obs and meas posangle", stats->sampleMedian); 
 
-    psastroModelFitBoresite (Lo, Mo, Po);
+    psVector *params = psastroModelFitBoresite (Xo, Yo, Po);
+    if (params->n != 6) psAbort ("error");
+
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.X0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_X0]); 
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.Y0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_Y0]); 
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.RX", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_RX]); 
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.RY", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_RY]); 
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.T0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_T0]); 
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.P0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_P0]); 
+    psMetadataAddStr (output->fpa->concepts, PS_LIST_TAIL, "FPA.REF.CHIP", PS_META_REPLACE, "boresite parameter", refChip);
 
     return true;