Changeset 577 in Sophya for trunk/Poubelle/archTOI.old


Ignore:
Timestamp:
Nov 16, 1999, 2:20:39 PM (26 years ago)
Author:
ansari
Message:

SST

Location:
trunk/Poubelle/archTOI.old
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/Poubelle/archTOI.old/covsrt.c

    r556 r577  
    11#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}
    22
    3 void covsrt(float **covar, int ma, int ia[], int mfit)
     3void covsrt(double **covar, int ma, int ia[], int mfit)
    44{
    55        int i,j,k;
    6         float swap;
     6        double swap;
    77
    88        for (i=mfit+1;i<=ma;i++)

  • trunk/Poubelle/archTOI.old/gaussj.c

    r556 r577  
    44#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
    55
    6 void gaussj(float **a, int n, float **b, int m)
     6void gaussj(double **a, int n, double **b, int m)
    77{
    88        int *indxc,*indxr,*ipiv;
    99        int i,icol,irow,j,k,l,ll;
    10         float big,dum,pivinv,temp;
     10        double big,dum,pivinv,temp;
    1111
    1212        indxc=ivector(1,n);

  • trunk/Poubelle/archTOI.old/gondolageom.cc

    r556 r577  
    6767}
    6868
    69 void GondolaGeom::solveStars() {
     69// $CHECK$ do a higher order fit ?
     70int GondolaGeom::solveStars() {
     71  if (nstars<2) return -1;
    7072  staz /= nstars;
    7173  st   /= nstars;
    … …  
    8183  if (azimut > 360) azimut -= 360;
    8284  if (azimut < 0)   azimut += 360;
     85
     86  return 0;
    8387}
    8488

  • trunk/Poubelle/archTOI.old/gondolageom.h

    r556 r577  
    1515  // Set orientation through SST
    1616  void addStar(double deltasn, double az, double elv, double diod);
    17   void solveStars();
     17  int solveStars();
    1818 
    1919  void getPointing(double elv, double az, double& trueElv, double& trueAz);

  • trunk/Poubelle/archTOI.old/lfit.c

    r566 r577  
    22#include "nrutil.h"
    33
    4 void lfit(float x[], float y[], float sig[], int ndat, float a[], int ia[],
    5         int ma, float **covar, float *chisq, void (*funcs)(float, float [], int))
     4void lfit(double x[], double y[], double sig[], int ndat, double a[], int ia[],
     5        int ma, double **covar, double *chisq, void (*funcs)(double, double [], int))
    66{
    7         void covsrt(float **covar, int ma, int ia[], int mfit);
    8         void gaussj(float **a, int n, float **b, int m);
     7        void covsrt(double **covar, int ma, int ia[], int mfit);
     8        void gaussj(double **a, int n, double **b, int m);
    99        int i,j,k,l,m,mfit=0;
    10         float ym,wt,sum,sig2i,**beta,*afunc;
     10        double ym,wt,sum,sig2i,**beta,*afunc;
    1111
    1212        beta=matrix(1,ma,1,1);

  • trunk/Poubelle/archTOI.old/nrutil.c

    r556 r577  
    1818}*/
    1919
    20 float *vector(long nl, long nh)
    21 /* allocate a float vector with subscript range v[nl..nh] */
    22 {
    23         float *v;
    24 
    25         v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
     20double *vector(long nl, long nh)
     21/* allocate a double vector with subscript range v[nl..nh] */
     22{
     23        double *v;
     24
     25        v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
    2626        if (!v) nrerror("allocation failure in vector()");
    2727        return v-nl+NR_END;
    … …  
    6868}
    6969
    70 float **matrix(long nrl, long nrh, long ncl, long nch)
    71 /* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
    72 {
    73         long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
    74         float **m;
    75 
    76         /* allocate pointers to rows */
    77         m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
    78         if (!m) nrerror("allocation failure 1 in matrix()");
    79         m += NR_END;
    80         m -= nrl;
    81 
    82         /* allocate rows and set pointers to them */
    83         m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
    84         if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
    85         m[nrl] += NR_END;
    86         m[nrl] -= ncl;
    87 
    88         for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
    89 
    90         /* return pointer to array of pointers to rows */
    91         return m;
    92 }
    93 
    94 double **dmatrix(long nrl, long nrh, long ncl, long nch)
     70double **matrix(long nrl, long nrh, long ncl, long nch)
    9571/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
    9672{
    … …  
    11692}
    11793
     94double **dmatrix(long nrl, long nrh, long ncl, long nch)
     95/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
     96{
     97        long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
     98        double **m;
     99
     100        /* allocate pointers to rows */
     101        m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
     102        if (!m) nrerror("allocation failure 1 in matrix()");
     103        m += NR_END;
     104        m -= nrl;
     105
     106        /* allocate rows and set pointers to them */
     107        m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
     108        if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
     109        m[nrl] += NR_END;
     110        m[nrl] -= ncl;
     111
     112        for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
     113
     114        /* return pointer to array of pointers to rows */
     115        return m;
     116}
     117
    118118int **imatrix(long nrl, long nrh, long ncl, long nch)
    119119/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
    … …  
    141141}
    142142
    143 float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
     143double **submatrix(double **a, long oldrl, long oldrh, long oldcl, long oldch,
    144144        long newrl, long newcl)
    145145/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
    146146{
    147147        long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
    148         float **m;
     148        double **m;
    149149
    150150        /* allocate array of pointers to rows */
    151         m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
     151        m=(double **) malloc((size_t) ((nrow+NR_END)*sizeof(double*)));
    152152        if (!m) nrerror("allocation failure in submatrix()");
    153153        m += NR_END;
    … …  
    161161}
    162162
    163 float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
    164 /* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
     163double **convert_matrix(double *a, long nrl, long nrh, long ncl, long nch)
     164/* allocate a double matrix m[nrl..nrh][ncl..nch] that points to the matrix
    165165declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
    166166and ncol=nch-ncl+1. The routine should be called with the address
    … …  
    168168{
    169169        long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
    170         float **m;
    171 
    172         /* allocate pointers to rows */
    173         m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
     170        double **m;
     171
     172        /* allocate pointers to rows */
     173        m=(double **) malloc((size_t) ((nrow+NR_END)*sizeof(double*)));
    174174        if (!m) nrerror("allocation failure in convert_matrix()");
    175175        m += NR_END;
    … …  
    183183}
    184184
    185 float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
    186 /* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
     185double ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
     186/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
    187187{
    188188        long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
    189         float ***t;
     189        double ***t;
    190190
    191191        /* allocate pointers to pointers to rows */
    192         t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
     192        t=(double ***) malloc((size_t)((nrow+NR_END)*sizeof(double**)));
    193193        if (!t) nrerror("allocation failure 1 in f3tensor()");
    194194        t += NR_END;
    … …  
    196196
    197197        /* allocate pointers to rows and set pointers to them */
    198         t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
     198        t[nrl]=(double **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double*)));
    199199        if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
    200200        t[nrl] += NR_END;
    … …  
    202202
    203203        /* allocate rows and set pointers to them */
    204         t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
     204        t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(double)));
    205205        if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
    206206        t[nrl][ncl] += NR_END;
    … …  
    218218}
    219219
    220 void free_vector(float *v, long nl, long nh)
    221 /* free a float vector allocated with vector() */
     220void free_vector(double *v, long nl, long nh)
     221/* free a double vector allocated with vector() */
    222222{
    223223        free((FREE_ARG) (v+nl-NR_END));
    … …  
    248248}
    249249
    250 void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
    251 /* free a float matrix allocated by matrix() */
     250void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
     251/* free a double matrix allocated by matrix() */
    252252{
    253253        free((FREE_ARG) (m[nrl]+ncl-NR_END));
    … …  
    269269}
    270270
    271 void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
     271void free_submatrix(double **b, long nrl, long nrh, long ncl, long nch)
    272272/* free a submatrix allocated by submatrix() */
    273273{
    … …  
    275275}
    276276
    277 void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
     277void free_convert_matrix(double **b, long nrl, long nrh, long ncl, long nch)
    278278/* free a matrix allocated by convert_matrix() */
    279279{
    … …  
    281281}
    282282
    283 void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
     283void free_f3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
    284284        long ndl, long ndh)
    285 /* free a float f3tensor allocated by f3tensor() */
     285/* free a double f3tensor allocated by f3tensor() */
    286286{
    287287        free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
    … …  
    309309}
    310310
    311 float *vector(nl,nh)
    312 long nh,nl;
    313 /* allocate a float vector with subscript range v[nl..nh] */
    314 {
    315         float *v;
    316 
    317         v=(float *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(float)));
     311double *vector(nl,nh)
     312long nh,nl;
     313/* allocate a double vector with subscript range v[nl..nh] */
     314{
     315        double *v;
     316
     317        v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double)));
    318318        if (!v) nrerror("allocation failure in vector()");
    319319        return v-nl+NR_END;
    … …  
    364364}
    365365
    366 float **matrix(nrl,nrh,ncl,nch)
    367 long nch,ncl,nrh,nrl;
    368 /* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
    369 {
    370         long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
    371         float **m;
    372 
    373         /* allocate pointers to rows */
    374         m=(float **) malloc((unsigned int)((nrow+NR_END)*sizeof(float*)));
    375         if (!m) nrerror("allocation failure 1 in matrix()");
    376         m += NR_END;
    377         m -= nrl;
    378 
    379         /* allocate rows and set pointers to them */
    380         m[nrl]=(float *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float)));
    381         if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
    382         m[nrl] += NR_END;
    383         m[nrl] -= ncl;
    384 
    385         for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
    386 
    387         /* return pointer to array of pointers to rows */
    388         return m;
    389 }
    390 
    391 double **dmatrix(nrl,nrh,ncl,nch)
     366double **matrix(nrl,nrh,ncl,nch)
    392367long nch,ncl,nrh,nrl;
    393368/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
    … …  
    414389}
    415390
     391double **dmatrix(nrl,nrh,ncl,nch)
     392long nch,ncl,nrh,nrl;
     393/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
     394{
     395        long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
     396        double **m;
     397
     398        /* allocate pointers to rows */
     399        m=(double **) malloc((unsigned int)((nrow+NR_END)*sizeof(double*)));
     400        if (!m) nrerror("allocation failure 1 in matrix()");
     401        m += NR_END;
     402        m -= nrl;
     403
     404        /* allocate rows and set pointers to them */
     405        m[nrl]=(double *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double)));
     406        if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
     407        m[nrl] += NR_END;
     408        m[nrl] -= ncl;
     409
     410        for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
     411
     412        /* return pointer to array of pointers to rows */
     413        return m;
     414}
     415
    416416int **imatrix(nrl,nrh,ncl,nch)
    417417long nch,ncl,nrh,nrl;
    … …  
    440440}
    441441
    442 float **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
    443 float **a;
     442double **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
     443double **a;
    444444long newcl,newrl,oldch,oldcl,oldrh,oldrl;
    445445/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
    446446{
    447447        long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
    448         float **m;
     448        double **m;
    449449
    450450        /* allocate array of pointers to rows */
    451         m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));
     451        m=(double **) malloc((unsigned int) ((nrow+NR_END)*sizeof(double*)));
    452452        if (!m) nrerror("allocation failure in submatrix()");
    453453        m += NR_END;
    … …  
    461461}
    462462
    463 float **convert_matrix(a,nrl,nrh,ncl,nch)
    464 float *a;
    465 long nch,ncl,nrh,nrl;
    466 /* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
     463double **convert_matrix(a,nrl,nrh,ncl,nch)
     464double *a;
     465long nch,ncl,nrh,nrl;
     466/* allocate a double matrix m[nrl..nrh][ncl..nch] that points to the matrix
    467467declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
    468468and ncol=nch-ncl+1. The routine should be called with the address
    … …  
    470470{
    471471        long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
    472         float **m;
    473 
    474         /* allocate pointers to rows */
    475         m=(float **) malloc((unsigned int) ((nrow+NR_END)*sizeof(float*)));
     472        double **m;
     473
     474        /* allocate pointers to rows */
     475        m=(double **) malloc((unsigned int) ((nrow+NR_END)*sizeof(double*)));
    476476        if (!m) nrerror("allocation failure in convert_matrix()");
    477477        m += NR_END;
    … …  
    485485}
    486486
    487 float ***f3tensor(nrl,nrh,ncl,nch,ndl,ndh)
     487double ***f3tensor(nrl,nrh,ncl,nch,ndl,ndh)
    488488long nch,ncl,ndh,ndl,nrh,nrl;
    489 /* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
     489/* allocate a double 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
    490490{
    491491        long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
    492         float ***t;
     492        double ***t;
    493493
    494494        /* allocate pointers to pointers to rows */
    495         t=(float ***) malloc((unsigned int)((nrow+NR_END)*sizeof(float**)));
     495        t=(double ***) malloc((unsigned int)((nrow+NR_END)*sizeof(double**)));
    496496        if (!t) nrerror("allocation failure 1 in f3tensor()");
    497497        t += NR_END;
    … …  
    499499
    500500        /* allocate pointers to rows and set pointers to them */
    501         t[nrl]=(float **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(float*)));
     501        t[nrl]=(double **) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(double*)));
    502502        if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
    503503        t[nrl] += NR_END;
    … …  
    505505
    506506        /* allocate rows and set pointers to them */
    507         t[nrl][ncl]=(float *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(float)));
     507        t[nrl][ncl]=(double *) malloc((unsigned int)((nrow*ncol*ndep+NR_END)*sizeof(double)));
    508508        if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
    509509        t[nrl][ncl] += NR_END;
    … …  
    522522
    523523void free_vector(v,nl,nh)
    524 float *v;
    525 long nh,nl;
    526 /* free a float vector allocated with vector() */
     524double *v;
     525long nh,nl;
     526/* free a double vector allocated with vector() */
    527527{
    528528        free((FREE_ARG) (v+nl-NR_END));
    … …  
    562562
    563563void free_matrix(m,nrl,nrh,ncl,nch)
    564 float **m;
    565 long nch,ncl,nrh,nrl;
    566 /* free a float matrix allocated by matrix() */
     564double **m;
     565long nch,ncl,nrh,nrl;
     566/* free a double matrix allocated by matrix() */
    567567{
    568568        free((FREE_ARG) (m[nrl]+ncl-NR_END));
    … …  
    589589
    590590void free_submatrix(b,nrl,nrh,ncl,nch)
    591 float **b;
     591double **b;
    592592long nch,ncl,nrh,nrl;
    593593/* free a submatrix allocated by submatrix() */
    … …  
    597597
    598598void free_convert_matrix(b,nrl,nrh,ncl,nch)
    599 float **b;
     599double **b;
    600600long nch,ncl,nrh,nrl;
    601601/* free a matrix allocated by convert_matrix() */
    … …  
    605605
    606606void free_f3tensor(t,nrl,nrh,ncl,nch,ndl,ndh)
    607 float ***t;
     607double ***t;
    608608long nch,ncl,ndh,ndl,nrh,nrl;
    609 /* free a float f3tensor allocated by f3tensor() */
     609/* free a double f3tensor allocated by f3tensor() */
    610610{
    611611        free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));

  • trunk/Poubelle/archTOI.old/nrutil.h

    r556 r577  
    22#define _NR_UTILS_H_
    33
    4 static float sqrarg;
     4static double sqrarg;
    55#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
    66
    … …  
    1616        (dminarg1) : (dminarg2))
    1717
    18 static float maxarg1,maxarg2;
     18static double maxarg1,maxarg2;
    1919#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
    2020        (maxarg1) : (maxarg2))
    2121
    22 static float minarg1,minarg2;
     22static double minarg1,minarg2;
    2323#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
    2424        (minarg1) : (minarg2))
    … …  
    4545
    4646void nrerror(char error_text[]);
    47 float *vector(long nl, long nh);
     47double *vector(long nl, long nh);
    4848int *ivector(long nl, long nh);
    4949unsigned char *cvector(long nl, long nh);
    5050unsigned long *lvector(long nl, long nh);
    5151double *dvector(long nl, long nh);
    52 float **matrix(long nrl, long nrh, long ncl, long nch);
     52double **matrix(long nrl, long nrh, long ncl, long nch);
    5353double **dmatrix(long nrl, long nrh, long ncl, long nch);
    5454int **imatrix(long nrl, long nrh, long ncl, long nch);
    55 float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
     55double **submatrix(double **a, long oldrl, long oldrh, long oldcl, long oldch,
    5656        long newrl, long newcl);
    57 float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch);
    58 float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
    59 void free_vector(float *v, long nl, long nh);
     57double **convert_matrix(double *a, long nrl, long nrh, long ncl, long nch);
     58double ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
     59void free_vector(double *v, long nl, long nh);
    6060void free_ivector(int *v, long nl, long nh);
    6161void free_cvector(unsigned char *v, long nl, long nh);
    6262void free_lvector(unsigned long *v, long nl, long nh);
    6363void free_dvector(double *v, long nl, long nh);
    64 void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
     64void free_matrix(double **m, long nrl, long nrh, long ncl, long nch);
    6565void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
    6666void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch);
    67 void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch);
    68 void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch);
    69 void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
     67void free_submatrix(double **b, long nrl, long nrh, long ncl, long nch);
     68void free_convert_matrix(double **b, long nrl, long nrh, long ncl, long nch);
     69void free_f3tensor(double ***t, long nrl, long nrh, long ncl, long nch,
    7070        long ndl, long ndh);
    7171
    … …  
    7474
    7575void nrerror();
    76 float *vector();
    77 float **matrix();
    78 float **submatrix();
    79 float **convert_matrix();
    80 float ***f3tensor();
     76double *vector();
     77double **matrix();
     78double **submatrix();
     79double **convert_matrix();
     80double ***f3tensor();
    8181double *dvector();
    8282double **dmatrix();

  • trunk/Poubelle/archTOI.old/starmatcher.cc

    r555 r577  
    88#include "archparam.h"
    99#include "gondolageom.h"
     10#include "polfitclip.h"
     11
     12#define STARDUMP
    1013
    1114extern "C" {
    … …  
    1417#include "nrutil.h"
    1518
    16 void lfit(float x[], float y[], float sig[], int ndat, float a[], int ia[],
    17         int ma, float **covar, float *chisq, void (*funcs)(float, float [], int));
    18 
    19 void polfunc(float x, float afunc[], int ma);
    20 void sinfunc(float x, float afunc[], int ma);
    21 }
    22 
    23 void polfunc(float x, float afunc[], int ma) {
     19void lfit(double x[], double y[], double sig[], int ndat, double a[], int ia[],
     20        int ma, double **covar, double *chisq, void (*funcs)(double, double [], int));
     21
     22void polfunc(double x, double afunc[], int ma);
     23void sinfunc(double x, double afunc[], int ma);
     24}
     25
     26void polfunc(double x, double afunc[], int ma) {
    2427  afunc[1] = 1;
    2528  for (int i=2; i<=ma; i++)
    … …  
    2730}
    2831
    29 void sinfunc(float x, float afunc[], int /*ma*/) {
     32void sinfunc(double x, double afunc[], int /*ma*/) {
    3033  afunc[1] = cos(x);
    3134  afunc[2] = sin(x);
    … …  
    3437
    3538
    36 float polval(float x, float a[], int ma);
    37 
    38 float polval(float x, float a[], int ma) {
    39   float r = a[ma];
     39double polval(double x, double a[], int ma);
     40
     41double polval(double x, double a[], int ma) {
     42  double r = a[ma];
    4043  for (int i=ma-1; i>0; i--) {
    4144    r = r*x+a[i];
    … …  
    186189static ofstream pendstream("pendul.dat");
    187190#endif
     191static ofstream logstream("starmatch.log");
    188192
    189193void StarMatcher::dataFeed(SSTEtoile const& x) {
    … …  
    295299 
    296300  // new set of matched stars... Clean, and get parameters...
    297   // We don't want more than 20 seconds of data
     301  // We don't want more than 30 seconds of data
    298302 
    299303  if (matchStars.empty()) return;
    … …  
    303307  deque<matchStar>::iterator it;
    304308  for (it = matchStars.begin();  it!=matchStars.end(); it++) {
    305      if ((snEnd - (*it).SN)*archParam.acq.perEch < 20)
     309     if ((snEnd - (*it).SN)*archParam.acq.perEch < 30 ||
     310         (*it).seq > lastCleanSave)
    306311       break;
     312  }
     313  if (it != matchStars.begin()) {
     314    it--;
    307315  }
    308316  if (it != matchStars.begin()) {
    … …  
    314322  int nskip=0;
    315323  for (it = matchStars.begin();  it!=matchStars.end(); it++) {
    316      if ((snEnd - (*it).SN)*archParam.acq.perEch < 7)
     324     if ((snEnd - (*it).SN)*archParam.acq.perEch < 5)
    317325       break;
    318326     nskip++;
    … …  
    334342    if (burstLen >= 4) {
    335343      for (deque<matchStar>::iterator it2=lastIt; it2 != it1; it2++) {
     344        //if ((*it2).ok)
     345        //  logstream << "  kill " << (*it2).seq << " " << setprecision(11) << (*it2).SN << " burst" << '\n';
    336346        (*it2).ok=false;
    337347      }
    … …  
    343353  // we fit the data to a polynomial, with clipping...
    344354 
    345   float* sn   =  ::vector(1, matchStars.size());
    346   float* elv0 =  ::vector(1, matchStars.size());
    347   float* azi  =  ::vector(1, matchStars.size());
    348   float* sig  =  ::vector(1, matchStars.size());
    349   float* ae   =  ::vector(1, 3);
    350   float* aa   =  ::vector(1, 3);
     355  //double* sn   =  ::vector(1, matchStars.size());
     356  double* elv0 =  ::vector(1, matchStars.size());
     357  double* azi  =  ::vector(1, matchStars.size());
     358  double* sig  =  ::vector(1, matchStars.size());
     359  //double* ae   =  ::vector(1, 3);
     360  double* aa   =  ::vector(1, 3);
    351361  int*  ia    = ivector(1, 3);
    352   float** cov =  matrix(1, 3, 1, 3);
     362  double** cov =  matrix(1, 3, 1, 3);
    353363  int ndata;
    354364 
    355   long sn0 = matchStars.front().SN;
    356   long snmin;
    357   long snmax;
    358   for (int i=0; i<4; i++) {
    359     ndata = 0;
    360     snmin = 99999999;
    361     snmax = -99999999;
    362     for (deque<matchStar>::iterator it1 = it ;  it1!=matchStars.end(); it1++) {
    363       matchStar s = (*it1);
    364       if (!s.ok) continue;
    365       double delv, daz;
    366       if (i) {
    367         delv = polval(s.SN-sn0, ae, 3)-(s.elvGSC - s.nDiode*1.41/45.);
    368         daz  = polval(s.SN-sn0, aa, 3)- s.azGSC;
    369         if (daz>=180) daz -= 360;
    370         if (daz<-180) daz += 360;
     365  //long sn0 = matchStars.front().SN;
     366  long sn0 = (*it).SN;
     367  PolFitClip2 fitElvAz(matchStars.size(), 2); fitElvAz.setClip(0.1,0,2,3);
     368  ndata = 0;
     369
     370  double oldAz = -1;
     371  for (deque<matchStar>::iterator it1 = it ;  it1!=matchStars.end(); it1++) {
     372    matchStar s1 = (*it1);
     373    if (!s1.ok) continue;
     374    double az = s1.azGSC;
     375    if (ndata > 0 && az - oldAz >  180) az -= 360;
     376    if (ndata > 0 && az - oldAz < -180) az += 360;
     377    fitElvAz.addData(s1.SN-sn0, s1.elvGSC - s1.nDiode*1.41/45., az);
     378    oldAz = az;
     379    ndata++;
     380  }
     381   
     382  double celv[3], caz[3];
     383  if (fitElvAz.doFit(celv,caz)) return;
     384  //if (fitElvAz.doFit()) return;
     385 
     386  //logstream << "*** Fit sig=" << fitElvAz.getSigmaY() << " " << fitElvAz.getSigmaZ()
     387  //          << " n  =" << fitElvAz.getNData() << " " << fitElvAz.getNDataUsed()
     388  //          << " SN :" << fitElvAz.getXMin() + sn0 << " - " << fitElvAz.getXMax() + sn0 << '\n';
     389  //logstream << " sn0 = " << sn0 << "; snmin =" << fitElvAz.getXMin() + sn0 << "; snmax ="
     390  //          << fitElvAz.getXMax() + sn0  << '\n';
     391  //logstream << " fitelv[x_] := " << celv[2] << " x^2 + " << celv[1] << " x + " << celv[0] << '\n';
     392  //logstream << " fitaz[x_]  := " <<  caz[2] << " x^2 + " <<  caz[1] << " x + " <<  caz[0] << '\n';
     393 
     394  //if (fitElvAz.getSigmaY() > 0.05 || fitElvAz.getSigmaZ() > 0.05) return;
     395  if (fitElvAz.getNDataUsed() < 5 ||
     396      (double)fitElvAz.getNDataUsed()/fitElvAz.getNData() < .5) return;
     397     
     398  double dcutElv = fitElvAz.getSigmaY()*3;
     399  double dcutAz  = fitElvAz.getSigmaZ()*3;
     400 
     401  if (dcutElv < 0.05) dcutElv = 0.05;
     402  if (dcutAz  < 0.05) dcutAz  = 0.05;
     403
     404  // don't kill borders of fit....
     405  //if (matchStars.end() - it > 6)
     406  //  for (deque<matchStar>::iterator it1 = it+3 ;  it1!=matchStars.end()-3; it1++) {
     407  for (deque<matchStar>::iterator it1 = it ;  it1!=matchStars.end(); it1++) {
     408      matchStar sss = (*it1);
     409      if (!sss.ok) continue;
     410      if (fabs(fitElvAz.valueY(sss.SN-sn0)-(sss.elvGSC - sss.nDiode*1.41/45.)) > dcutElv) {
     411        (*it1).ok = false;
     412        //logstream << "  kill " << sss.seq << " " << setprecision(11) << sss.SN << " "
     413        //          << fitElvAz.valueY(sss.SN-sn0)-(sss.elvGSC - sss.nDiode*1.41/45.) << '\n';
     414        continue;
    371415      }
    372       double dcutelv=0.2;
    373       double dcutaz =0.4;
    374       if (i>=2) {
    375          dcutelv = 0.05;
    376          dcutaz  = 0.1;
     416      double daz = fitElvAz.valueZ(sss.SN-sn0) - sss.azGSC;
     417      if (daz>=180) daz -= 360;
     418      if (daz<-180) daz += 360;
     419      if (fabs(daz) > dcutAz) (*it1).ok = false;
     420      if (!(*it1).ok) {
     421        //logstream << "  kill " << sss.seq << " " << setprecision(11) << sss.SN << " "
     422        //          << fitElvAz.valueY(sss.SN-sn0)-(sss.elvGSC - sss.nDiode*1.41/45.) << " "
     423        //          << daz << '\n';
    377424      }
    378       if (i>=3) {
    379          dcutelv = 0.02;
    380          dcutaz  = 0.03;
    381       }
    382       if (i == 0 || ((fabs(delv)<dcutelv && fabs(daz)<dcutaz) && i<3)) {
    383           ndata++;
    384           sn  [ndata] = s.SN-sn0;
    385           elv0[ndata] = s.elvGSC - s.nDiode*1.41/45.;
    386           azi [ndata] = s.azGSC; // $CHECK$  ou ajuster difference avec az galcross ?
    387           if (ndata>1 && azi[ndata] - azi[ndata-1] > 180) azi[ndata] -= 360;
    388           if (ndata>1 && azi[ndata] - azi[ndata-1] < -180) azi[ndata] += 360;
    389           sig [ndata] = 0.01;
    390           if (s.SN-sn0 > snmax) snmax = s.SN-sn0;
    391           if (s.SN-sn0 < snmin) snmin = s.SN-sn0;
    392       }
    393       if ((fabs(delv)>=dcutelv || fabs(daz)>=dcutaz) && i==3) {
    394           (*it1).ok = false;
    395       }
    396     }
    397     if (i==3) break;
    398     if (ndata<5) return; // on ne peut rien faire
    399     ia[1]   =  ia[2] = ia[3] = 1;
    400     float chi2;
    401     try{
    402       lfit(sn, elv0, sig, ndata, ae, ia, 3, cov, &chi2, polfunc);
    403       lfit(sn, azi,  sig, ndata, aa, ia, 3, cov, &chi2, polfunc);
    404     } catch(string s) {
    405       return;
    406     }
    407   }
    408  
    409   for (deque<matchStar>::iterator it1 = it ;  it1!=matchStars.end(); it1++) {
     425  }
     426   
     427  bool gotNewStars = false;
     428  for (deque<matchStar>::iterator it1 = matchStars.begin() ;  it1!=it; it1++) {
    410429    if ((*it1).ok && (*it1).seq > lastCleanSave) {
     430      gotNewStars = true;
    411431      lastCleanSave = (*it1).seq;
    412432#ifdef STARDUMP
    … …  
    425445  }
    426446 
     447  if (!gotNewStars) return;
     448 
    427449  // On a des etoiles nettoyees, on va trouver amplitude et phase du
    428450  // signal en elevation, ce qui va nous donner les deux angles d'Euler
    … …  
    430452 
    431453  // Il faut avoir une periode entiere ou pas loin, sinon on ne peut
    432   // rien dire simplement....
    433  
    434   it = matchStars.end(); it--;
    435   if ((((*it).SN) - (*matchStars.begin()).SN)*archParam.acq.perEch < 17) return;
    436  
    437   long snmid = (((*it).SN) + (*matchStars.begin()).SN)/2;
     454  // rien dire simplement.... -> we want to run on the last 18 seconds of
     455  // data before the last fully cleaned star (it).
     456 
     457  deque<matchStar>::iterator itstart;
     458 
     459  for (itstart = matchStars.begin();  itstart != it; itstart++) {
     460     if (((*it).SN - (*itstart).SN)*archParam.acq.perEch < 19)
     461       break;
     462  }
     463 
     464  if (((*it).SN - (*itstart).SN)*archParam.acq.perEch < 15) return;
     465
     466 
     467 // it = matchStars.end(); it--;
     468 // if (((*it).SN - matchStars.front().SN)*archParam.acq.perEch < 17) return;
     469 
     470  // $CHECK$ utiliser plutot le SN moyen/median de tous les points effectivement utilises.
     471  long snmid = ((*it).SN + (*itstart).SN)/2;
    438472   
    439473  ndata=0; 
    440  
    441   for (deque<matchStar>::iterator it1 = matchStars.begin() ;  it1!=matchStars.end(); it1++) {
    442     if (!(*it1).ok) continue;
     474  double snmean = 0;
     475
     476  logstream << "PendFit : " <<  setprecision(11) << (*itstart).SN << '-' << (*it).SN <<  " "
     477            << setprecision(4)
     478            << ((*it).SN - (*itstart).SN)*archParam.acq.perEch << " " ;
     479             
     480  for (deque<matchStar>::iterator it1 = itstart ;  it1!=it; it1++) {
     481    matchStar st = *it1;
     482    if (!st.ok) continue;
    443483    ndata++;
    444     azi[ndata]  = (*it1).azGSC * 3.1415926/180;
    445     elv0[ndata] = (*it1).elvGSC - (*it1).nDiode*1.41/45.;
     484    snmean += st.SN;
     485    azi[ndata]  = st.azGSC * 3.1415926/180;
     486    elv0[ndata] = st.elvGSC - st.nDiode*1.41/45.;
    446487    sig[ndata]  = 0.01;
    447488  }
    448    ia[1]   =  ia[2] = 1;
    449    ia[3] = 0;
    450    aa[3] = GondolaGeom::elevSST0;// do not fit elv0
     489  if (ndata) snmean /= ndata;
     490 
     491  ia[1]   =  ia[2] = 1;
     492  ia[3] = 0;
     493  aa[3] = GondolaGeom::elevSST0;// do not fit elv0
    451494
    452495  if (ndata<5) return;
    453   float chi2;
     496  double chi2;
    454497  try {
    455498    lfit(azi, elv0,  sig, ndata, aa, ia, 3, cov, &chi2, sinfunc);
    456   } catch(string s) {
     499  } catch(string st) {
    457500    return;
    458501  }
    459502
    460   double c = aa[1];
    461   double s = aa[2];
     503  double cc = aa[1];
     504  double ss = aa[2];
     505 
     506  logstream << setprecision(11) << snmean << setprecision(4)
     507            << " cs=" << cc << " " << ss << " chi2r=" << chi2/ndata
     508            << " cov " << cov[1][1] << " " << cov[2][2] << '\n';
    462509 
    463510  // Get rid of bad fits. The cuts are rather ad hoc
    464511 
    465512  //if (aa[3] < 39.64 || aa[3] > 39.68) return;
    466   if (chi2/ndata > 4) return;
     513  if (chi2/ndata > 9) return;
    467514  if (cov[1][1] > 0.0001) return;
    468515  if (cov[2][2] > 0.0001) return;
    469516
    470   double ampl = sqrt(c*c+s*s);
    471   double phase = atan2(c,s)/(3.1415926/180);
    472  
    473 #ifdef STARDUMP
    474   pendstream << snmid << " " << ampl << " " << phase << " " << ndata << " " << chi2/ndata << 
    475        " "  << cov[1][1] << " " << cov[2][2] << '\n';
    476 #endif
     517  double ampl = sqrt(cc*cc+ss*ss);
     518  double phase = atan2(cc,ss)/(3.1415926/180);
     519 
    477520 
    478521  pendulInfo info;
    479   info.SN          =   snmid;
    480   info.azPendul    =   phase;
     522  info.SN          =   snmean;
     523  info.azPendul    =   180-phase;
     524  if (info.azPendul > 360) info.azPendul -= 360;
     525  if (info.azPendul <   0) info.azPendul += 360;
    481526  info.angPendul   =   ampl;
    482527  pendulInfos[info.SN] = info;
    483528
     529#ifdef STARDUMP
     530  pendstream << setprecision(11) << snmean << " "
     531             << setprecision(4) << ampl << " " << info.azPendul << " " << ndata << " " << chi2/ndata <<  " "
     532             << cov[1][1] << " " << cov[2][2] << '\n';
     533#endif
    484534  /*
    485535  double snum = (matchStars.front().SN + matchStars.back().SN)/2-sn0;
    … …  
    497547//  }
    498548 
    499   free_vector(sn, 1, matchStars.size());
     549  //free_vector(sn, 1, matchStars.size());
    500550  free_vector(elv0, 1, matchStars.size());
    501551  free_vector(azi, 1, matchStars.size());
    502552  free_vector(sig, 1, matchStars.size());
    503   free_vector(ae, 1, 3);
     553  //free_vector(ae, 1, 3);
    504554  free_vector(aa, 1, 3);
    505555  free_ivector(ia, 1, matchStars.size());
    506556  free_matrix(cov, 1, 3, 1, 3);
    507557}
     558
     559
     560// $CHECK$ do a polynomial fit with several points...
     561int StarMatcher::getPendulInfo(double sampleNum, pendulInfo& info) {
     562
     563  PolFitClip2 fitPendul(30,2);
     564
     565  map<double, pendulInfo>::iterator i = pendulInfos.lower_bound(sampleNum);
     566  if (i == pendulInfos.begin() && (*i).second.SN >= sampleNum) return -1;
     567  if (i == pendulInfos.end()   && (*i).second.SN <= sampleNum) return -1;
     568 
     569  if ((*i).second.SN > sampleNum) i--;
     570
     571  int nn=0;
     572  double aziprev=0, azicur=0, azi0=0;
     573  for (map<double, pendulInfo>::iterator ii=i; ii != pendulInfos.begin(); ii--) {
     574    nn++;
     575    pendulInfo inf1 = (*ii).second;
     576    aziprev = azicur;
     577    azicur = inf1.azPendul;
     578    if (nn==1) azi0 = azicur;
     579    if (nn>1 && azicur - aziprev > 180)  azicur -= 360;
     580    if (nn>1 && azicur - aziprev < -180) azicur += 360;
     581    fitPendul.addData(inf1.SN, inf1.angPendul, azicur);
     582    if (nn>=5) break;
     583  }
     584 
     585  azicur = azi0;
     586  if (i != pendulInfos.end()) i++;
     587  for (map<double, pendulInfo>::iterator ii=i; ii != pendulInfos.end(); ii++) {
     588    nn++;
     589    pendulInfo inf1 = (*ii).second;
     590    aziprev = azicur;
     591    azicur = inf1.azPendul;
     592    if (nn>1 && azicur - aziprev > 180)  azicur -= 360;
     593    if (nn>1 && azicur - aziprev < -180) azicur += 360;
     594    fitPendul.addData(inf1.SN, inf1.angPendul, azicur);
     595    if (nn>=10) break;
     596  }
     597 
     598  if (fitPendul.doFit()) return -1;
     599 
     600  info.SN = sampleNum;
     601  info.azPendul   = fitPendul.valueZ(sampleNum);
     602  if (info.azPendul > 360) info.azPendul -= 360;
     603  if (info.azPendul <   0) info.azPendul += 360;
     604  info.angPendul  = fitPendul.valueY(sampleNum);
     605  return 0;
     606}
     607
     608#if 0
     609
     610int StarMatcher::getPendulInfo(double sampleNum, pendulInfo& info) {
     611  static double* sn    = ::vector(1, 100);
     612  static double* aziP  = ::vector(1, 100);
     613  static double* ampP  = ::vector(1, 100);
     614  static double* sig   = ::vector(1, 100);
     615  static double* aAzi  = ::vector(1, 3);
     616  static double* aAmp  = ::vector(1, 3);
     617  static int*   ia    = ::ivector(1,3);
     618  static double** cov  = ::matrix(1,3,1,3);
     619  int ndata = 0;
     620  map<double, pendulInfo>::iterator i = pendulInfos.lower_bound(sampleNum);
     621  if (i == pendulInfos.begin() && (*i).second.SN >= sampleNum) return -1;
     622  if (i == pendulInfos.end()   && (*i).second.SN <= sampleNum) return -1;
     623 
     624  if ((*i).second.SN > sampleNum) i--;
     625
     626  int nn=0;
     627  for (map<double, pendulInfo>::iterator ii=i; ii != pendulInfos.begin(); ii--) {
     628    nn++;
     629    ndata++;
     630    pendulInfo inf1 = (*ii).second;
     631    sn[ndata]   = inf1.SN;
     632    ampP[ndata] = inf1.angPendul;
     633    aziP[ndata] = inf1.azPendul;
     634    int prev = ndata-1;
     635    if (ndata>1 && aziP[ndata] - aziP[prev] > 180)  aziP[ndata] -= 360;
     636    if (ndata>1 && aziP[ndata] - aziP[prev] < -180) aziP[ndata] += 360;
     637    sig[ndata]  = 1;
     638    if (nn>=50) break;
     639  }
     640 
     641  nn=0;
     642  if (i != pendulInfos.end()) i++;
     643  for (map<double, pendulInfo>::iterator ii=i; ii != pendulInfos.end(); ii++) {
     644    nn++;
     645    ndata++;
     646    pendulInfo inf1 = (*ii).second;
     647    sn[ndata]   = inf1.SN;
     648    ampP[ndata] = inf1.angPendul;
     649    aziP[ndata] = inf1.azPendul;
     650    int prev = ndata-1;
     651    if (nn==1) prev=1;
     652    if (ndata>1 && aziP[ndata] - aziP[prev] > 180)  aziP[ndata] -= 360;
     653    if (ndata>1 && aziP[ndata] - aziP[prev] < -180) aziP[ndata] += 360;
     654    sig[ndata]  = 1;
     655    if (nn>=50) break;
     656  }
     657 
     658  if (ndata < 3) return -1;
     659 
     660  ia[1] = ia[2] = ia[3] = 1;
     661  double chi2;
     662  try {
     663    lfit(sn, aziP,  sig, ndata, aAzi, ia, 3, cov, &chi2, polfunc);
     664    lfit(sn, ampP,  sig, ndata, aAmp, ia, 3, cov, &chi2, polfunc);
     665  } catch(string st) {
     666    return -1;
     667  }
     668
     669  info.SN = sampleNum;
     670  info.azPendul   = polval(sampleNum, aAzi, 3);
     671  if (info.azPendul > 360) info.azPendul -= 360;
     672  if (info.azPendul <   0) info.azPendul += 360;
     673  info.angPendul  = polval(sampleNum, aAmp, 3);
     674  return 0;
     675}
     676
     677#endif
     678
     679#if 0
    508680
    509681int StarMatcher::getPendulInfo(double sampleNum, pendulInfo& info) {
    … …  
    527699}
    528700
     701#endif
     702
    529703
    530704double StarMatcher::getValue(long sampleNum, TOI const& toi) {
    … …  
    546720  if ((*i).second.SN > sampleNum) i--;
    547721 
     722  // $CHECK$ if time spent here, can keep a GondolaGeom object for several
     723  // samples...
    548724  GondolaGeom geom;
    549725  geom.setEarthPos((*i).second.lon,(*i).second.lat);
    550726  geom.setTSid((*i).second.ts);
    551  
     727  geom.setPendulation(pendul.azPendul, pendul.angPendul);
     728 
     729  int ns=0;
    552730  for (map<double, posInfo>::iterator it=i; it != posInfos.end(); it++) {
    553731    posInfo s = (*it).second;
    554732    double delsn = s.SN - sampleNum;
    555     if (delsn * archParam.acq.perEch > 1) break;
     733    ns++;
     734    //if (delsn * archParam.acq.perEch > 1 && ns > 4) break;
     735    if (delsn * archParam.acq.perEch > 5) break;
    556736    geom.addStar(delsn, s.azStar, s.elvStar, s.diodStar);
    557737  }
    558738 
    559739  if (i != posInfos.begin()) i--;
     740  ns = 0;
    560741  for (map<double, posInfo>::iterator it=i; it != posInfos.begin(); it--) {
    561742    posInfo s = (*it).second;
    562743    double delsn = s.SN - sampleNum;
    563     if (-delsn * archParam.acq.perEch > 1) break;
     744    ns++;
     745    //if (-delsn * archParam.acq.perEch > 1 && ns > 4) break;
     746    if (-delsn * archParam.acq.perEch > 5) break;
    564747    geom.addStar(delsn, s.azStar, s.elvStar, s.diodStar);
    565748  }
    566749
    567   geom.solveStars();
     750  if (geom.solveStars()) return -99999;
    568751 
    569752  if (toi.name == azimuthAxis) return geom.getAzimutAxis();
    … …  
    608791  if (i == pendulInfos.begin()) return true;
    609792  i--;
    610   return (sampleNum > (*i).second.SN);
     793  return (sampleNum+4000> (*i).second.SN);
    611794}
    612795
    … …  
    627810
    628811void StarMatcher::propagateLowBound(TOI const& toi, long sampleNum) {
     812 // we want to keep some past information to interpolate...
     813 // keep 1000 sampleNums (easier than a number of posinfos...)
     814 
     815 sampleNum -= 4000;
     816
    629817  map<double, posInfo>::iterator i = posInfos.begin();
    630818  while (i != posInfos.end() && (*i).first < sampleNum) i++;
Note: See TracChangeset for help on using the changeset viewer.