Main Page | Modules | Data Structures | File List | Data Fields | Globals | Related Pages

calib/source/wcalib.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00018 /* =========================================================== */
00019 
00020 #include <IFU_io.h>
00021 #include <IFU_math.h>
00022 #include <snifs.h>
00023 #include <gsl/gsl_statistics.h>
00024 #include <gsl/gsl_math.h>            /* gsl_isnan() */
00025 #include <gsl/gsl_chebyshev.h>
00026 
00027 #define SWAP(a,b,tmp) (tmp)=(a); (a)=(b); (b)=(tmp);
00028 
00029 #define SIGINIT 2.                  
00030 
00031 #define CRITER_SIG_1    0.1
00032 #define CRITER_SIG_2    0.2
00033 #define CRITER_NREJ_1   0.2
00034 #define CRITER_NREJ_2   0.5
00035 #define CRITER_DIFF_1   3
00036 #define CRITER_DIFF_2   2            
00037 
00038 /* global variable needed for MultiGauss */
00039 
00040 int glnlines;
00041 double glLbdaStart, glLbdaStep;
00042 double *glArray;
00043 
00044 int glncon;
00045 
00048 typedef struct ArcLine {
00049 
00050   double          lambda;    
00051   struct ArcLine *next;      
00052 
00053 } ArcLine;
00054 
00057 typedef struct ArcGroup {
00058 
00059   int              nlines;            
00060   double           lbda_inf,lbda_sup; 
00061   ArcLine         *lines;             
00062   struct ArcGroup *next;              
00063 
00064 } ArcGroup;
00065 
00066 /* === Doxygen Comment ======================================= */
00073 /* =========================================================== */
00074 
00075 void addArcLine(ArcLine *list[], double lambda)
00076 {
00077   ArcLine *new,*courant,*precedent;
00078 
00079   /* Create the new line */
00080   new = (ArcLine *)malloc(sizeof(ArcLine));
00081   new->lambda = lambda;
00082   new->next   = NULL;
00083 
00084   if (*list == NULL) /* Initialize a new list */
00085     *list = new;
00086   else {
00087     /* Update list: Add reference to new line in preceeding line */
00088     courant = precedent = *list;
00089     while (courant != NULL) { /* Loop over the lines */
00090       precedent = courant;
00091       courant   = courant->next;
00092     }
00093     precedent->next = new;
00094   }
00095 }
00096 
00097 /* === Doxygen Comment ======================================= */
00103 /* =========================================================== */
00104 
00105 void delArcLines(ArcLine list[])
00106 {
00107   /* Recursive deletion */
00108   if (list->next != NULL) delArcLines(list->next);
00109   free(list);
00110 }
00111 
00112 /* === Doxygen Comment ======================================= */
00119 /* =========================================================== */
00120 
00121 ArcGroup *newArcGroup(ArcGroup *list[], double lambda)
00122 {
00123   ArcGroup *new,*courant,*precedent;
00124 
00125   /* Create a brand new group with the single line */
00126   new = (ArcGroup *)malloc(sizeof(ArcGroup));
00127   new->nlines   = 1;
00128   new->lbda_inf = 0;
00129   new->lbda_sup = 0;
00130   new->lines    = NULL;
00131   new->next     = NULL;
00132   addArcLine(&(new->lines),lambda);
00133 
00134   if (*list == NULL) /* Initialize a new list */
00135     *list = new;
00136   else {
00137     /* Update list: Add reference to new group in preceeding group */
00138     courant = precedent = *list;
00139     while (courant != NULL) { /* Loop over the groups */
00140       precedent = courant;
00141       courant   = courant->next;
00142     }
00143     precedent->next = new;
00144   }
00145 
00146   return(new);
00147 }
00148 
00149 /* === Doxygen Comment ======================================= */
00156 /* =========================================================== */
00157 
00158 void expandArcGroup(ArcGroup *group, double lambda)
00159 { 
00160   addArcLine(&(group->lines),lambda);
00161   group->nlines++;
00162 }
00163 
00164 /* === Doxygen Comment ======================================= */
00170 /* =========================================================== */
00171 
00172 void delArcGroups(ArcGroup list[])
00173 {
00174   /* Recursive deletion */
00175   if (list->next != NULL) delArcGroups(list->next);
00176   delArcLines(list->lines);
00177   free(list);
00178 }
00179 
00180 /* === Doxygen Comment ======================================= */
00188 /* =========================================================== */
00189 
00190 void boundArcGroup(ArcGroup *group, double halfwidth, double lstep)
00191 {
00192   ArcLine *courant, *precedent;
00193 
00194   courant = precedent = group->lines;  
00195   group->lbda_inf = courant->lambda - (halfwidth*lstep);
00196 
00197   while (courant != NULL) {                  /* Loop over lines of group */
00198     precedent = courant;
00199     courant = courant->next;
00200   }
00201   group->lbda_sup = precedent->lambda + (halfwidth*lstep);
00202 }
00203 
00204 /* === Doxygen Comment ======================================= */
00214 /* =========================================================== */
00215 
00216 ArcGroup *findArcGroups(int nlbda, double *lbda, double halfwidth, double lstep)
00217 {
00218   int i;
00219   ArcGroup *group=NULL, *last=NULL;
00220 
00221   if (nlbda==1) last = newArcGroup(&group,lbda[0]);
00222 
00223   for (i=1; i<nlbda; i++) {
00224     if (lbda[i]-lbda[i-1] < 2*halfwidth*lstep) {
00225       if (last==NULL) {
00226         last = newArcGroup(&group,lbda[i-1]);
00227         expandArcGroup(last,lbda[i]);
00228       }
00229       else 
00230         expandArcGroup(last,lbda[i]);
00231     }
00232     else {
00233       if (last==NULL) {
00234         last = newArcGroup(&group,lbda[i-1]);
00235         boundArcGroup(last,halfwidth,lstep);
00236         last = newArcGroup(&group,lbda[i]);
00237         boundArcGroup(last,halfwidth,lstep);
00238       }
00239       else {
00240         boundArcGroup(last,halfwidth,lstep); 
00241         last = newArcGroup(&group,lbda[i]);
00242       }
00243     }
00244   }
00245   boundArcGroup(last,halfwidth,lstep); 
00246 
00247   return(group);
00248 }
00249 
00250 /* === Doxygen Comment ======================================= */
00260 /* =========================================================== */
00261 
00262 void setGslCheb(gsl_cheb_series *cheb,
00263                 double *coef, long ncoef, double xmin, double xmax)
00264 {
00265   memcpy(cheb->c, coef, ncoef*sizeof(double));
00266   cheb->order = ncoef-1;   /* BEWARE order vs. ncoef */
00267   cheb->a = xmin;
00268   cheb->b = xmax;
00269 }
00270 
00271 /* === Doxygen Comment ======================================= */
00281 /* =========================================================== */
00282 
00283 int fillSpectrumArray(SPECTRUM *spectre, double linf, double lsup,
00284                       gsl_cheb_series *continuum)
00285 {
00286   int ninf,nsup,nb,i,j;
00287   double *array, lbda;
00288 
00289   ninf = MAX(pixel_spec(spectre,linf),0);
00290   nsup = MIN(pixel_spec(spectre,lsup),spectre->npts-1);
00291   glLbdaStart = coord_spec(spectre, ninf);
00292   if ((nb = nsup-ninf + 1) < 5) return(nb);  /* not between bounds */
00293 
00294   array = (double *)malloc(nb*sizeof(double));
00295   for (i=ninf,j=0; i<=nsup; i++,j++) {
00296     array[j] = RD_spec(spectre,i);
00297     lbda = coord_spec(spectre,i);
00298     if (lbda >= continuum->a && lbda <= continuum->b) 
00299       array[j] -= gsl_cheb_eval(continuum, lbda);
00300   }
00301   glArray = array; 
00302   return(nb);
00303 }
00304 
00305 /* === Doxygen Comment ======================================= */
00324 /* =========================================================== */
00325 
00326 /* === Doxygen Comment ======================================= */
00339 /* =========================================================== */
00340 int nllsqfit_constrained(int npar, int npts, int nconst, 
00341                          double *par, double *bl, double *bu, 
00342                          void (*objfun)(), double *rnorm)
00343 {
00344   int dqed();
00345   int nvars = npar;             /* Parameters of the adjusted model */
00346   int mequa = npts;             /* Points to be adjusted */
00347   int mcon = nconst;            /* Linear constraints, not including bounds */
00348 
00349   int blflag, buflag;
00350   int i, *ind, ldfj, wsizes[2], *iw, iopt[24], status;
00351   double *fj, *w, ropt[4];
00352 
00353   if ((nvars+mcon) >= mequa) return(-1);
00354 
00355   /* Simple bounds and linear constraints */
00356   ind = (int *)malloc((nvars+mcon)*sizeof(int));
00357   for (i=0; i<(nvars+mcon); i++) {
00358     blflag = (bl != NULL && bl[i] != -MAXDOUBLE);
00359     buflag = (bu != NULL && bu[i] !=  MAXDOUBLE);
00360     if (blflag && buflag) ind[i] = 3;        /* Both bounds */
00361     else {
00362       if      (blflag) ind[i] = 1;           /* Only lower bound */
00363       else if (buflag) ind[i] = 2;           /* Only upper bound */
00364       else             ind[i] = 4;           /* No bound */
00365     }
00366   }
00367 
00368   ldfj = mcon+mequa+1;
00369   fj = (double *)malloc(ldfj*(nvars+1)*sizeof(double)); /* f & deriv. */  
00370 
00371   /* Options */
00372   iopt[0] = 2;  iopt[1] = 150;               /* Max. iterations = 150 */
00373   iopt[2] = 14; iopt[3] = 1;              /* Suppress quadratic model */
00374   iopt[4] = 1;  iopt[5] = 0;           /* Set ot 1 for Mega-verbosity */
00375   iopt[6] = 17; iopt[7] = 1;          /* Stop after a full model step */
00376   iopt[8] = 99;                                     /* End of options */
00377 
00378   /* Blanck call to dqed to get appropriate working array sizes */
00379   wsizes[0] = 1;
00380   wsizes[1] = 1;
00381   dqed(objfun,&mequa,&nvars,&mcon,ind,bl,bu,par,fj,&ldfj,rnorm,
00382        &status,iopt,ropt,wsizes,NULL);
00383 
00384   /* Allocate working arrays according to dqed needs */
00385   w  = (double *)malloc(wsizes[0]*sizeof(double));
00386   iw =    (int *)malloc(wsizes[1]*sizeof(int));
00387   iw[0] = wsizes[0];
00388   iw[1] = wsizes[1];
00389 
00390   dqed(objfun,&mequa,&nvars,&mcon,ind,bl,bu,par,fj,&ldfj,rnorm,
00391        &status,iopt,ropt,iw,w);
00392   
00393   free(ind); free(fj); free(iw); free(w);
00394   //if (status) print_msg("DQED status: %d",status);
00395   if (status >=2 && status <=7) status = 0;
00396   return(status);
00397 }
00398 
00399 /* === Doxygen Comment ======================================= */
00407 /* =========================================================== */
00408 
00409 void MultiGauss(double par[], double fj[], int *ldfj, 
00410                 int *igo, int iopt[], double ropt[])
00411 {
00412   int i,k,l, npts=*ldfj-1;
00413   double lbda,fsum,dx,dx2,expo,f;
00414 
00415   for (lbda=glLbdaStart, i=0; i<npts; i++, lbda+=glLbdaStep) {
00416     for (fsum=k=l=0; k<glnlines; k++, l+=3) {
00417       dx = (par[l+1]-lbda)/par[l+2];
00418       dx2 = dx*dx;
00419       expo = exp(-0.5*dx2);
00420       f = par[l]*expo;
00421       fsum += f;
00422       fj[  l  *(*ldfj)+i] = expo;              /* df/dIk */
00423       fj[(l+1)*(*ldfj)+i] = -f*dx/par[l+2];    /* df/dLk */
00424       fj[(l+2)*(*ldfj)+i] =  f*dx2/par[l+2];   /* df/dSk */
00425     }
00426     fj[(glnlines*3)**ldfj+i] = fsum - glArray[i]; /* valeur de f - measured */
00427   }
00428 }
00429 
00430 /* === Doxygen Comment ======================================= */
00444 /* =========================================================== */
00445 
00446 void fitLineSet(long NbPixel, double Wbound, double *refLbda, 
00447                 double *fitLbda, double *fitSig, double *fitInt, int *fitFlag,
00448                 double sigmin,double sigmax)
00449 { 
00450   long npar;
00451   int i, j, k, ifail, sorted = 1;
00452   double *par, *bl, *bu, residual, lbdaend, tmp;
00453 
00454   npar = 3*glnlines;                         /* I,L,S per gaussian */
00455   par= (double *)malloc(npar*sizeof(double));
00456   bl = (double *)malloc(npar*sizeof(double)); 
00457   bu = (double *)malloc(npar*sizeof(double)); 
00458 
00459   lbdaend = glLbdaStart + (NbPixel-1) * glLbdaStep;
00460   
00461   for (i=k=0; i<glnlines; i++, k+=3) {
00462     j = (int)((refLbda[i]-glLbdaStart)/glLbdaStep + 0.5);   
00463     par[k]   = ABS(glArray[j]);              /* Initial intensity */
00464     par[k+1] = refLbda[i];                   /* Initial wavelength [AA] */
00465     par[k+2] = SIGINIT;                      /* Initial sigma [AA] */
00466     /* Simple bounds */
00467     bl[k]   = 0.;          bu[k]   = MAXDOUBLE;
00468     bl[k+1] = glLbdaStart; bu[k+1] = lbdaend;
00469     bl[k+2] = sigmin;      bu[k+2] = sigmax;
00470   }
00471 
00472   ifail = nllsqfit_constrained(npar,NbPixel,0,par,bl,bu,
00473                                MultiGauss,&residual);
00474 
00475   for (i=k=0; i<glnlines; i++, k+=3) { 
00476     fitInt[i]  = par[k];        /* Final intensity */
00477     fitLbda[i] = par[k+1];      /* Adjusted wavelength */
00478     fitSig[i]  = par[k+2];      /* Adjusted sigma */
00479   }
00480   
00481   /* Sort the lines if needed (no one enforced X1 < X2 < X3... yet) */
00482   for (i=0; i<glnlines-1; i++) 
00483     sorted = sorted && (fitLbda[i] <= fitLbda[i+1]);
00484   while(!sorted) {
00485     if (DEBUG) print_msg("Sorting line group %.1f-%.1f", 
00486                          refLbda[0], refLbda[glnlines-1]);
00487     for (sorted=1, i=0; i<glnlines-1; i++) {
00488       if (fitLbda[i] > fitLbda[i+1]) {                   /* Swap elements */
00489         SWAP(fitLbda[i],fitLbda[i+1],tmp);
00490         SWAP(fitSig[i], fitSig[i+1], tmp);
00491         SWAP(fitInt[i], fitInt[i+1], tmp);
00492       }
00493       sorted = sorted && (fitLbda[i]<=fitLbda[i+1]);     /* Check sort */
00494     }
00495   }
00496 
00497   /* Adjustment flag */
00498   for (i=k=0; i<glnlines; i++, k+=3) { 
00499     if (fitInt[i] == 0) 
00500       fitFlag[i] = 98;          /* Discard if intensity=0 */
00501     else if ( Wbound > 0 && 
00502               ( (fitLbda[i] - Wbound*fitSig[i]) < glLbdaStart || 
00503                 (fitLbda[i] + Wbound*fitSig[i]) > lbdaend ) ) 
00504       fitFlag[i] = 99;          /* Discard if too close to the edges */
00505     else 
00506       fitFlag[i] = ifail;
00507   }
00508   free(bl); free(bu); free(par);
00509 }
00510 
00511 /* === Doxygen Comment ======================================= */
00521 /* =========================================================== */
00522 
00523 void MultiGauss_constrained(double par[], double fj[], int *ldfj, 
00524                             int *igo, int iopt[], double ropt[])
00525 {
00526   int ncon = glncon;
00527   int npts = *ldfj-1-ncon;
00528   int i,k,l,p;
00529   double lbda,fsum,dx,dx2,expo,f;
00530 
00531   for (i=0; i<ncon; i++) {
00532     /* Constraint derivatives */
00533     for (k=0; k<glnlines; k++) { 
00534       l = 3*k;
00535       fj[  l  *(*ldfj)+i] = 0;  /* No constraint on line intensity */
00536       if (k==i+1)         /* Constraint on fitLbda[k+1]-fitLbda[k] */
00537         fj[(l+1)*(*ldfj)+i] = 1;
00538       else if (k==i)
00539         fj[(l+1)*(*ldfj)+i] = -1;
00540       else
00541         fj[(l+1)*(*ldfj)+i] = 0;
00542       fj[(l+2)*(*ldfj)+i] = 0;  /* No constraint on line sigma */
00543     }
00544     /* Constraints */
00545     fj[(glnlines*3)*(*ldfj)+i] = par[3*(i+1)+1] - par[3*i+1];
00546   }
00547 
00548   for (p=0; p<npts; p++) {
00549     i = p + ncon;
00550     /* Objective function derivatives */
00551     lbda = glLbdaStart + p*glLbdaStep;
00552     for (fsum=k=0; k<glnlines; k++) {         /* k: line index */
00553       l = 3*k;                                /* l: parameter index */
00554       dx = (par[l+1]-lbda)/par[l+2];
00555       dx2 = dx*dx;
00556       expo = exp(-0.5*dx2);
00557       f = par[l]*expo;
00558       fsum += f;
00559       fj[  l  *(*ldfj)+i] = expo;              /* df/dIk */
00560       fj[(l+1)*(*ldfj)+i] =-f*dx/par[l+2];     /* df/dLk */
00561       fj[(l+2)*(*ldfj)+i] = f*dx2/par[l+2];    /* df/dSk */
00562     }
00563     /* Objective function */
00564     fj[(glnlines*3)*(*ldfj)+i] = fsum - glArray[p]; /* model - measured */
00565   }
00566 }
00567 
00568 /* === Doxygen Comment ======================================= */
00589 /* =========================================================== */
00590 
00591 void fitLineSet_constrained(long NbPixel, double Wbound, 
00592                             double *refLbda, double *fitLbda, 
00593                             double *fitSig, double *fitInt, int *fitFlag,
00594                             double sigmin,double sigmax)
00595 { 
00596   long npar, ncon;
00597   int i, j, k, ifail, sorted = 1;
00598   double *par, *bl, *bu, residual, lbdaend;
00599 
00600   npar = 3*glnlines;                        /* I,L,S per gaussian */
00601   ncon = glnlines-1;                        /* Nb of linear constraints */
00602   glncon = ncon;
00603 
00604   par = (double *)malloc(npar*sizeof(double));
00605   bl  = (double *)malloc((npar+ncon)*sizeof(double)); 
00606   bu  = (double *)malloc((npar+ncon)*sizeof(double)); 
00607 
00608   lbdaend = glLbdaStart + (NbPixel-1) * glLbdaStep;
00609   
00610   for (i=k=0; i<glnlines; i++, k+=3) {
00611     j = (int)((refLbda[i]-glLbdaStart)/glLbdaStep + 0.5);   
00612     par[k]   = ABS(glArray[j]);              /* Initial intensity */
00613     par[k+1] = refLbda[i];                   /* Initial wavelength [AA] */
00614     par[k+2] = SIGINIT;                      /* Initial sigma [AA] */
00615     /* Simple bounds */
00616     bl[k]   = 0.;          bu[k]   = MAXDOUBLE;
00617     bl[k+1] = glLbdaStart; bu[k+1] = lbdaend;
00618     bl[k+2] = sigmin;      bu[k+2] = sigmax;
00619   }
00620 
00621   /* Constraints */
00622   for (i=0; i<ncon; i++, k++) {
00623     /* Simple wavelength ordering constraint */
00624     //bl[k] = 0; bu[k] = MAXDOUBLE;
00625     /* Line separation in 50-150% of reference separation */
00626     bl[k] = 0.5*(refLbda[i+1]-refLbda[i]);
00627     bu[k] = 1.5*(refLbda[i+1]-refLbda[i]);
00628   }
00629 
00630   ifail = nllsqfit_constrained(npar,NbPixel,ncon,par,bl,bu,
00631                                MultiGauss_constrained,&residual);
00632 
00633   for (i=k=0; i<glnlines; i++, k+=3) { 
00634     fitInt[i]  = par[k];        /* Final intensity */
00635     fitLbda[i] = par[k+1];      /* Adjusted wavelength */
00636     fitSig[i]  = par[k+2];      /* Adjusted sigma */
00637   }
00638   
00639   /* Check the line ordering */
00640   for (i=0; i<glnlines-1; i++) sorted = sorted && (fitLbda[i]<=fitLbda[i+1]);
00641   if (!sorted && DEBUG) 
00642     print_warning("Argh!!! Line group %.1f-%.1f is not sorted!!!",
00643                   refLbda[0],refLbda[glnlines-1]);
00644 
00645   /* Adjustment flag */
00646   for (i=k=0; i<glnlines; i++, k+=3) { 
00647     if (!sorted) 
00648       fitFlag[i] = 97;          /* Discard lines from unsorted groups */
00649     else if (fitInt[i] == 0) 
00650       fitFlag[i] = 98;          /* Discard if intensity=0 */
00651     else if ( Wbound > 0 && 
00652               ( (fitLbda[i] - Wbound*fitSig[i]) < glLbdaStart || 
00653                 (fitLbda[i] + Wbound*fitSig[i]) > lbdaend ) ) 
00654       fitFlag[i] = 99;          /* Discard if too close to the edges */
00655     else 
00656       fitFlag[i] = ifail;
00657   }
00658   free(bl); free(bu); free(par);
00659 }
00660 
00661 /* === Doxygen Comment ======================================= */
00672 /* =========================================================== */
00673 
00674 void saveGauss(SPECTRUM *spectre, 
00675                double *fitLbda, double *fitSig, double *fitInt, 
00676                int ng, gsl_cheb_series *continuum)
00677 {
00678   int i,k;
00679   double lbda,dx,fsum;
00680 
00681   for (i=0; i<spectre->npts; i++) {
00682     lbda = coord_spec(spectre, i);
00683     for (fsum=k=0; k<ng; k++) {
00684       dx = (fitLbda[k]-lbda)/fitSig[k];
00685       fsum += fitInt[k] * exp(-0.5*dx*dx);
00686     }
00687     if (lbda >= continuum->a && lbda <= continuum->b) 
00688       fsum += gsl_cheb_eval(continuum, lbda);
00689 
00690     WR_spec(spectre,i,fsum);
00691   }
00692 }
00693 
00694 /* === Doxygen Comment ======================================= */
00701 /* =========================================================== */
00702 
00703 void saveCalib(SPECTRUM *spectre, gsl_cheb_series *wcalib)
00704 {
00705   int i;
00706 
00707   for (i=0; i<spectre->npts; i++) 
00708     WR_spec(spectre,i, gsl_cheb_eval(wcalib, coord_spec(spectre, i)));
00709 }
00710 
00711 /* === Doxygen Comment ======================================= */
00733 /* =========================================================== */
00734 
00735 int main(int argc,char **argv)
00736 {
00737   TABLE cubetable, dbg_table;
00738   TIGERfile cube;
00739   SPECTRUM spec, dbg_spec;
00740 
00741   gsl_cheb_series *continuum, *wcalib;
00742 
00743   ArcGroup *refArcGroup,*courant;
00744 
00745   char **arglabel,**argval, **namecol_c;
00746   char *refname, *cubename, *dbg_tablename, *dbg_specname;
00747 
00748   double *errcal, *Weightcal, *fitLbda,*refLbda,*fitInt,*fitSig,*fitError;
00749   double *coefcont, *coefcal;
00750   double dbl_desc[4];
00751   double err_mean, err_std, err_min, err_max;
00752   double halfwidth, Wbound, sigmax,sigmin, lborne1, lborne2, delta_borne;
00753   double start_polcal, end_polcal;
00754   double error, error_norm, perc_nrej, vardouble;
00755   double siginf_cont=-3, sigsup_cont=0.5; /* Continuum rejection (orig: -3,1)*/
00756   double siginf_cal=-2,  sigsup_cal=2;
00757 
00758   float *x0, *y0, *flbda, valfloat;
00759 
00760   long nb_error=0, nflag1, nflag2, nflag3;
00761   long NbPixel, varlong;
00762   long ncoefcont, ncoefcal;
00763 
00764   int *fitFlag, *nolens, *id_c, *il;
00765   int degcont, degcal;
00766   int singlelens, i, j, k, k_min, nbligne, nlbda, jstart, class;
00767   int col_no, ntotlines, nrej, npts_cal, ndiff, flag, status;
00768   int id_no,id_lref,id_x,id_i,id_sig, id_dx, id_fit, id_ifail;
00769   int id_deg, id_xmin, id_xmax, id_nrej, id_npts, id_x0, id_y0;
00770   int id_error, id_flag, id_weightcal;
00771   int nbgoodcal, niter, itermax_cont=10, itermax_cal=20;
00772 
00773   set_purpose("Wavelength calibration (computation)");
00774   set_arglist("-in none -reftable LbdaRef -N|halfwidth 6 "
00775               "-sigmax 2.0 -sigmin 0.5 "
00776               "-lens|lensno 0 -dbgtable dbg_wcalib -dbgspec dbg_wcalfit "
00777               "-degcont 3 -degcal 3 -W|edge 1");
00778 
00779   init_snifs("$Name:  $","$Revision: 1.27 $");
00780   init_session(argv,argc,&arglabel,&argval);
00781 
00782   if (DEBUG) {
00783     print_msg("$Id: wcalib.c,v 1.27 2005/10/25 00:08:02 ycopin Exp $");
00784     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00785   }
00786 
00787   cubename = argval[0];                      /* Input datacube */
00788   refname  = argval[1];                      /* Ref. wavelength table */
00789   get_argval(2,"%lf",&halfwidth);            /* Half-width line window */
00790   get_argval(3,"%lf",&sigmax);
00791   get_argval(4,"%lf",&sigmin);
00792 
00793   get_argval(5,"%d",&singlelens);            /* Single lens calibration */
00794   if (DEBUG) {                               /* Debug output */
00795     dbg_tablename = argval[6];
00796     if (singlelens) dbg_specname  = argval[7];
00797   }
00798 
00799   get_argval(8,"%d",&degcont);
00800   get_argval(9,"%d",&degcal);
00801   if (degcal+1 > 9) {    /* There's not enough room in column name */ 
00802     print_error("Degree of calibration is too high! (<9)");
00803     exit_session(ERR_BAD_PARAM);
00804   }
00805 
00806   get_argval(10,"%lf",&Wbound); /* minimum distance (in sigma) of a line to
00807                                    the edges of the fitting domain */
00808 
00809   /*======================================================================*/
00810   /* Ouverture du cube TIGER */
00811   /*======================================================================*/
00812 
00813   print_msg("Opening datacube %s", cubename);
00814   if (open_tiger_frame(&cube,cubename,"I")<0){
00815     print_error("Cannot open datacube %s",cubename);
00816     exit_session(ERR_OPEN);
00817   }
00818 
00819   switch ((class = read_file_class(&cube))) {
00820   case SUPER_CLASS:
00821   case RAW_CAL_CUBE: break;
00822   default:
00823     print_error("Wrong fclass (%d) for %s (should be "
00824                 "Raw Calibration Datacube)",class,cube.name);
00825 
00826     print_msg("DEBUG: Pb with datacube FCLASS, keep going...");
00827     //      exit_session(ERR_BAD_TYPE);
00828   }
00829   nolens=(int *)malloc(cube.nbspec*sizeof(int)  );
00830   il =   (int *)malloc(cube.nbspec*sizeof(int)  );
00831   x0 = (float *)malloc(cube.nbspec*sizeof(float));
00832   y0 = (float *)malloc(cube.nbspec*sizeof(float));
00833 
00834   /* get_lenses_no(&cube,nolens); */
00835   get_lenses_coord(&cube,"X0","Y0",nolens,x0,y0,il);
00836 
00837   if (!singlelens) {
00838     /*======================================================================*/
00839     /* Ouverture de la table associee au cube */
00840     /*======================================================================*/
00841 
00842     if (open_table(&cubetable,cube.table_name,"IO")<0){
00843       close_tiger_frame(&cube);
00844       print_error("Cannot open datacube's table",cube.table_name);
00845       exit_session(ERR_OPEN);
00846     }
00847     switch ((class = read_file_class(&cubetable))) {
00848     case SUPER_CLASS:
00849     case TBL_FILE: break;
00850     case TBL_FILE_CAL:
00851       print_warning("Computating new calibration for table %s", cubetable.name);
00852       break;
00853     default:
00854       print_error("Wrong fclass (%d) for %s (should be a "
00855                   "[Calibrated] table)",class,cubetable.name);
00856       exit_session(ERR_BAD_TYPE);
00857     }
00858     col_no = get_col_ref(&cubetable, LAB_COL_NO);
00859   } 
00860   else print_warning("DEBUG: table %s untouched", cube.table_name);
00861 
00862   /*======================================================================*/
00863   /* retrieve filter bounds in the datacube descriptor */
00864   /*======================================================================*/
00865 
00866   if (RD_desc(&cube,LRANGE4,DOUBLE,4,dbl_desc) < 0) {
00867     print_error("Cannot read descriptor %s from datacube %s",LRANGE4,cube.name);
00868 
00869     print_warning("DEPRECATED:\n"
00870                   "Previous versions of extract_spec (bef. 10/11/04, v1.31)\n"
00871                   "were writing an invalidly-long keyword "
00872                   "LBDARANGn (n=1...4).\n"
00873                   "See if you can find it back and rename it to LRANGEn.\n");
00874     exit_session(ERR_NODESC);
00875   }
00876   else {
00877     lborne1 = dbl_desc[0];
00878     lborne2 = dbl_desc[1];
00879     delta_borne = lborne2 - lborne1;
00880     print_msg("Spectral domain: %.1f-%.1f AA",lborne1,lborne2);
00881   }
00882 
00883   /*======================================================================*/
00884   /* convert sigma bounds to AA */
00885   /*======================================================================*/
00886   sigmin *= cube.step;
00887   sigmax *= cube.step;
00888   if (DEBUG)
00889     print_msg("Spectral resolution sigma range: %.2f-%.2f AA",sigmin,sigmax);
00890 
00891   /*======================================================================*/
00892   /* Ouverture de la table des raies d'emission de reference */
00893   /*======================================================================*/
00894 
00895   print_msg("Reading wavelength reference table %s", refname);
00896   if ((nlbda = read_lbdaref_table(refname,"LAMBDA","CALIB",
00897                                   &flbda,lborne1,lborne2)) < 2) {
00898     print_error("Not enough lines read from %s",refname);
00899     exit_session(ERR_NODATA);
00900   }
00901 
00902   /* Convert to DOUBLE */
00903   refLbda = (double *)malloc(nlbda*sizeof(double));
00904   for (i=0; i<nlbda; i++) refLbda[i] = flbda[i];
00905   free(flbda);
00906 
00907   /*======================================================================*/
00908   /* Regroup peaks according to proximity */
00909   /*======================================================================*/
00910    
00911   refArcGroup = findArcGroups(nlbda, refLbda, halfwidth, cube.step);
00912 
00913   /* Degree -> Nb of coeff. */
00914   ncoefcont = degcont + 1;
00915   ncoefcal  = degcal  + 1;
00916 
00917   /*======================================================================*/
00918   /* allocation memoire */
00919   /*======================================================================*/
00920   fitLbda   = (double *)malloc(nlbda * sizeof(double));
00921   fitError  = (double *)malloc(nlbda * sizeof(double));
00922   Weightcal = (double *)calloc(nlbda,  sizeof(double));
00923   fitInt    = (double *)malloc(nlbda * sizeof(double));
00924   fitSig    = (double *)malloc(nlbda * sizeof(double));
00925   fitFlag   =    (int *)malloc(nlbda * sizeof(int));
00926   coefcont  = (double *)malloc(ncoefcont * sizeof(double));
00927   coefcal   = (double *)malloc(ncoefcal  * sizeof(double));
00928   errcal    = (double *)malloc(cube.nbspec * sizeof(double));
00929 
00930   continuum = gsl_cheb_alloc(degcont);
00931   wcalib    = gsl_cheb_alloc(degcal);
00932 
00933   /*======================================================================*/
00934   /* creation de la table de debug et des colonnes */
00935   /*======================================================================*/
00936 
00937   if (DEBUG) {
00938     print_msg("Creating debug table %s", dbg_tablename);
00939     if (create_table(&dbg_table,dbg_tablename,-1,7,'Q',"Wcalib table")<0) {
00940       print_error("Cannot create table %s",dbg_tablename);
00941       exit_session(ERR_CREAT);
00942     }
00943     else {
00944       id_no  = create_col(&dbg_table,LAB_COL_NO,LONG, 'O',"I4",  "Lens Number");
00945       id_x0  = create_col(&dbg_table,"X0",      FLOAT,'O',"F9.6","px");
00946       id_y0  = create_col(&dbg_table,"Y0",      FLOAT,'O',"F9.6","px");
00947       id_lref= create_col(&dbg_table,"REF",     FLOAT,'O',"F9.6","AA");
00948       id_x   = create_col(&dbg_table,"X",       FLOAT,'O',"F9.6","AA");
00949       id_i   = create_col(&dbg_table,"I",       FLOAT,'O',"F9.6","e-");
00950       id_sig = create_col(&dbg_table,"SIG",     FLOAT,'O',"F9.6","AA");
00951       id_dx  = create_col(&dbg_table,"DX",      FLOAT,'O',"F9.6","AA");
00952       id_fit = create_col(&dbg_table,"XFIT",    FLOAT,'O',"F9.6","AA");
00953       id_ifail     = create_col(&dbg_table,"IFAIL",LONG, 'O',"I2","");
00954       id_weightcal = create_col(&dbg_table,"WCAL", FLOAT,'O',"F3.1","number");
00955     }
00956   }
00957 
00958   k_min = -1;
00959 
00960   if (singlelens) {
00961 
00962     /*======================================================================*/
00963     /* Une seule lentille */
00964     /*======================================================================*/
00965 
00966     j=0;
00967     jstart=0;
00968     nbgoodcal = 0;
00969     
00970     print_msg("Starting computation for debug lens %d",singlelens);
00971     
00972     get_tiger_spec(&cube,&spec,NULL,singlelens);
00973 
00974     /*======================================================================*/
00975     /* Ajustement du continu sous les raies: pour cela on met
00976        une rejection importante positive et un nombre d'iterations
00977        de 10 */
00978     /*======================================================================*/
00979 
00980     fit_poly_rej_nag(&spec, ncoefcont, siginf_cont, sigsup_cont, 
00981                      itermax_cont, coefcont, &niter, &error, &nrej,
00982                      spec.start, spec.end, &status);
00983 
00984     /* Convert continuum to GSL Chebyshev series */
00985     setGslCheb(continuum, coefcont, ncoefcont, spec.start, spec.end);
00986 
00987     courant = refArcGroup; 
00988     k = 0; 
00989     ntotlines = 0;
00990     glLbdaStep = spec.step;
00991 
00992     print_msg("Scanning the calibration lines set");
00993     while (courant != NULL) {
00994 
00995       /*====================================================================*/
00996       /* attention: fillSpectrumArray inclut la soustraction du continu
00997          qui a ete ajuste ci-dessus */
00998       /*====================================================================*/
00999 
01000       NbPixel = fillSpectrumArray(&spec, courant->lbda_inf, courant->lbda_sup, 
01001                                   continuum);
01002 
01003       if (NbPixel < 5) {                     /* Skip small line groups */
01004         k += courant->nlines;
01005         courant = courant->next;
01006         continue;
01007       }
01008       glnlines   = courant->nlines;
01009       ntotlines += glnlines;
01010 
01011       /*====================================================================*/
01012       /* ajustement des gaussiennes */
01013       /*====================================================================*/
01014 
01015       if (k_min < 0) k_min = k;
01016 
01017       //fitLineSet(NbPixel, Wbound, &(refLbda[k]), &(fitLbda[k]), &(fitSig[k]),
01018       //           &(fitInt[k]), &(fitFlag[k]), sigmin, sigmax);
01019       fitLineSet_constrained(NbPixel, Wbound, &(refLbda[k]), &(fitLbda[k]), 
01020                              &(fitSig[k]), &(fitInt[k]), &(fitFlag[k]), 
01021                              sigmin, sigmax);
01022 
01023       for (i=0; i<glnlines; i++) {
01024         fitError[k+i] = fitLbda[k+i] - refLbda[k+i];
01025 
01026         /*==================================================================*/
01027         /* si l'ajustement de cette raie n'a pas marche on elimine */
01028         /*==================================================================*/
01029 
01030         if (fitFlag[k+i] <= 1)
01031           Weightcal[k+i] = 1.;
01032         else 
01033           Weightcal[k+i] = 1.e-99;
01034 
01035         if (DEBUG) {
01036           print_msg("Line %2d (%.2f, %.2f-%.2f): dl=%.2f (S=%.2f,I=%g,flag=%d)",
01037                     j+1, refLbda[k+i], 
01038                     glLbdaStart, glLbdaStart+(NbPixel-1)*glLbdaStep, 
01039                     fitError[k+i], fitSig[k+i], fitInt[k+i], fitFlag[k+i]);
01040 
01041           WR_tbl(&dbg_table,j,id_no,&singlelens);
01042           WR_tbl(&dbg_table,j,id_lref,(float[]){refLbda[k+i]});
01043           WR_tbl(&dbg_table,j,id_x,   (float[]){fitLbda[k+i]});
01044           WR_tbl(&dbg_table,j,id_i,   (float[]){fitInt[k+i]});
01045           WR_tbl(&dbg_table,j,id_sig, (float[]){fitSig[k+i]});
01046           WR_tbl(&dbg_table,j,id_dx,  (float[]){fitError[k+i]});
01047           WR_tbl(&dbg_table,j,id_ifail,&(fitFlag[k+i]));
01048           j++;
01049         }
01050       }
01051 
01052       k += courant->nlines;
01053       courant=courant->next;
01054     }
01055 
01056     if (DEBUG) {
01057       print_msg("Saving fit result as debug spectrum %s",dbg_specname);
01058       disable_erase_flag();
01059       create_spec(&dbg_spec,dbg_specname,spec.npts,spec.start,
01060                   spec.step,DOUBLE,"essai","AA");
01061 
01062       saveGauss(&dbg_spec,fitLbda+k_min,fitSig+k_min,fitInt+k_min, 
01063                 ntotlines, continuum);
01064       write_file_class(&dbg_spec, DEBUG_FILE);
01065       close_spec(&dbg_spec);
01066     }
01067 
01068     /*======================================================================*/
01069     /* On determine les bornes a passer au polynome */
01070     /* ceci au cas ou le centre d'une raie ajustee tomberait hors spectre */
01071     /*======================================================================*/
01072     start_polcal = MIN(lborne1, refLbda[0]);
01073     end_polcal = MAX(lborne2, refLbda[nlbda-1]);
01074 
01075     /*======================================================================*/
01076     /* Ajustement du polynome de calibration sur les raies ajustees */
01077     /*======================================================================*/
01078     if (ntotlines < ncoefcal) {
01079       flag = 3;
01080       print_error("Spectrum %d cannot be calibrated: not enough lines to fit",
01081                   singlelens);
01082       exit_session(UNKNOWN);
01083     }
01084 
01085     fit_poly_rej_nag_tab(refLbda+k_min,fitError+k_min,Weightcal+k_min,
01086                          ntotlines, 
01087                          ncoefcal,siginf_cal,sigsup_cal,itermax_cal,coefcal,
01088                          &niter, &error, &nrej, start_polcal, end_polcal, 
01089                          &status);
01090 
01091     /* Convert calibration to GSL Chebyshev series */
01092     setGslCheb(wcalib, coefcal, ncoefcal, start_polcal, end_polcal);
01093 
01094     if (DEBUG) {
01095       for (i=0; i<ntotlines; i++) {
01096         if (Weightcal[i+k_min] < 1.e-20) valfloat = 0.;
01097         else                             valfloat = Weightcal[i+k_min];
01098         WR_tbl(&dbg_table,i,id_weightcal,&valfloat);
01099       }
01100     }
01101 
01102     /*======================================================================*/
01103     /*  Calcul du flag 0, 1 ou 2 du spectre courant */
01104     /*======================================================================*/
01105 
01106     switch(status) {
01107     case -1:
01108       print_error("Spectrum %d not calibrated: bounds error",singlelens);
01109       exit_session(UNKNOWN);
01110     case -2:
01111       flag = 3;
01112       print_error("Spectrum %d not calibrated: lack of lines to fit",
01113                   singlelens);
01114       exit_session(UNKNOWN);
01115     default:
01116       nbgoodcal = nbgoodcal + 1;
01117       if (DEBUG) 
01118         for (i=0; i<ntotlines; i++) {
01119           valfloat = gsl_cheb_eval(wcalib, refLbda[i+k_min]);
01120           WR_tbl(&dbg_table,i,id_fit,&valfloat);
01121         }
01122       break;
01123 
01124     }
01125 
01126     error_norm = error / glLbdaStep;         /* In pixels */
01127     perc_nrej  = (double)nrej / ntotlines;
01128     npts_cal   = ntotlines - nrej;
01129     ndiff      = npts_cal - ncoefcal;
01130     
01131     if (ndiff < 0) {
01132       flag = 3;
01133       print_error("Spectrum %d not calibrated: not enough lines to fit",
01134                   singlelens);
01135       exit_session(UNKNOWN);
01136     }
01137 
01138     flag = 0;
01139     if ((error_norm > CRITER_SIG_2) || 
01140         (perc_nrej > CRITER_NREJ_2) || 
01141         (ndiff < CRITER_DIFF_2)) {
01142       flag = 2;
01143     }
01144     else if ((error_norm > CRITER_SIG_1) || 
01145              (perc_nrej > CRITER_NREJ_1) || 
01146              (ndiff < CRITER_DIFF_1)) {
01147       flag = 1;
01148     }
01149 
01150     /*======================================================================*/
01151     /* Ecriture des resultats dans la table */
01152     /*======================================================================*/
01153     print_msg("%d lines in total: %d used, %d rejected, residual: %f",
01154               ntotlines, npts_cal, nrej, error_norm);
01155     switch (flag) {
01156     case 0: print_msg("Fit quality: EXCELLENT"); break;
01157     case 1: print_msg("Fit quality: GOOD");      break;
01158     case 2: print_msg("Fit quality: FAIR");      break;
01159     case 3: print_msg("Fit quality: POOR");      break;
01160     }
01161 
01162     free_spec_mem(&spec);
01163 
01164   }
01165   else {
01166 
01167     /*======================================================================*/
01168     /* Toutes les lentilles */
01169     /*======================================================================*/
01170 
01171     nbligne = 0;
01172 
01173     /*======================================================================*/
01174     /* Creation des colonnes de sortie pour la table associee au cube */
01175     /*======================================================================*/
01176 
01177     id_xmin = create_col(&cubetable,"XMIN",DOUBLE,'R',"F9.6","AA");
01178     id_xmax = create_col(&cubetable,"XMAX",DOUBLE,'R',"F9.6","AA");
01179     id_deg  = create_col(&cubetable,"DEGREE",INT,'R',"I3","number");
01180 
01181     id_c      =   (int *)malloc(ncoefcal*sizeof(int));
01182     namecol_c = (char **)malloc(ncoefcal*sizeof(char *));
01183     for (i=0; i<ncoefcal; i++) {
01184       namecol_c[i] = (char *)malloc(3 * sizeof(char)); /* ncoefcal up to 9 */
01185       sprintf(&(namecol_c[i][0]),"C%d",i);
01186       id_c[i] = create_col(&cubetable,namecol_c[i],DOUBLE,'R',"F9.6","number");
01187     }
01188     id_nrej = create_col(&cubetable,"NREJ", INT,'R',"I4","number");
01189     id_error= create_col(&cubetable,"ERROR",DOUBLE,'R',"F9.6","px");
01190     id_npts = create_col(&cubetable,"NFIT",INT,'R',"I4","number");
01191     id_flag = create_col(&cubetable,"FLAG",INT,'R',"I4","number");
01192 
01193     print_msg("Starting computation for all lenses");
01194 
01195     reset_print_progress();
01196     for (i=0; i<cube.nbspec; i++) {          
01197       get_tiger_spec(&cube,&spec,NULL,nolens[i]);
01198 
01199       print_progress("Calibration in progress",(i+1.)/cube.nbspec*100.,1.0);
01200 
01201       /*====================================================================*/
01202       /* Ajustement du continu sous les raies: pour cela on met
01203          une rejection importante positive et un nombre d'iterations
01204          de 10 */
01205       /*====================================================================*/
01206       fit_poly_rej_nag(&spec, ncoefcont, siginf_cont, sigsup_cont, 
01207                        itermax_cont, coefcont, &niter, &error, &nrej, 
01208                        spec.start, spec.end, &status);
01209 
01210       /* Convert continuum to GSL Chebyshev series */
01211       setGslCheb(continuum, coefcont, ncoefcont, spec.start, spec.end);
01212 
01213       courant = refArcGroup; 
01214       k = 0; 
01215       ntotlines = 0;
01216       glLbdaStep = spec.step;
01217 
01218       while (courant != NULL) {
01219         /*==================================================================*/
01220         /* ATTENTION: fillSpectrumArray inclus la soustraction du continu
01221            qui a ete ajuste ci-dessus */
01222         /*==================================================================*/
01223         NbPixel = fillSpectrumArray(&spec, courant->lbda_inf,courant->lbda_sup,
01224                                     continuum);
01225 
01226         if (NbPixel < 5) {
01227           k += courant->nlines;
01228           courant = courant->next;
01229           continue;
01230         }
01231         glnlines   = courant->nlines;
01232         ntotlines += glnlines;
01233 
01234         /*==================================================================*/
01235         /* ajustement des gaussiennes */
01236         /*==================================================================*/
01237           
01238         if (k_min < 0) k_min = k;
01239 
01240         //fitLineSet(NbPixel, Wbound, &(refLbda[k]), &(fitLbda[k]), &(fitSig[k]),
01241         //           &(fitInt[k]), &(fitFlag[k]), sigmin, sigmax);
01242         fitLineSet_constrained(NbPixel, Wbound, &(refLbda[k]), &(fitLbda[k]), 
01243                                &(fitSig[k]), &(fitInt[k]), &(fitFlag[k]), 
01244                                sigmin, sigmax);
01245 
01246         for (j=0; j<glnlines; j++) { 
01247           fitError[k+j] = fitLbda[k+j] - refLbda[k+j];
01248 
01249           /*================================================================*/
01250           /* si l'ajustement de cette raie n'a pas marche on elimine */
01251           /*================================================================*/
01252           if (fitFlag[k+j] <= 1)
01253             Weightcal[k+j] = 1.;
01254           else
01255             Weightcal[k+j] = 1.e-99;
01256         }
01257 
01258         if (DEBUG) {
01259           for (j=0; j<glnlines; j++) {
01260             WR_tbl(&dbg_table,nbligne,id_no,  &(nolens[i]));
01261             WR_tbl(&dbg_table,nbligne,id_x0,  &(x0[i]));
01262             WR_tbl(&dbg_table,nbligne,id_y0,  &(y0[i]));
01263             WR_tbl(&dbg_table,nbligne,id_lref,&(refLbda[k+j]));
01264             WR_tbl(&dbg_table,nbligne,id_x,   &(fitLbda[k+j]));
01265             WR_tbl(&dbg_table,nbligne,id_i,   &(fitInt[k+j]));
01266             WR_tbl(&dbg_table,nbligne,id_sig, &(fitSig[k+j]));
01267             WR_tbl(&dbg_table,nbligne,id_dx,  &(fitError[k+j]));
01268             WR_tbl(&dbg_table,nbligne,id_ifail,&(fitFlag[k+i]));
01269             nbligne++;
01270           }
01271         }
01272         k += courant->nlines;
01273 
01274         courant=courant->next;
01275       }    
01276        
01277       /*====================================================================*/
01278       /* On determine les bornes a passer au polynome */
01279       /*====================================================================*/
01280       start_polcal = MIN(lborne1, refLbda[0]);
01281       end_polcal = MAX(lborne2, refLbda[nlbda-1]);
01282 
01283       /*====================================================================*/
01284       /* Ajustement du polynome sur les raies ajustees */
01285       /*====================================================================*/
01286       if (ntotlines < ncoefcal) {
01287         flag = 3;
01288         WR_tbl(&cubetable, il[i], id_flag, &flag);
01289         print_warning("Spectrum %d not calibrated: "
01290                       "not enough lines to fit (%d<%d)",
01291                       nolens[i],ntotlines,ncoefcal);
01292         continue;
01293       }
01294 
01295       fit_poly_rej_nag_tab(refLbda+k_min, fitError+k_min, Weightcal+k_min, 
01296                            ntotlines, ncoefcal, 
01297                            siginf_cal, sigsup_cal, itermax_cal, coefcal, 
01298                            &niter, &error, &nrej, start_polcal, end_polcal,
01299                            &status);
01300       
01301       /* Convert calibration to GSL Chebyshev series */
01302       setGslCheb(wcalib, coefcal, ncoefcal, start_polcal, end_polcal);
01303 
01304       /*====================================================================*/
01305       /*        Calcul du flag 0, 1 ou 2 du spectre courant */
01306       /*====================================================================*/
01307       switch(status) {
01308       case -1:
01309         flag = 3;
01310         WR_tbl(&cubetable, il[i], id_flag, &flag);
01311         print_warning("Spectrum %d not calibrated: error on bounds",
01312                       nolens[i]);
01313         continue;
01314       case -2:
01315         flag = 3;
01316         WR_tbl(&cubetable, il[i], id_flag, &flag);
01317         print_warning("Spectrum %d not calibrated: too many lines rejected",
01318                       nolens[i]);
01319         continue;
01320       default:
01321         nbgoodcal = nbgoodcal + 1;
01322         break;
01323 
01324       }
01325 
01326       error_norm = error / glLbdaStep;
01327       perc_nrej  = (double)nrej / ntotlines;
01328       npts_cal   = ntotlines - nrej;
01329       ndiff      = npts_cal - ncoefcal;
01330 
01331       if (ndiff < 0) {
01332         flag = 3;
01333         nflag3++;
01334         WR_tbl(&cubetable, il[i], id_flag, &flag);
01335         print_warning("Spectrum %d not calibrated: lack of lines to fit",
01336                       nolens[i]);
01337         continue;
01338       }
01339 
01340       flag = 0;
01341       if ((error_norm > CRITER_SIG_2) ||
01342           (perc_nrej  > CRITER_NREJ_2) ||
01343           (ndiff      < CRITER_DIFF_2)) {
01344         flag = 2;
01345         nflag2++;
01346       }
01347       else if ((error_norm > CRITER_SIG_1) || 
01348                (perc_nrej  > CRITER_NREJ_1) ||
01349                (ndiff      < CRITER_DIFF_1)) {
01350         flag = 1;
01351         nflag1++;
01352       }
01353       errcal[nb_error++] = error_norm;
01354 
01355       /*====================================================================*/
01356       /* Ecriture des resultats dans la table */
01357       /*====================================================================*/
01358       WR_tbl(&cubetable, il[i], id_xmin, &start_polcal);
01359       WR_tbl(&cubetable, il[i], id_xmax, &end_polcal);
01360       WR_tbl(&cubetable, il[i], id_nrej, &nrej);
01361       WR_tbl(&cubetable, il[i], id_npts, &npts_cal);
01362       WR_tbl(&cubetable, il[i], id_error,&error_norm);
01363       WR_tbl(&cubetable, il[i], id_flag, &flag);
01364       varlong = degcal;
01365       WR_tbl(&cubetable, il[i], id_deg, &varlong);
01366       for (j=0; j<ncoefcal; j++) {
01367         vardouble = (wcalib->c)[j];
01368         WR_tbl(&cubetable, il[i], id_c[j], &vardouble);
01369       }
01370       free_spec_mem(&spec);
01371     }
01372 
01373     /* Statistics on error */
01374     err_mean = gsl_stats_mean(errcal,1,nb_error);
01375     err_std  = gsl_stats_sd_m(errcal,1,nb_error,err_mean);
01376     err_std  = gsl_isnan(err_std)? 0:err_std; /* Discard NaN */
01377     gsl_stats_minmax(&err_min,&err_max,errcal,1,nb_error);
01378 
01379     print_msg("Calibration Mean Error: %5.3f +- %5.3f px [%5.3f - %5.3f]",
01380               err_mean, err_std, err_min, err_max);
01381     /* Use C99 compounds literals (supported from GCC 3.1) */
01382     WR_desc(&cubetable, WLMEANERR,FLOAT, 1, (float[]){err_mean});
01383     WR_desc(&cubetable, WLSTDERR, FLOAT, 1, (float[]){err_std});
01384     WR_desc(&cubetable, WLMINERR, FLOAT, 1, (float[]){err_min});
01385     WR_desc(&cubetable, WLMAXERR, FLOAT, 1, (float[]){err_max});
01386     WR_desc(&cubetable, WLMAXDEGRE, INT, 1, &ncoefcal);
01387 
01388     /* write type of table (wavelength calibrated table) */
01389     write_file_class(&cubetable, TBL_FILE_CAL);
01390   }
01391 
01392   /* delArcGroups(refArcGroup); */
01393   free(refLbda);
01394 
01395   close_tiger_frame(&cube);
01396   if (!singlelens) close_table(&cubetable);
01397   if (DEBUG) {
01398     write_file_class(&dbg_table, DEBUG_FILE);
01399     close_table(&dbg_table);
01400   }
01401 
01402   /* Nombre total de spectres etalonnes avec succes, en mode non-DEBUG */
01403 
01404   if (!DEBUG) 
01405     print_msg("-> %d spectra successfully calibrated",nbgoodcal);
01406 
01407   /* Liberation de la memoire */
01408   free(nolens);
01409   free(fitLbda); free(fitError); free(Weightcal); 
01410   free(fitInt); free(fitSig); free(fitFlag);
01411   free(coefcont); free(coefcal);
01412 
01413   gsl_cheb_free(continuum);
01414 
01415   free(glArray);
01416 
01417   exit_session(OK);
01418   return(OK);
01419 }    
01420 
01421 /* 
01422 ;;; Local Variables: ***
01423 ;;; eval: (add-to-list 'c-font-lock-extra-types "ArcLine") ***
01424 ;;; eval: (add-to-list 'c-font-lock-extra-types "ArcGroup") ***
01425 ;;; End: ***
01426 */

Generated on Wed Oct 26 23:59:39 2005 for Snifs by doxygen 1.3.5