Index: trunk/Poubelle/archTOI.old/covsrt.c
===================================================================
--- trunk/Poubelle/archTOI.old/covsrt.c	(revision 566)
+++ trunk/Poubelle/archTOI.old/covsrt.c	(revision 577)
@@ -1,8 +1,8 @@
 #define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}
 
-void covsrt(float **covar, int ma, int ia[], int mfit)
+void covsrt(double **covar, int ma, int ia[], int mfit)
 {
 	int i,j,k;
-	float swap;
+	double swap;
 
 	for (i=mfit+1;i<=ma;i++)
Index: trunk/Poubelle/archTOI.old/gaussj.c
===================================================================
--- trunk/Poubelle/archTOI.old/gaussj.c	(revision 566)
+++ trunk/Poubelle/archTOI.old/gaussj.c	(revision 577)
@@ -4,9 +4,9 @@
 #define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
 
-void gaussj(float **a, int n, float **b, int m)
+void gaussj(double **a, int n, double **b, int m)
 {
 	int *indxc,*indxr,*ipiv;
 	int i,icol,irow,j,k,l,ll;
-	float big,dum,pivinv,temp;
+	double big,dum,pivinv,temp;
 
 	indxc=ivector(1,n);
Index: trunk/Poubelle/archTOI.old/gondolageom.cc
===================================================================
--- trunk/Poubelle/archTOI.old/gondolageom.cc	(revision 566)
+++ trunk/Poubelle/archTOI.old/gondolageom.cc	(revision 577)
@@ -67,5 +67,7 @@
 }
 
-void GondolaGeom::solveStars() {
+// $CHECK$ do a higher order fit ?
+int GondolaGeom::solveStars() {
+  if (nstars<2) return -1;
   staz /= nstars;
   st   /= nstars;
@@ -81,4 +83,6 @@
   if (azimut > 360) azimut -= 360;
   if (azimut < 0)   azimut += 360;
+
+  return 0;
 }
 
Index: trunk/Poubelle/archTOI.old/gondolageom.h
===================================================================
--- trunk/Poubelle/archTOI.old/gondolageom.h	(revision 566)
+++ trunk/Poubelle/archTOI.old/gondolageom.h	(revision 577)
@@ -15,5 +15,5 @@
   // Set orientation through SST
   void addStar(double deltasn, double az, double elv, double diod);
-  void solveStars();
+  int solveStars();
   
   void getPointing(double elv, double az, double& trueElv, double& trueAz);
Index: trunk/Poubelle/archTOI.old/lfit.c
===================================================================
--- trunk/Poubelle/archTOI.old/lfit.c	(revision 566)
+++ trunk/Poubelle/archTOI.old/lfit.c	(revision 577)
@@ -2,11 +2,11 @@
 #include "nrutil.h"
 
-void lfit(float x[], float y[], float sig[], int ndat, float a[], int ia[],
-	int ma, float **covar, float *chisq, void (*funcs)(float, float [], int))
+void lfit(double x[], double y[], double sig[], int ndat, double a[], int ia[],
+	int ma, double **covar, double *chisq, void (*funcs)(double, double [], int))
 {
-	void covsrt(float **covar, int ma, int ia[], int mfit);
-	void gaussj(float **a, int n, float **b, int m);
+	void covsrt(double **covar, int ma, int ia[], int mfit);
+	void gaussj(double **a, int n, double **b, int m);
 	int i,j,k,l,m,mfit=0;
-	float ym,wt,sum,sig2i,**beta,*afunc;
+	double ym,wt,sum,sig2i,**beta,*afunc;
 
 	beta=matrix(1,ma,1,1);
Index: trunk/Poubelle/archTOI.old/nrutil.c
===================================================================
--- trunk/Poubelle/archTOI.old/nrutil.c	(revision 566)
+++ trunk/Poubelle/archTOI.old/nrutil.c	(revision 577)
@@ -18,10 +18,10 @@
 }*/
 
-float *vector(long nl, long nh)
-/* allocate a float vector with subscript range v[nl..nh] */
-{
-	float *v;
-
-	v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
+double *vector(long nl, long nh)
+/* allocate a double vector with subscript range v[nl..nh] */
+{
+	double *v;
+
+	v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
 	if (!v) nrerror("allocation failure in vector()");
 	return v-nl+NR_END;
@@ -68,29 +68,5 @@
 }
 
-float **matrix(long nrl, long nrh, long ncl, long nch)
-/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
-{
-	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
-	float **m;
-
-	/* allocate pointers to rows */
-	m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
-	if (!m) nrerror("allocation failure 1 in matrix()");
-	m += NR_END;
-	m -= nrl;
-
-	/* allocate rows and set pointers to them */
-	m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
-	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
-	m[nrl] += NR_END;
-	m[nrl] -= ncl;
-
-	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
-
-	/* return pointer to array of pointers to rows */
-	return m;
-}
-
-double **dmatrix(long nrl, long nrh, long ncl, long nch)
+double **matrix(long nrl, long nrh, long ncl, long nch)
 /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
 {
@@ -116,4 +92,28 @@
 }
 
+double **dmatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	double **m;
+
+	/* allocate pointers to rows */
+	m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
 int **imatrix(long nrl, long nrh, long ncl, long nch)
 /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
@@ -141,13 +141,13 @@
 }
 
-float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
+double **submatrix(double **a, long oldrl, long oldrh, long oldcl, long oldch,
 	long newrl, long newcl)
 /* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
 {
 	long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
-	float **m;
+	double **m;
 
 	/* allocate array of pointers to rows */
-	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
+	m=(double **) malloc((size_t) ((nrow+NR_END)*sizeof(double*)));
 	if (!m) nrerror("allocation failure in submatrix()");
 	m += NR_END;
@@ -161,6 +161,6 @@
 }
 
-float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
-/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
+double **convert_matrix(double *a, long nrl, long nrh, long ncl, long nch)
+/* allocate a double matrix m[nrl..nrh][ncl..nch] that points to the matrix
 declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
 and ncol=nch-ncl+1. The routine should be called with the address
@@ -168,8 +168,8 @@
 {
 	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
-	float **m;
-
-	/* allocate pointers to rows */
-	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
+	double **m;
+
+	/* allocate pointers to rows */
+	m=(double **) malloc((size_t) ((nrow+NR_END)*sizeof(double*)));
 	if (!m) nrerror("allocation failure in convert_matrix()");
 	m += NR_END;
@@ -183,12 +183,12 @@
 }
 
-float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
-/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+double ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
+/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
 {
 	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
-	float ***t;
+	double ***t;
 
 	/* allocate pointers to pointers to rows */
-	t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
+	t=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double**)));
 	if (!t) nrerror("allocation failure 1 in f3tensor()");
 	t += NR_END;
@@ -196,5 +196,5 @@
 
 	/* allocate pointers to rows and set pointers to them */
-	t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
+	t[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double*)));
 	if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
 	t[nrl] += NR_END;
@@ -202,5 +202,5 @@
 
 	/* allocate rows and set pointers to them */
-	t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
+	t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(double)));
 	if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
 	t[nrl][ncl] += NR_END;
@@ -218,6 +218,6 @@
 }
 
-void free_vector(float *v, long nl, long nh)
-/* free a float vector allocated with vector() */
+void free_vector(double *v, long nl, long nh)
+/* free a double vector allocated with vector() */
 {
 	free((FREE_ARG) (v+nl-NR_END));
@@ -248,6 +248,6 @@
 }
 
-void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
-/* free a float matrix allocated by matrix() */
+void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
+/* free a double matrix allocated by matrix() */
 {
 	free((FREE_ARG) (m[nrl]+ncl-NR_END));
@@ -269,5 +269,5 @@
 }
 
-void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
+void free_submatrix(double **b, long nrl, long nrh, long ncl, long nch)
 /* free a submatrix allocated by submatrix() */
 {
@@ -275,5 +275,5 @@
 }
 
-void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
+void free_convert_matrix(double **b, long nrl, long nrh, long ncl, long nch)
 /* free a matrix allocated by convert_matrix() */
 {
@@ -281,7 +281,7 @@
 }
 
-void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
+void free_f3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
 	long ndl, long ndh)
-/* free a float f3tensor allocated by f3tensor() */
+/* free a double f3tensor allocated by f3tensor() */
 {
 	free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
@@ -309,11 +309,11 @@
 }
 
-float *vector(nl,nh)
-long nh,nl;
-/* allocate a float vector with subscript range v[nl..nh] */
-{
-	float *v;
-
-	v=(float *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(float)));
+double *vector(nl,nh)
+long nh,nl;
+/* allocate a double vector with subscript range v[nl..nh] */
+{
+	double *v;
+
+	v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double)));
 	if (!v) nrerror("allocation failure in vector()");
 	return v-nl+NR_END;
@@ -364,30 +364,5 @@
 }
 
-float **matrix(nrl,nrh,ncl,nch)
-long nch,ncl,nrh,nrl;
-/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
-{
-	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
-	float **m;
-
-	/* allocate pointers to rows */
-	m=(float **) malloc((unsigned int)((nrow+NR_END)*sizeof(float*)));
-	if (!m) nrerror("allocation failure 1 in matrix()");
-	m += NR_END;
-	m -= nrl;
-
-	/* allocate rows and set pointers to them */
-	m[nrl]=(float *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float)));
-	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
-	m[nrl] += NR_END;
-	m[nrl] -= ncl;
-
-	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
-
-	/* return pointer to array of pointers to rows */
-	return m;
-}
-
-double **dmatrix(nrl,nrh,ncl,nch)
+double **matrix(nrl,nrh,ncl,nch)
 long nch,ncl,nrh,nrl;
 /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
@@ -414,4 +389,29 @@
 }
 
+double **dmatrix(nrl,nrh,ncl,nch)
+long nch,ncl,nrh,nrl;
+/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	double **m;
+
+	/* allocate pointers to rows */
+	m=(double **) malloc((unsigned int)((nrow+NR_END)*sizeof(double*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(double *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
 int **imatrix(nrl,nrh,ncl,nch)
 long nch,ncl,nrh,nrl;
@@ -440,14 +440,14 @@
 }
 
-float **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
-float **a;
+double **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
+double **a;
 long newcl,newrl,oldch,oldcl,oldrh,oldrl;
 /* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
 {
 	long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
-	float **m;
+	double **m;
 
 	/* allocate array of pointers to rows */
-	m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));
+	m=(double **) malloc((unsigned int) ((nrow+NR_END)*sizeof(double*)));
 	if (!m) nrerror("allocation failure in submatrix()");
 	m += NR_END;
@@ -461,8 +461,8 @@
 }
 
-float **convert_matrix(a,nrl,nrh,ncl,nch)
-float *a;
-long nch,ncl,nrh,nrl;
-/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
+double **convert_matrix(a,nrl,nrh,ncl,nch)
+double *a;
+long nch,ncl,nrh,nrl;
+/* allocate a double matrix m[nrl..nrh][ncl..nch] that points to the matrix
 declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
 and ncol=nch-ncl+1. The routine should be called with the address
@@ -470,8 +470,8 @@
 {
 	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
-	float **m;
-
-	/* allocate pointers to rows */
-	m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));
+	double **m;
+
+	/* allocate pointers to rows */
+	m=(double **) malloc((unsigned int) ((nrow+NR_END)*sizeof(double*)));
 	if (!m)	nrerror("allocation failure in convert_matrix()");
 	m += NR_END;
@@ -485,13 +485,13 @@
 }
 
-float ***f3tensor(nrl,nrh,ncl,nch,ndl,ndh)
+double ***f3tensor(nrl,nrh,ncl,nch,ndl,ndh)
 long nch,ncl,ndh,ndl,nrh,nrl;
-/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
 {
 	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
-	float ***t;
+	double ***t;
 
 	/* allocate pointers to pointers to rows */
-	t=(float ***) malloc((unsigned int)((nrow+NR_END)*sizeof(float**)));
+	t=(double ***) malloc((unsigned int)((nrow+NR_END)*sizeof(double**)));
 	if (!t) nrerror("allocation failure 1 in f3tensor()");
 	t += NR_END;
@@ -499,5 +499,5 @@
 
 	/* allocate pointers to rows and set pointers to them */
-	t[nrl]=(float **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float*)));
+	t[nrl]=(double **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double*)));
 	if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
 	t[nrl] += NR_END;
@@ -505,5 +505,5 @@
 
 	/* allocate rows and set pointers to them */
-	t[nrl][ncl]=(float *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(float)));
+	t[nrl][ncl]=(double *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(double)));
 	if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
 	t[nrl][ncl] += NR_END;
@@ -522,7 +522,7 @@
 
 void free_vector(v,nl,nh)
-float *v;
-long nh,nl;
-/* free a float vector allocated with vector() */
+double *v;
+long nh,nl;
+/* free a double vector allocated with vector() */
 {
 	free((FREE_ARG) (v+nl-NR_END));
@@ -562,7 +562,7 @@
 
 void free_matrix(m,nrl,nrh,ncl,nch)
-float **m;
-long nch,ncl,nrh,nrl;
-/* free a float matrix allocated by matrix() */
+double **m;
+long nch,ncl,nrh,nrl;
+/* free a double matrix allocated by matrix() */
 {
 	free((FREE_ARG) (m[nrl]+ncl-NR_END));
@@ -589,5 +589,5 @@
 
 void free_submatrix(b,nrl,nrh,ncl,nch)
-float **b;
+double **b;
 long nch,ncl,nrh,nrl;
 /* free a submatrix allocated by submatrix() */
@@ -597,5 +597,5 @@
 
 void free_convert_matrix(b,nrl,nrh,ncl,nch)
-float **b;
+double **b;
 long nch,ncl,nrh,nrl;
 /* free a matrix allocated by convert_matrix() */
@@ -605,7 +605,7 @@
 
 void free_f3tensor(t,nrl,nrh,ncl,nch,ndl,ndh)
-float ***t;
+double ***t;
 long nch,ncl,ndh,ndl,nrh,nrl;
-/* free a float f3tensor allocated by f3tensor() */
+/* free a double f3tensor allocated by f3tensor() */
 {
 	free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
Index: trunk/Poubelle/archTOI.old/nrutil.h
===================================================================
--- trunk/Poubelle/archTOI.old/nrutil.h	(revision 566)
+++ trunk/Poubelle/archTOI.old/nrutil.h	(revision 577)
@@ -2,5 +2,5 @@
 #define _NR_UTILS_H_
 
-static float sqrarg;
+static double sqrarg;
 #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
 
@@ -16,9 +16,9 @@
         (dminarg1) : (dminarg2))
 
-static float maxarg1,maxarg2;
+static double maxarg1,maxarg2;
 #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
         (maxarg1) : (maxarg2))
 
-static float minarg1,minarg2;
+static double minarg1,minarg2;
 #define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
         (minarg1) : (minarg2))
@@ -45,27 +45,27 @@
 
 void nrerror(char error_text[]);
-float *vector(long nl, long nh);
+double *vector(long nl, long nh);
 int *ivector(long nl, long nh);
 unsigned char *cvector(long nl, long nh);
 unsigned long *lvector(long nl, long nh);
 double *dvector(long nl, long nh);
-float **matrix(long nrl, long nrh, long ncl, long nch);
+double **matrix(long nrl, long nrh, long ncl, long nch);
 double **dmatrix(long nrl, long nrh, long ncl, long nch);
 int **imatrix(long nrl, long nrh, long ncl, long nch);
-float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
+double **submatrix(double **a, long oldrl, long oldrh, long oldcl, long oldch,
 	long newrl, long newcl);
-float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);
-float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
-void free_vector(float *v, long nl, long nh);
+double **convert_matrix(double *a, long nrl, long nrh, long ncl, long nch);
+double ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
+void free_vector(double *v, long nl, long nh);
 void free_ivector(int *v, long nl, long nh);
 void free_cvector(unsigned char *v, long nl, long nh);
 void free_lvector(unsigned long *v, long nl, long nh);
 void free_dvector(double *v, long nl, long nh);
-void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
+void free_matrix(double **m, long nrl, long nrh, long ncl, long nch);
 void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
 void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
-void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);
-void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);
-void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
+void free_submatrix(double **b, long nrl, long nrh, long ncl, long nch);
+void free_convert_matrix(double **b, long nrl, long nrh, long ncl, long nch);
+void free_f3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
 	long ndl, long ndh);
 
@@ -74,9 +74,9 @@
 
 void nrerror();
-float *vector();
-float **matrix();
-float **submatrix();
-float **convert_matrix();
-float ***f3tensor();
+double *vector();
+double **matrix();
+double **submatrix();
+double **convert_matrix();
+double ***f3tensor();
 double *dvector();
 double **dmatrix();
Index: trunk/Poubelle/archTOI.old/starmatcher.cc
===================================================================
--- trunk/Poubelle/archTOI.old/starmatcher.cc	(revision 566)
+++ trunk/Poubelle/archTOI.old/starmatcher.cc	(revision 577)
@@ -8,4 +8,7 @@
 #include "archparam.h"
 #include "gondolageom.h"
+#include "polfitclip.h"
+
+#define STARDUMP
 
 extern "C" {
@@ -14,12 +17,12 @@
 #include "nrutil.h"
 
-void lfit(float x[], float y[], float sig[], int ndat, float a[], int ia[],
-	int ma, float **covar, float *chisq, void (*funcs)(float, float [], int));
-
-void polfunc(float x, float afunc[], int ma);
-void sinfunc(float x, float afunc[], int ma);
-}
-
-void polfunc(float x, float afunc[], int ma) {
+void lfit(double x[], double y[], double sig[], int ndat, double a[], int ia[],
+	int ma, double **covar, double *chisq, void (*funcs)(double, double [], int));
+
+void polfunc(double x, double afunc[], int ma);
+void sinfunc(double x, double afunc[], int ma);
+}
+
+void polfunc(double x, double afunc[], int ma) {
   afunc[1] = 1;
   for (int i=2; i<=ma; i++)
@@ -27,5 +30,5 @@
 }
 
-void sinfunc(float x, float afunc[], int /*ma*/) {
+void sinfunc(double x, double afunc[], int /*ma*/) {
   afunc[1] = cos(x);
   afunc[2] = sin(x);
@@ -34,8 +37,8 @@
 
 
-float polval(float x, float a[], int ma);
-
-float polval(float x, float a[], int ma) {
-  float r = a[ma];
+double polval(double x, double a[], int ma);
+
+double polval(double x, double a[], int ma) {
+  double r = a[ma];
   for (int i=ma-1; i>0; i--) {
     r = r*x+a[i];
@@ -186,4 +189,5 @@
 static ofstream pendstream("pendul.dat");
 #endif
+static ofstream logstream("starmatch.log");
 
 void StarMatcher::dataFeed(SSTEtoile const& x) {
@@ -295,5 +299,5 @@
   
   // new set of matched stars... Clean, and get parameters...
-  // We don't want more than 20 seconds of data
+  // We don't want more than 30 seconds of data
   
   if (matchStars.empty()) return;
@@ -303,6 +307,10 @@
   deque<matchStar>::iterator it;
   for (it = matchStars.begin();  it!=matchStars.end(); it++) {
-     if ((snEnd - (*it).SN)*archParam.acq.perEch < 20) 
+     if ((snEnd - (*it).SN)*archParam.acq.perEch < 30 ||
+         (*it).seq > lastCleanSave) 
        break;
+  }
+  if (it != matchStars.begin()) {
+    it--;
   }
   if (it != matchStars.begin()) {
@@ -314,5 +322,5 @@
   int nskip=0;
   for (it = matchStars.begin();  it!=matchStars.end(); it++) {
-     if ((snEnd - (*it).SN)*archParam.acq.perEch < 7) 
+     if ((snEnd - (*it).SN)*archParam.acq.perEch < 5) 
        break;
      nskip++;
@@ -334,4 +342,6 @@
     if (burstLen >= 4) {
       for (deque<matchStar>::iterator it2=lastIt; it2 != it1; it2++) {
+        //if ((*it2).ok) 
+        //  logstream << "  kill " << (*it2).seq << " " << setprecision(11) << (*it2).SN << " burst" << '\n';
         (*it2).ok=false;
       }
@@ -343,70 +353,80 @@
   // we fit the data to a polynomial, with clipping...
   
-  float* sn   =  ::vector(1, matchStars.size());
-  float* elv0 =  ::vector(1, matchStars.size());
-  float* azi  =  ::vector(1, matchStars.size());
-  float* sig  =  ::vector(1, matchStars.size());
-  float* ae   =  ::vector(1, 3);
-  float* aa   =  ::vector(1, 3);
+  //double* sn   =  ::vector(1, matchStars.size());
+  double* elv0 =  ::vector(1, matchStars.size());
+  double* azi  =  ::vector(1, matchStars.size());
+  double* sig  =  ::vector(1, matchStars.size());
+  //double* ae   =  ::vector(1, 3);
+  double* aa   =  ::vector(1, 3);
   int*  ia    = ivector(1, 3);
-  float** cov =  matrix(1, 3, 1, 3);
+  double** cov =  matrix(1, 3, 1, 3);
   int ndata;
   
-  long sn0 = matchStars.front().SN;
-  long snmin;
-  long snmax;
-  for (int i=0; i<4; i++) {
-    ndata = 0;
-    snmin = 99999999;
-    snmax = -99999999;
-    for (deque<matchStar>::iterator it1 = it ;  it1!=matchStars.end(); it1++) {
-      matchStar s = (*it1);
-      if (!s.ok) continue;
-      double delv, daz;
-      if (i) {
-        delv = polval(s.SN-sn0, ae, 3)-(s.elvGSC - s.nDiode*1.41/45.);
-        daz  = polval(s.SN-sn0, aa, 3)- s.azGSC;
-        if (daz>=180) daz -= 360;
-        if (daz<-180) daz += 360;
+  //long sn0 = matchStars.front().SN;
+  long sn0 = (*it).SN;
+  PolFitClip2 fitElvAz(matchStars.size(), 2); fitElvAz.setClip(0.1,0,2,3);
+  ndata = 0;
+
+  double oldAz = -1;
+  for (deque<matchStar>::iterator it1 = it ;  it1!=matchStars.end(); it1++) {
+    matchStar s1 = (*it1);
+    if (!s1.ok) continue;
+    double az = s1.azGSC;
+    if (ndata > 0 && az - oldAz >  180) az -= 360;
+    if (ndata > 0 && az - oldAz < -180) az += 360;
+    fitElvAz.addData(s1.SN-sn0, s1.elvGSC - s1.nDiode*1.41/45., az);
+    oldAz = az;
+    ndata++;
+  }
+   
+  double celv[3], caz[3];
+  if (fitElvAz.doFit(celv,caz)) return;
+  //if (fitElvAz.doFit()) return;
+  
+  //logstream << "*** Fit sig=" << fitElvAz.getSigmaY() << " " << fitElvAz.getSigmaZ()
+  //          << " n  =" << fitElvAz.getNData() << " " << fitElvAz.getNDataUsed() 
+  //          << " SN :" << fitElvAz.getXMin() + sn0 << " - " << fitElvAz.getXMax() + sn0 << '\n';
+  //logstream << " sn0 = " << sn0 << "; snmin =" << fitElvAz.getXMin() + sn0 << "; snmax =" 
+  //          << fitElvAz.getXMax() + sn0  << '\n';
+  //logstream << " fitelv[x_] := " << celv[2] << " x^2 + " << celv[1] << " x + " << celv[0] << '\n';
+  //logstream << " fitaz[x_]  := " <<  caz[2] << " x^2 + " <<  caz[1] << " x + " <<  caz[0] << '\n';
+  
+  //if (fitElvAz.getSigmaY() > 0.05 || fitElvAz.getSigmaZ() > 0.05) return;
+  if (fitElvAz.getNDataUsed() < 5 ||
+      (double)fitElvAz.getNDataUsed()/fitElvAz.getNData() < .5) return;
+      
+  double dcutElv = fitElvAz.getSigmaY()*3;
+  double dcutAz  = fitElvAz.getSigmaZ()*3;
+  
+  if (dcutElv < 0.05) dcutElv = 0.05;
+  if (dcutAz  < 0.05) dcutAz  = 0.05;
+
+  // don't kill borders of fit....
+  //if (matchStars.end() - it > 6)
+  //  for (deque<matchStar>::iterator it1 = it+3 ;  it1!=matchStars.end()-3; it1++) {
+  for (deque<matchStar>::iterator it1 = it ;  it1!=matchStars.end(); it1++) {
+      matchStar sss = (*it1);
+      if (!sss.ok) continue;
+      if (fabs(fitElvAz.valueY(sss.SN-sn0)-(sss.elvGSC - sss.nDiode*1.41/45.)) > dcutElv) {
+        (*it1).ok = false;
+        //logstream << "  kill " << sss.seq << " " << setprecision(11) << sss.SN << " "
+        //          << fitElvAz.valueY(sss.SN-sn0)-(sss.elvGSC - sss.nDiode*1.41/45.) << '\n';
+        continue;
       }
-      double dcutelv=0.2;
-      double dcutaz =0.4;
-      if (i>=2) {
-         dcutelv = 0.05;
-         dcutaz  = 0.1;
+      double daz = fitElvAz.valueZ(sss.SN-sn0) - sss.azGSC;
+      if (daz>=180) daz -= 360;
+      if (daz<-180) daz += 360;
+      if (fabs(daz) > dcutAz) (*it1).ok = false;
+      if (!(*it1).ok) {
+        //logstream << "  kill " << sss.seq << " " << setprecision(11) << sss.SN << " "
+        //          << fitElvAz.valueY(sss.SN-sn0)-(sss.elvGSC - sss.nDiode*1.41/45.) << " "
+        //          << daz << '\n';
       }
-      if (i>=3) {
-         dcutelv = 0.02;
-         dcutaz  = 0.03;
-      }
-      if (i == 0 || ((fabs(delv)<dcutelv && fabs(daz)<dcutaz) && i<3)) {
-          ndata++;
-          sn  [ndata] = s.SN-sn0;
-          elv0[ndata] = s.elvGSC - s.nDiode*1.41/45.;
-          azi [ndata] = s.azGSC; // $CHECK$  ou ajuster difference avec az galcross ?
-          if (ndata>1 && azi[ndata] - azi[ndata-1] > 180) azi[ndata] -= 360;
-          if (ndata>1 && azi[ndata] - azi[ndata-1] < -180) azi[ndata] += 360;
-          sig [ndata] = 0.01;
-          if (s.SN-sn0 > snmax) snmax = s.SN-sn0;
-          if (s.SN-sn0 < snmin) snmin = s.SN-sn0;
-      }
-      if ((fabs(delv)>=dcutelv || fabs(daz)>=dcutaz) && i==3) {
-          (*it1).ok = false; 
-      }
-    }
-    if (i==3) break;
-    if (ndata<5) return; // on ne peut rien faire
-    ia[1]   =  ia[2] = ia[3] = 1;
-    float chi2;
-    try{
-      lfit(sn, elv0, sig, ndata, ae, ia, 3, cov, &chi2, polfunc);
-      lfit(sn, azi,  sig, ndata, aa, ia, 3, cov, &chi2, polfunc);
-    } catch(string s) {
-      return;
-    }
-  }
-  
-  for (deque<matchStar>::iterator it1 = it ;  it1!=matchStars.end(); it1++) {
+  }
+    
+  bool gotNewStars = false;
+  for (deque<matchStar>::iterator it1 = matchStars.begin() ;  it1!=it; it1++) {
     if ((*it1).ok && (*it1).seq > lastCleanSave) {
+      gotNewStars = true;
       lastCleanSave = (*it1).seq;
 #ifdef STARDUMP
@@ -425,4 +445,6 @@
   }
   
+  if (!gotNewStars) return;
+  
   // On a des etoiles nettoyees, on va trouver amplitude et phase du
   // signal en elevation, ce qui va nous donner les deux angles d'Euler
@@ -430,56 +452,84 @@
   
   // Il faut avoir une periode entiere ou pas loin, sinon on ne peut
-  // rien dire simplement....
-  
-  it = matchStars.end(); it--;
-  if ((((*it).SN) - (*matchStars.begin()).SN)*archParam.acq.perEch < 17) return;
-  
-  long snmid = (((*it).SN) + (*matchStars.begin()).SN)/2;
+  // rien dire simplement.... -> we want to run on the last 18 seconds of
+  // data before the last fully cleaned star (it).
+  
+  deque<matchStar>::iterator itstart;
+  
+  for (itstart = matchStars.begin();  itstart != it; itstart++) {
+     if (((*it).SN - (*itstart).SN)*archParam.acq.perEch < 19) 
+       break;
+  }
+  
+  if (((*it).SN - (*itstart).SN)*archParam.acq.perEch < 15) return;
+
+  
+ // it = matchStars.end(); it--;
+ // if (((*it).SN - matchStars.front().SN)*archParam.acq.perEch < 17) return;
+  
+  // $CHECK$ utiliser plutot le SN moyen/median de tous les points effectivement utilises.
+  long snmid = ((*it).SN + (*itstart).SN)/2;
     
   ndata=0;  
-  
-  for (deque<matchStar>::iterator it1 = matchStars.begin() ;  it1!=matchStars.end(); it1++) {
-    if (!(*it1).ok) continue;
+  double snmean = 0;
+
+  logstream << "PendFit : " <<  setprecision(11) << (*itstart).SN << '-' << (*it).SN <<  " "
+            << setprecision(4)
+            << ((*it).SN - (*itstart).SN)*archParam.acq.perEch << " " ;
+              
+  for (deque<matchStar>::iterator it1 = itstart ;  it1!=it; it1++) {
+    matchStar st = *it1;
+    if (!st.ok) continue;
     ndata++;
-    azi[ndata]  = (*it1).azGSC * 3.1415926/180;
-    elv0[ndata] = (*it1).elvGSC - (*it1).nDiode*1.41/45.;
+    snmean += st.SN;
+    azi[ndata]  = st.azGSC * 3.1415926/180;
+    elv0[ndata] = st.elvGSC - st.nDiode*1.41/45.;
     sig[ndata]  = 0.01;
   }
-   ia[1]   =  ia[2] = 1;
-   ia[3] = 0;
-   aa[3] = GondolaGeom::elevSST0;// do not fit elv0
+  if (ndata) snmean /= ndata;
+  
+  ia[1]   =  ia[2] = 1;
+  ia[3] = 0;
+  aa[3] = GondolaGeom::elevSST0;// do not fit elv0
 
   if (ndata<5) return;
-  float chi2;
+  double chi2;
   try {
     lfit(azi, elv0,  sig, ndata, aa, ia, 3, cov, &chi2, sinfunc);
-  } catch(string s) {
+  } catch(string st) {
     return;
   }
 
-  double c = aa[1];
-  double s = aa[2];
+  double cc = aa[1];
+  double ss = aa[2];
+  
+  logstream << setprecision(11) << snmean << setprecision(4)
+            << " cs=" << cc << " " << ss << " chi2r=" << chi2/ndata 
+            << " cov " << cov[1][1] << " " << cov[2][2] << '\n';
   
   // Get rid of bad fits. The cuts are rather ad hoc
   
   //if (aa[3] < 39.64 || aa[3] > 39.68) return;
-  if (chi2/ndata > 4) return;
+  if (chi2/ndata > 9) return;
   if (cov[1][1] > 0.0001) return;
   if (cov[2][2] > 0.0001) return;
 
-  double ampl = sqrt(c*c+s*s);
-  double phase = atan2(c,s)/(3.1415926/180);
-  
-#ifdef STARDUMP
-  pendstream << snmid << " " << ampl << " " << phase << " " << ndata << " " << chi2/ndata <<  
-       " "  << cov[1][1] << " " << cov[2][2] << '\n';
-#endif
+  double ampl = sqrt(cc*cc+ss*ss);
+  double phase = atan2(cc,ss)/(3.1415926/180);
+  
   
   pendulInfo info;
-  info.SN          =   snmid;
-  info.azPendul    =   phase;
+  info.SN          =   snmean;
+  info.azPendul    =   180-phase;
+  if (info.azPendul > 360) info.azPendul -= 360;
+  if (info.azPendul <   0) info.azPendul += 360;
   info.angPendul   =   ampl;
   pendulInfos[info.SN] = info;
 
+#ifdef STARDUMP
+  pendstream << setprecision(11) << snmean << " " 
+             << setprecision(4) << ampl << " " << info.azPendul << " " << ndata << " " << chi2/ndata <<  " "
+             << cov[1][1] << " " << cov[2][2] << '\n';
+#endif
   /*
   double snum = (matchStars.front().SN + matchStars.back().SN)/2-sn0;
@@ -497,13 +547,135 @@
 //  }
   
-  free_vector(sn, 1, matchStars.size());
+  //free_vector(sn, 1, matchStars.size());
   free_vector(elv0, 1, matchStars.size());
   free_vector(azi, 1, matchStars.size());
   free_vector(sig, 1, matchStars.size());
-  free_vector(ae, 1, 3);
+  //free_vector(ae, 1, 3);
   free_vector(aa, 1, 3);
   free_ivector(ia, 1, matchStars.size());
   free_matrix(cov, 1, 3, 1, 3);
 }
+
+
+// $CHECK$ do a polynomial fit with several points...
+int StarMatcher::getPendulInfo(double sampleNum, pendulInfo& info) {
+
+  PolFitClip2 fitPendul(30,2);
+
+  map<double, pendulInfo>::iterator i = pendulInfos.lower_bound(sampleNum);
+  if (i == pendulInfos.begin() && (*i).second.SN >= sampleNum) return -1;
+  if (i == pendulInfos.end()   && (*i).second.SN <= sampleNum) return -1;
+  
+  if ((*i).second.SN > sampleNum) i--;
+
+  int nn=0;
+  double aziprev=0, azicur=0, azi0=0;
+  for (map<double, pendulInfo>::iterator ii=i; ii != pendulInfos.begin(); ii--) {
+    nn++;
+    pendulInfo inf1 = (*ii).second;
+    aziprev = azicur;
+    azicur = inf1.azPendul;
+    if (nn==1) azi0 = azicur;
+    if (nn>1 && azicur - aziprev > 180)  azicur -= 360;
+    if (nn>1 && azicur - aziprev < -180) azicur += 360;
+    fitPendul.addData(inf1.SN, inf1.angPendul, azicur);
+    if (nn>=5) break;
+  }
+  
+  azicur = azi0;
+  if (i != pendulInfos.end()) i++;
+  for (map<double, pendulInfo>::iterator ii=i; ii != pendulInfos.end(); ii++) {
+    nn++;
+    pendulInfo inf1 = (*ii).second;
+    aziprev = azicur;
+    azicur = inf1.azPendul;
+    if (nn>1 && azicur - aziprev > 180)  azicur -= 360;
+    if (nn>1 && azicur - aziprev < -180) azicur += 360;
+    fitPendul.addData(inf1.SN, inf1.angPendul, azicur);
+    if (nn>=10) break;
+  }
+  
+  if (fitPendul.doFit()) return -1;
+  
+  info.SN = sampleNum;
+  info.azPendul   = fitPendul.valueZ(sampleNum);
+  if (info.azPendul > 360) info.azPendul -= 360;
+  if (info.azPendul <   0) info.azPendul += 360;
+  info.angPendul  = fitPendul.valueY(sampleNum);
+  return 0;
+}
+
+#if 0
+
+int StarMatcher::getPendulInfo(double sampleNum, pendulInfo& info) {
+  static double* sn    = ::vector(1, 100);
+  static double* aziP  = ::vector(1, 100);
+  static double* ampP  = ::vector(1, 100);
+  static double* sig   = ::vector(1, 100);
+  static double* aAzi  = ::vector(1, 3);
+  static double* aAmp  = ::vector(1, 3);
+  static int*   ia    = ::ivector(1,3);
+  static double** cov  = ::matrix(1,3,1,3);
+  int ndata = 0;
+  map<double, pendulInfo>::iterator i = pendulInfos.lower_bound(sampleNum);
+  if (i == pendulInfos.begin() && (*i).second.SN >= sampleNum) return -1;
+  if (i == pendulInfos.end()   && (*i).second.SN <= sampleNum) return -1;
+  
+  if ((*i).second.SN > sampleNum) i--;
+
+  int nn=0;
+  for (map<double, pendulInfo>::iterator ii=i; ii != pendulInfos.begin(); ii--) {
+    nn++;
+    ndata++;
+    pendulInfo inf1 = (*ii).second;
+    sn[ndata]   = inf1.SN;
+    ampP[ndata] = inf1.angPendul;
+    aziP[ndata] = inf1.azPendul;
+    int prev = ndata-1;
+    if (ndata>1 && aziP[ndata] - aziP[prev] > 180)  aziP[ndata] -= 360;
+    if (ndata>1 && aziP[ndata] - aziP[prev] < -180) aziP[ndata] += 360;
+    sig[ndata]  = 1;
+    if (nn>=50) break;
+  }
+  
+  nn=0;
+  if (i != pendulInfos.end()) i++;
+  for (map<double, pendulInfo>::iterator ii=i; ii != pendulInfos.end(); ii++) {
+    nn++;
+    ndata++;
+    pendulInfo inf1 = (*ii).second;
+    sn[ndata]   = inf1.SN;
+    ampP[ndata] = inf1.angPendul;
+    aziP[ndata] = inf1.azPendul;
+    int prev = ndata-1;
+    if (nn==1) prev=1;
+    if (ndata>1 && aziP[ndata] - aziP[prev] > 180)  aziP[ndata] -= 360;
+    if (ndata>1 && aziP[ndata] - aziP[prev] < -180) aziP[ndata] += 360;
+    sig[ndata]  = 1;
+    if (nn>=50) break;
+  }
+  
+  if (ndata < 3) return -1;
+  
+  ia[1] = ia[2] = ia[3] = 1;
+  double chi2;
+  try {
+    lfit(sn, aziP,  sig, ndata, aAzi, ia, 3, cov, &chi2, polfunc);
+    lfit(sn, ampP,  sig, ndata, aAmp, ia, 3, cov, &chi2, polfunc);
+  } catch(string st) {
+    return -1;
+  }
+
+  info.SN = sampleNum;
+  info.azPendul   = polval(sampleNum, aAzi, 3);
+  if (info.azPendul > 360) info.azPendul -= 360;
+  if (info.azPendul <   0) info.azPendul += 360;
+  info.angPendul  = polval(sampleNum, aAmp, 3);
+  return 0;
+}
+
+#endif
+
+#if 0
 
 int StarMatcher::getPendulInfo(double sampleNum, pendulInfo& info) {
@@ -527,4 +699,6 @@
 }
 
+#endif
+
 
 double StarMatcher::getValue(long sampleNum, TOI const& toi) {
@@ -546,24 +720,33 @@
   if ((*i).second.SN > sampleNum) i--;
   
+  // $CHECK$ if time spent here, can keep a GondolaGeom object for several
+  // samples...
   GondolaGeom geom;
   geom.setEarthPos((*i).second.lon,(*i).second.lat);
   geom.setTSid((*i).second.ts);
-  
+  geom.setPendulation(pendul.azPendul, pendul.angPendul);
+  
+  int ns=0;
   for (map<double, posInfo>::iterator it=i; it != posInfos.end(); it++) {
     posInfo s = (*it).second;
     double delsn = s.SN - sampleNum;
-    if (delsn * archParam.acq.perEch > 1) break;
+    ns++;
+    //if (delsn * archParam.acq.perEch > 1 && ns > 4) break;
+    if (delsn * archParam.acq.perEch > 5) break;
     geom.addStar(delsn, s.azStar, s.elvStar, s.diodStar);
   }
   
   if (i != posInfos.begin()) i--;
+  ns = 0;
   for (map<double, posInfo>::iterator it=i; it != posInfos.begin(); it--) {
     posInfo s = (*it).second;
     double delsn = s.SN - sampleNum;
-    if (-delsn * archParam.acq.perEch > 1) break;
+    ns++;
+    //if (-delsn * archParam.acq.perEch > 1 && ns > 4) break;
+    if (-delsn * archParam.acq.perEch > 5) break;
     geom.addStar(delsn, s.azStar, s.elvStar, s.diodStar);
   }
 
-  geom.solveStars();
+  if (geom.solveStars()) return -99999;
   
   if (toi.name == azimuthAxis) return geom.getAzimutAxis();
@@ -608,5 +791,5 @@
   if (i == pendulInfos.begin()) return true;
   i--;
-  return (sampleNum > (*i).second.SN);
+  return (sampleNum+4000> (*i).second.SN);
 }
 
@@ -627,4 +810,9 @@
 
 void StarMatcher::propagateLowBound(TOI const& toi, long sampleNum) {
+ // we want to keep some past information to interpolate...
+ // keep 1000 sampleNums (easier than a number of posinfos...)
+ 
+ sampleNum -= 4000;
+
   map<double, posInfo>::iterator i = posInfos.begin();
   while (i != posInfos.end() && (*i).first < sampleNum) i++;
