Index: /trunk/psastro/src/Makefile.am
===================================================================
--- /trunk/psastro/src/Makefile.am	(revision 15886)
+++ /trunk/psastro/src/Makefile.am	(revision 15887)
@@ -27,4 +27,6 @@
 	psastroModelDataLoad.c      \
 	psastroModelAnalysis.c      \
+	psastroModelBoresite.c      \
+	psastroModelFitBoresite.c   \
 	psastroCleanup.c
 
Index: /trunk/psastro/src/psastroModelAnalysis.c
===================================================================
--- /trunk/psastro/src/psastroModelAnalysis.c	(revision 15886)
+++ /trunk/psastro/src/psastroModelAnalysis.c	(revision 15887)
@@ -45,7 +45,8 @@
 
     // temp data vectors
-    psVector *posZero = psVectorAlloc (files->n, PS_TYPE_F32);
-    psVector *Lo      = psVectorAlloc (files->n, PS_TYPE_F32);
-    psVector *Mo      = psVectorAlloc (files->n, PS_TYPE_F32);
+    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);
 
     int n = 0;
@@ -71,4 +72,5 @@
 	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];
@@ -79,4 +81,5 @@
     Lo->n = n;
     Mo->n = n;
+    Po->n = n;
     posZero->n = n;
 
@@ -86,4 +89,6 @@
     fprintf (stderr, "pos zero %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
 
+    psastroModelFitBoresite (Lo, Mo, Po);
+
     return true;
 }
Index: /trunk/psastro/src/psastroModelBoresite.c
===================================================================
--- /trunk/psastro/src/psastroModelBoresite.c	(revision 15887)
+++ /trunk/psastro/src/psastroModelBoresite.c	(revision 15887)
@@ -0,0 +1,101 @@
+# include "psastroStandAlone.h"
+
+// chisq = sum ((L_obs - L_fit(t))^2 + (M_obs - M_fit(t))^2)
+// coord[0] = measured L or measured M
+// coord[1] =          0 or          1
+psF32 psastroModelBoresite (psVector *deriv, const psVector *params, const psVector *coord) {
+
+    psF32 *PAR = params->data.F32;
+
+    float csphi = cos(PAR[PAR_P0]);
+    float snphi = sin(PAR[PAR_P0]);
+
+    float dtheta = coord->data.F32[0] - PAR[PAR_T0];
+    float cstht = cos(dtheta);
+    float sntht = sin(dtheta);
+
+    // value is L
+    if (coord->data.F32[1] == 0) {
+
+	float value = PAR[PAR_L0] + PAR[PAR_RL]*cstht*csphi + PAR[PAR_RM]*sntht*snphi;
+
+	if (deriv) {
+	    psF32 *dPAR = deriv->data.F32;
+	    dPAR[PAR_L0] = 1.0;
+	    dPAR[PAR_M0] = 0.0;
+	    dPAR[PAR_RL] = +cstht*csphi;
+	    dPAR[PAR_RM] = +sntht*snphi;
+	    dPAR[PAR_P0] = -PAR[PAR_RL]*cstht*snphi + PAR[PAR_RM]*sntht*csphi;
+	    dPAR[PAR_T0] =  PAR[PAR_RL]*sntht*csphi - PAR[PAR_RM]*cstht*snphi;
+	}
+	return (value);
+    }  
+
+    // value is M
+    if (coord->data.F32[1] == 1) {
+	float value = PAR[PAR_M0] + PAR[PAR_RM]*sntht*csphi - PAR[PAR_RL]*cstht*snphi;
+
+	if (deriv) {
+	    psF32 *dPAR = deriv->data.F32;
+	    dPAR[PAR_L0]  = 0.0;
+	    dPAR[PAR_M0]  = 1.0;
+	    dPAR[PAR_RL]  = -cstht*snphi;
+	    dPAR[PAR_RM]  = +sntht*csphi;
+	    dPAR[PAR_P0]  = -PAR[PAR_RM]*sntht*snphi - PAR[PAR_RL]*cstht*csphi;
+	    dPAR[PAR_T0]  = -PAR[PAR_RM]*sntht*csphi - PAR[PAR_RL]*cstht*snphi;
+	}
+	return (value);
+    }  
+    psAbort ("programming error: invalid coordinate");
+}
+
+# undef PAR_L0
+# undef PAR_RL
+# undef PAR_PHI
+# define PAR_L0  0
+# define PAR_RL  1
+# define PAR_PHI 2
+
+psF32 psastroModelBoresiteL (psVector *deriv, const psVector *params, const psVector *coord) {
+
+    psF32 *PAR = params->data.F32;
+
+    float csphi = cos(PAR[PAR_PHI]);
+    float snphi = sin(PAR[PAR_PHI]);
+
+    float theta = coord->data.F32[0];
+    float cstht = cos(theta);
+    float sntht = sin(theta);
+
+    float value = PAR[PAR_L0] + PAR[PAR_RL]*cstht*csphi - PAR[PAR_RL]*sntht*snphi;
+
+    if (deriv) {
+	psF32 *dPAR = deriv->data.F32;
+	dPAR[PAR_L0]  = 1.0;
+	dPAR[PAR_RL]  = cstht*csphi - sntht*snphi;
+	dPAR[PAR_PHI] = -PAR[PAR_RL]*cstht*snphi - PAR[PAR_RL]*sntht*csphi;
+    }
+    return (value);
+}
+
+psF32 psastroModelBoresiteM (psVector *deriv, const psVector *params, const psVector *coord) {
+
+    psF32 *PAR = params->data.F32;
+
+    float csphi = cos(PAR[PAR_PHI]);
+    float snphi = sin(PAR[PAR_PHI]);
+
+    float theta = coord->data.F32[0];
+    float cstht = cos(theta);
+    float sntht = sin(theta);
+
+    float value = PAR[PAR_L0] + PAR[PAR_RL]*sntht*csphi + PAR[PAR_RL]*cstht*snphi;
+
+    if (deriv) {
+	psF32 *dPAR = deriv->data.F32;
+	dPAR[PAR_L0]  = 1.0;
+	dPAR[PAR_RL]  = +sntht*csphi + cstht*snphi;
+	dPAR[PAR_PHI] = -PAR[PAR_RL]*sntht*snphi + PAR[PAR_RL]*cstht*csphi;
+    }
+    return (value);
+}
Index: /trunk/psastro/src/psastroModelFitBoresite.c
===================================================================
--- /trunk/psastro/src/psastroModelFitBoresite.c	(revision 15887)
+++ /trunk/psastro/src/psastroModelFitBoresite.c	(revision 15887)
@@ -0,0 +1,210 @@
+# include "psastroStandAlone.h"
+
+// we now have a set of observed L,M values.  fit these to the boresite model
+bool psastroModelFitBoresite (psVector *Lo, psVector *Mo, psVector *Po) {
+
+    assert (Lo->n > 2);
+    assert (Lo->n == Mo->n);
+    assert (Lo->n == Po->n);
+
+    // arrays to hold the data to be fitted
+    psArray *x = psArrayAlloc(2*Lo->n);
+    psVector *y = psVectorAlloc(2*Lo->n, PS_TYPE_F32);
+
+    int n = 0;
+    for (int i = 0; i < Lo->n; i++) {
+
+	psVector *coord = NULL;
+
+	// L coordinate value
+	coord = psVectorAlloc (2, PS_TYPE_F32);
+	coord->data.F32[1] = 0.0;
+	coord->data.F32[0] = Po->data.F32[i];
+	x->data[n] = coord;
+	y->data.F32[n] = Lo->data.F32[i];
+	n++;
+	
+	// M coordinate value
+	coord = psVectorAlloc (2, PS_TYPE_F32);
+	coord->data.F32[1] = 1.0;
+	coord->data.F32[0] = Po->data.F32[i];
+	x->data[n] = coord;
+	y->data.F32[n] = Mo->data.F32[i];
+	n++;
+    }	
+    assert (x->n == n);
+    assert (y->n == n);
+
+    psVector *params = psVectorAlloc (6, PS_TYPE_F32);
+    
+    // create the minimization constraints
+    psMinConstraint *constraint = psMinConstraintAlloc();
+
+    // XXX for now, no parameter masks, skip checkLimits
+    // constraint->checkLimits = psastroModelBoresiteLimits;
+
+    // make an initial guess:
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_MAX | PS_STAT_MIN);
+
+    // center (Lo) = mean(Lo), RL = range / 2
+    psVectorStats (stats, Lo, NULL, NULL, 0);
+    params->data.F32[PAR_L0] = stats->sampleMean;
+    params->data.F32[PAR_RL] = (stats->max - stats->min) / 2.0;
+
+    // center (Mo) = mean(Mo), RM = range / 2
+    psVectorStats (stats, Mo, NULL, NULL, 0);
+    params->data.F32[PAR_M0] = stats->sampleMean;
+    params->data.F32[PAR_RM] = (stats->max - stats->min) / 2.0;
+
+    params->data.F32[PAR_P0] = 0.0;
+    params->data.F32[PAR_T0] = 0.0;
+
+    psMinimization *myMin = psMinimizationAlloc (15, 0.001);
+    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
+    
+    fprintf (stderr, "guess values:\n");
+    fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
+    fprintf (stderr, "Mo:  %f\n", params->data.F32[PAR_M0]);
+    fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
+    fprintf (stderr, "RM:  %f\n", params->data.F32[PAR_RM]);
+    fprintf (stderr, "P0:  %f\n", params->data.F32[PAR_P0]);
+    fprintf (stderr, "T0:  %f\n", params->data.F32[PAR_T0]);
+
+    // XXX skip the weights for now
+    psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psastroModelBoresite);
+
+    fprintf (stderr, "fitted values:\n");
+    fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
+    fprintf (stderr, "Mo:  %f\n", params->data.F32[PAR_M0]);
+    fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
+    fprintf (stderr, "RM:  %f\n", params->data.F32[PAR_RM]);
+    fprintf (stderr, "P0:  %f\n", params->data.F32[PAR_P0]);
+    fprintf (stderr, "T0:  %f\n", params->data.F32[PAR_T0]);
+
+    return true;
+}
+
+# if (0)
+# undef PAR_L0
+# undef PAR_RL
+# undef PAR_PHI
+# define PAR_L0  0
+# define PAR_RL  1
+# define PAR_PHI 2
+
+// we now have a set of observed L,M values.  fit these to the boresite model
+bool psastroModelFitBoresite_L (psVector *Lo, psVector *Mo, psVector *Po) {
+
+    assert (Lo->n > 2);
+    assert (Lo->n == Mo->n);
+    assert (Lo->n == Po->n);
+
+    // arrays to hold the data to be fitted
+    psArray *x = psArrayAlloc(Lo->n);
+    psVector *y = psVectorAlloc(Lo->n, PS_TYPE_F32);
+
+    for (int i = 0; i < Lo->n; i++) {
+
+	psVector *coord = NULL;
+
+	// L coordinate value
+	coord = psVectorAlloc (1, PS_TYPE_F32);
+	coord->data.F32[0] = Po->data.F32[i];
+	x->data[i] = coord;
+	y->data.F32[i] = Lo->data.F32[i];
+    }	
+
+    psVector *params = psVectorAlloc (3, PS_TYPE_F32);
+    
+    // create the minimization constraints
+    psMinConstraint *constraint = psMinConstraintAlloc();
+
+    // XXX for now, no parameter masks, skip checkLimits
+    // constraint->checkLimits = psastroModelBoresiteLimits;
+
+    // make an initial guess:
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_MAX | PS_STAT_MIN);
+
+    // center (Lo) = mean(Lo), RL = range / 2
+    psVectorStats (stats, Lo, NULL, NULL, 0);
+    params->data.F32[PAR_L0] = stats->sampleMean;
+    params->data.F32[PAR_RL] = (stats->max - stats->min) / 2.0;
+    params->data.F32[PAR_PHI] = 0.0;
+
+    psMinimization *myMin = psMinimizationAlloc (15, 0.001);
+    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
+    
+    fprintf (stderr, "guess values:\n");
+    fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
+    fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
+    fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]);
+
+    // XXX skip the weights for now
+    psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psastroModelBoresiteL);
+
+    fprintf (stderr, "fitted values:\n");
+    fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
+    fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
+    fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]);
+
+    return true;
+}
+
+// we now have a set of observed L,M values.  fit these to the boresite model
+bool psastroModelFitBoresite (psVector *Lo, psVector *Mo, psVector *Po) {
+
+    assert (Lo->n > 2);
+    assert (Lo->n == Mo->n);
+    assert (Lo->n == Po->n);
+
+    // arrays to hold the data to be fitted
+    psArray *x = psArrayAlloc(Lo->n);
+    psVector *y = psVectorAlloc(Lo->n, PS_TYPE_F32);
+
+    for (int i = 0; i < Lo->n; i++) {
+
+	psVector *coord = NULL;
+
+	// L coordinate value
+	coord = psVectorAlloc (1, PS_TYPE_F32);
+	coord->data.F32[0] = Po->data.F32[i];
+	x->data[i] = coord;
+	y->data.F32[i] = Mo->data.F32[i];
+    }	
+
+    psVector *params = psVectorAlloc (3, PS_TYPE_F32);
+    
+    // create the minimization constraints
+    psMinConstraint *constraint = psMinConstraintAlloc();
+
+    // XXX for now, no parameter masks, skip checkLimits
+    // constraint->checkLimits = psastroModelBoresiteLimits;
+
+    // make an initial guess:
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_MAX | PS_STAT_MIN);
+
+    // center (Lo) = mean(Lo), RL = range / 2
+    psVectorStats (stats, Mo, NULL, NULL, 0);
+    params->data.F32[PAR_L0] = stats->sampleMean;
+    params->data.F32[PAR_RL] = (stats->max - stats->min) / 2.0;
+    params->data.F32[PAR_PHI] = 0.0;
+
+    psMinimization *myMin = psMinimizationAlloc (15, 0.001);
+    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
+    
+    fprintf (stderr, "guess values:\n");
+    fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
+    fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
+    fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]);
+
+    // XXX skip the weights for now
+    psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psastroModelBoresiteM);
+
+    fprintf (stderr, "fitted values:\n");
+    fprintf (stderr, "Lo:  %f\n", params->data.F32[PAR_L0]);
+    fprintf (stderr, "RL:  %f\n", params->data.F32[PAR_RL]);
+    fprintf (stderr, "phi: %f\n", params->data.F32[PAR_PHI]);
+
+    return true;
+}
+# endif
Index: /trunk/psastro/src/psastroStandAlone.h
===================================================================
--- /trunk/psastro/src/psastroStandAlone.h	(revision 15886)
+++ /trunk/psastro/src/psastroStandAlone.h	(revision 15887)
@@ -22,3 +22,17 @@
 bool psastroModelAnalysis (pmConfig *config);
 
+bool psastroModelFitBoresite (psVector *Lo, psVector *Mo, psVector *Po);
+psF32 psastroModelBoresite (psVector *deriv, const psVector *params, const psVector *coord);
+
+psF32 psastroModelBoresiteL (psVector *deriv, const psVector *params, const psVector *coord);
+psF32 psastroModelBoresiteM (psVector *deriv, const psVector *params, const psVector *coord);
+
+// these are used to define the boresite model parameters
+# define PAR_L0  0  // Lo = params[0] 
+# define PAR_M0  1  // Mo = params[1] 
+# define PAR_RL  2  // RL = params[2] 
+# define PAR_RM  3  // RM = params[3] 
+# define PAR_P0  4  // P0 = params[4]
+# define PAR_T0  5  // phi = params[4]
+
 #endif