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 
00025 /* variable pour l'initialisation de la FWHM pour le fit gaussien */
00026 #define FWHMINIT (double)2.
00027 
00028 /* Criteres pour etablir la qualite du fit de calibration: flag */
00029 #define CRITER_SIG_1    (double)0.05
00030 #define CRITER_SIG_2    (double)0.10
00031 #define CRITER_NREJ_1   (double)0.2
00032 #define CRITER_NREJ_2   (double)0.5
00033 #define CRITER_DIFF_1   (int)3
00034 #define CRITER_DIFF_2   (int)2
00035 
00036 /* global variable needed for MultiGauss */
00037 
00038 double lbdastart;
00039 double lbdastep;
00040 int NbGauss;
00041 double *SPECTRE;
00042 
00045 typedef struct ArcLine
00046 {
00047   double          lambda;    
00048   struct ArcLine *next;      
00049 } ArcLine;
00050 
00053 typedef struct ArcGroup 
00054 {
00055   int              nb;                
00056   double           lbda_inf,lbda_sup; 
00057   ArcLine         *lines;             
00058   struct ArcGroup *next;              
00059 } ArcGroup;
00060 
00061 /* === Doxygen Comment ======================================= */
00068 /* =========================================================== */
00069 
00070 void ArcLine_Add(ArcLine **list,double lambda)
00071 {
00072   ArcLine *new,*courant,*precedent;
00073 
00074   /* Create the new line */
00075   new = (ArcLine *)malloc(sizeof(ArcLine));
00076   new->lambda = lambda;
00077   new->next   = NULL;
00078 
00079   if (*list == NULL) {
00080     /* Initialize a new list */
00081     *list = new;
00082   } 
00083   else {
00084     /* Update list: Add reference to new line in preceeding line */
00085     courant = precedent = *list;
00086     while (courant != NULL) { /* Loop over the lines */
00087       precedent = courant;
00088       courant   = courant->next;
00089     }
00090     precedent->next = new;
00091   }
00092 }
00093 
00094 /* === Doxygen Comment ======================================= */
00101 /* =========================================================== */
00102 
00103 ArcGroup *ArcGroup_New(ArcGroup **list,double lambda)
00104 {
00105   ArcGroup *new,*courant,*precedent;
00106 
00107   /* Create a brand new group with the single line */
00108   new = (ArcGroup *)malloc(sizeof(ArcGroup));
00109   new->nb       = 1;
00110   new->lbda_inf = 0;
00111   new->lbda_sup = 0;
00112   new->lines    = NULL;
00113   new->next     = NULL;
00114   ArcLine_Add(&(new->lines),lambda);
00115 
00116   if (*list == NULL) {
00117     /* Initialize a new list */
00118     *list = new;
00119   } 
00120   else {
00121     /* Update list: Add reference to new group in preceeding group */
00122     courant = precedent = *list;
00123     while (courant != NULL) { /* Loop over the groups */
00124       precedent = courant;
00125       courant   = courant->next;
00126     }
00127     precedent->next = new;
00128   }
00129 
00130   return(new);
00131 }
00132 
00133 /* === Doxygen Comment ======================================= */
00140 /* =========================================================== */
00141 
00142 void ArcGroup_Add(ArcGroup **list,double lambda)
00143 { 
00144   ArcLine_Add(&((*list)->lines),lambda);
00145   (*list)->nb++;
00146 }
00147 
00148 /* === Doxygen Comment ======================================= */
00154 /* =========================================================== */
00155 
00156 void ArcLine_Delete(ArcLine *list)
00157 {
00158   /* Recursive deletion */
00159   if (list->next != NULL) ArcLine_Delete(list->next);
00160 
00161   free(list);
00162 }
00163 
00164 /* === Doxygen Comment ======================================= */
00170 /* =========================================================== */
00171 
00172 void ArcGroup_Delete(ArcGroup *list)
00173 {
00174   if (list->next != NULL) {
00175     ArcGroup_Delete(list->next);
00176     ArcLine_Delete(list->lines);
00177     free(list);
00178   }
00179   else {
00180     ArcLine_Delete(list->lines);
00181     free(list);
00182   }
00183 }
00184 
00185 /* === Doxygen Comment ======================================= */
00196 /* =========================================================== */
00197 
00198 int Take_Ref(TABLE *table,double **ref,double borne1,double borne2)
00199 {
00200   int i,nb=0,j, imin, imax;
00201   double *tbl;
00202   float val=0.;
00203   int coll=get_col_ref(table,"LAMBDA");
00204   int colc=get_col_ref(table,"CALIB");
00205 
00206   for (i=0;(i<table->row) && (val < borne1);i++)
00207     RD_tbl(table,i,coll,&val);
00208 
00209   imin = i-1;
00210   for (;(i<table->row) && (val <= borne2);i++)
00211     RD_tbl(table,i,coll,&val);
00212 
00213   imax = i-1;
00214   nb = (imax -imin)+1;
00215   tbl = (double *)malloc(sizeof(double)*nb);
00216 
00217   for (j=0,i=imin;i<=imax;i++) {
00218     if (RD_tbl(table,i,colc,&val) != ERR_NODATA)
00219       if (val == 0) continue;
00220     if (RD_tbl(table,i,coll,&val) != ERR_NODATA)
00221       tbl[j++]=(double)val;
00222   }
00223 
00224   *ref = tbl;
00225   return(j);
00226 }
00227 
00228 /* === Doxygen Comment ======================================= */
00237 /* =========================================================== */
00238 
00239 void Eval_bounds(ArcGroup *groupe,double Nhw,double dlambda)
00240 {
00241   ArcLine *courant,*precedent;
00242 
00243   courant = precedent = groupe->lines;
00244   
00245   groupe->lbda_inf = courant->lambda - (Nhw*dlambda);
00246 
00247   while (courant!=NULL) {
00248     precedent=courant;
00249     courant=courant->next;
00250   }
00251   
00252   groupe->lbda_sup = precedent->lambda + (Nhw*dlambda);
00253 }
00254 
00255 /* === Doxygen Comment ======================================= */
00266 /* =========================================================== */
00267 
00268 ArcGroup *Join_Ref_Peak(int size,double *array,double Nhw,double dlambda)
00269 {
00270   int i;
00271   ArcGroup *groupe=NULL,*last=NULL;
00272 
00273   if (size==1) {
00274     last = ArcGroup_New(&groupe,array[0]);
00275   }
00276 
00277   for (i=1;i<size;i++) {
00278     if ((array[i]-array[i-1])<(2*Nhw*dlambda)) {
00279       if (last==NULL) {
00280         last = ArcGroup_New(&groupe,array[i-1]);
00281         ArcGroup_Add(&last,array[i]);
00282       }
00283       else {
00284         ArcGroup_Add(&last,array[i]);
00285       }
00286     }
00287     else {
00288       if (last==NULL) {
00289         last = ArcGroup_New(&groupe,array[i-1]);
00290         Eval_bounds(last,Nhw,dlambda);
00291         last = ArcGroup_New(&groupe,array[i]);
00292         Eval_bounds(last,Nhw,dlambda);
00293       }
00294       else {
00295         Eval_bounds(last,Nhw,dlambda); 
00296         last = ArcGroup_New(&groupe,array[i]);
00297       }
00298     }
00299   }
00300 
00301   Eval_bounds(last,Nhw,dlambda); 
00302   return(groupe);
00303 }
00304 
00305 /* === Doxygen Comment ======================================= */
00318 /* =========================================================== */
00319 
00320 int Get_Spectrum_Data(SPECTRUM *spectre,double inf, double sup,
00321                       double *coef, long ncoef, double xmin, double xmax)
00322 {
00323   int ninf,nsup,nb,i,j, status;
00324   double *array, lbda;
00325 
00326   ninf = MAX(0,pixel_spec(spectre,(float)inf));
00327   nsup = MIN(spectre->npts-1,pixel_spec(spectre,(float)sup));
00328 
00329   lbdastart = (double)coord_spec(spectre, ninf);
00330 
00331   nb = (nsup - ninf)+1;
00332 
00333   if (nb < 5) /* not between bounds */
00334     return(nb);
00335 
00336   array = (double *)malloc(sizeof(double)*nb);
00337 
00338   lbda = lbdastart;
00339 
00340   for (i=ninf,j=0;i<=nsup;i++,j++) {
00341     array[j]=(double)RD_spec(spectre,i);
00342 
00343     if ((lbda >= xmin) && (lbda <= xmax)) {
00344       array[j] -= (double)val_poly_nag(lbda, coef, ncoef, xmin, xmax, &status);
00345     }
00346 
00347     lbda += spectre->step;
00348   }
00349 
00350   SPECTRE = array; 
00351   return(nb);
00352 }
00353 
00354 /* === Doxygen Comment ======================================= */
00366 /* =========================================================== */
00367 
00368 void MultiGauss(double *par, double *fj, int *ldfj, int *igo, 
00369                 int *iopt, double *ropt)
00370 {
00371   int i, npts=*ldfj-1,k,l;
00372   double lbda,fsum,dx,dx2,expo,f;
00373 
00374   lbda = lbdastart;
00375   for (i=0; i<npts; i++) {
00376     fsum = 0;
00377     for (k=0, l=0; k<NbGauss; k++, l+=3) {
00378       dx = (par[l+1]-lbda)/par[l+2];
00379       dx2 = dx*dx;
00380       expo = exp(-0.5*dx2);
00381       f = par[l]*expo;
00382       fsum += f;
00383       fj[l**ldfj+i] = expo; /* d f/d Ik */
00384       fj[(l+1)**ldfj+i] = -f*dx/par[l+2]; /* d f/d Xk */
00385       fj[(l+2)**ldfj+i] =  f*dx2/par[l+2]; /* d f/d Sigmak */
00386     }
00387     fj[(NbGauss*3)**ldfj+i] = fsum - SPECTRE[i]; /* valeur de f - measured */
00388     lbda += lbdastep;
00389   }
00390 }
00391 
00392 /* === Doxygen Comment ======================================= */
00407 /* =========================================================== */
00408 
00409 void MultiGauss_nag(long *MODE, long *npts, long *nvar, long *ldfj,
00410                     double *par, double *f, double *fj, long *nstate,
00411                     long *iuser, double *ruser)
00412 {
00413   long i,k,l;
00414   double lbda,fsum,dx,dx2,expo,func;
00415 
00416   lbda = lbdastart;
00417   for (i=0; i<*npts; i++) {
00418     fsum = 0.;
00419     for (k=0, l=0; k<NbGauss; k++, l+=3) {
00420       dx = (par[l+1]-lbda)/par[l+2];
00421       dx2 = dx*dx;
00422       expo = exp(-0.5*dx2);
00423       func = par[l]*expo;
00424       fsum += func;
00425       if ((*MODE == 1) || (*MODE == 2)) {
00426         fj[l**ldfj+i] = expo; /* d f/d Ik */
00427         fj[(l+1)**ldfj+i] = -func*dx/par[l+2]; /* d f/d Xk */
00428         fj[(l+2)**ldfj+i] =  func*dx2/par[l+2]; /* d f/d Sigmak */
00429       }
00430     }
00431     if ((*MODE == 0) || (*MODE == 2)) {
00432       f[i] = fsum - SPECTRE[i]; /* valeur de f - measured */
00433     }
00434     lbda += lbdastep;
00435   }
00436 }
00437 
00438 /* === Doxygen Comment ======================================= */
00453 /* =========================================================== */
00454 
00455 void SaveGauss(SPECTRUM *spectre, double *X, double *SIGMA, double *I,
00456                int ng, double *coef, long ncoef, double xmin, double xmax)
00457 {
00458   int i,k, status;
00459   double lbda,dx,dx2,expo,fsum,f;
00460 
00461   for (i=0; i<spectre->npts; i++) {
00462     lbda = coord_spec(spectre, i);
00463     fsum = 0;
00464     for (k=0; k<ng; k++) {
00465       dx = (X[k]-lbda)/SIGMA[k];
00466       dx2 = dx*dx;
00467       expo = exp(-0.5*dx2);
00468       f = I[k]*expo;
00469       fsum += f;
00470     }
00471     if ((lbda >= xmin) && (lbda <= xmax)) {
00472       fsum += (double)val_poly_nag(lbda, coef, ncoef, xmin, xmax, &status);
00473     }
00474     WR_spec(spectre,i,fsum);
00475   }
00476 }
00477 
00478 /* === Doxygen Comment ======================================= */
00488 /* =========================================================== */
00489 
00490 void SaveCalib(SPECTRUM *spectre, double *coefcal, 
00491                long ncoefcal, double start, double end)
00492 {
00493   int i, status;
00494   double lbda,f;
00495 
00496   for (i=0; i<spectre->npts; i++) {
00497     lbda = coord_spec(spectre, i);
00498     f = (double)val_poly_nag(lbda, coefcal, ncoefcal, start, end, &status);
00499     WR_spec(spectre,i,f);
00500   }
00501 }
00502 
00503 /* === Doxygen Comment ======================================= */
00517 /* =========================================================== */
00518 
00519 void Detection_peak(long NbPixel, double Wbound, double *X,
00520                     double *LAMBDA,double *SIGMA,
00521                     double *I,int *FLAG,double siginf,double sigsup)
00522 { 
00523   long nbvar;
00524   int i, j, k;
00525   double *par, *bl, *bu, residual, lbdaend;
00526   int ifail;
00527 
00528   nbvar = (long)(3*NbGauss);
00529   bl = (double *)malloc(nbvar*sizeof(double)); 
00530   bu = (double *)malloc(nbvar*sizeof(double)); 
00531   par = (double *)malloc(nbvar*sizeof(double));
00532 
00533   lbdaend = lbdastart + (NbPixel-1) * lbdastep;
00534 
00535   for (i=0, k=0; i<NbGauss; i++, k+=3) {
00536     j = (int)((LAMBDA[i]-lbdastart)/lbdastep + 0.5);   
00537     par[k] = ABS(SPECTRE[j]); /* intensite intiale */
00538     bl[k] = 0.; bu[k] = MAXDOUBLE;
00539     par[k+1] = LAMBDA[i]; /* long onde initiale */
00540     bl[k+1] = lbdastart; bu[k+1] = lbdaend;
00541     par[k+2] = FWHMINIT/2.355; /* sigma initial */
00542     bl[k+2] = siginf; bu[k+2] = sigsup;
00543   }
00544   
00545   ifail = lsq_fit_with_bounds_nonag(nbvar,NbPixel,par,bl,bu,
00546                                     MultiGauss,&residual,0);
00547   /*
00548     ifail = lsq_fit_with_bounds_nag(nbvar,NbPixel,par,bl,bu,
00549     MultiGauss_nag,&residual);
00550   */
00551 
00552   for (i=0, k=0; i<NbGauss; i++, k+=3) { 
00553     I[i] = par[k]; /* intensite finale */
00554     X[i] = par[k+1]; /* long onde */
00555     SIGMA[i] = par[k+2]; /* sigma  */
00556 
00557     /* Parfois, I=0 => erreur dans l'ajustement */
00558     if (I[i] == 0.)
00559       FLAG[i] = 98;
00560     /* test pour savoir si la raie est vraiment en bord de champ */
00561     else if ((Wbound > 0) && 
00562              ( ((X[i] - Wbound * SIGMA[i]) < lbdastart) || 
00563                ((X[i] + Wbound * SIGMA[i]) > lbdaend) )){
00564       FLAG[i] = 99;
00565     }
00566     else {
00567       FLAG[i] = ifail;
00568     }
00569   }
00570   free(bl); free(bu); free(par);
00571 }
00572 
00573 /* === Doxygen Comment ======================================= */
00593 /* =========================================================== */
00594 
00595 int main(int argc,char **argv)
00596 {
00597   TABLE tassoc,reftable,table;
00598   TIGERfile cube;
00599   SPECTRUM spectre, resultat;
00600 
00601   char **arglabel,**argval;
00602   char cubename[lg_name+1],*refname,*dbgtable,*dbgspec, *inname;
00603   char **namecol_c;
00604   double *errcal;
00605   double err_mean, err_std, err_min, err_max;
00606   float rdesc, *x0, *y0, valfloat;
00607   long nb_error=0, nflag1, nflag2, nflag3;
00608 
00609   double Nhw, Wbound, sigmax,sigmin, lborne1, lborne2, delta_borne, *LRef;
00610   int lensno, i, k, k_min, j, nbligne, LRefSize, *nolens, jstart;
00611   int class;
00612   long NbPixel;
00613 
00614   ArcGroup *RefPeak,*courant;
00615 
00616   double vardouble, *X,*LAMBDA,*I,*SIGMA, *dX;
00617   double *Weightcal, *dXcal, *LRefcal, dbl_desc[4];
00618   long varlong;
00619   int *FLAG, NbTotGauss, nrej, npts_cal, ndiff, flag, status;
00620   int col_no, id_no,id_lref,id_x,id_i,id_sig, id_dx, id_fit, id_ifail;
00621   int id_deg, id_xmin, id_xmax, *id_c, id_nrej, id_npts, id_x0, id_y0;
00622   int val_status, *il, id_error, id_flag, id_weightcal, CAT;
00623   int nbgoodcal;
00624 
00625   double *coefcont, start_polcont, end_polcont;
00626   double start_polcal, end_polcal, *coefcal;
00627   double error, error_norm, perc_nrej;
00628   double siginf_cont=-3, sigsup_cont=0.5; /* Continuum rejection (orig: -3,1)*/
00629   double siginf_cal=-2, sigsup_cal=2;
00630   int niter, itermax_cont=10, itermax_cal=20;
00631 
00632   long ncoefcont, ncoefcal, npts_poly;
00633 
00634   set_purpose("Wavelength calibration (computation)");
00635   set_arglist("-in Cd -reftable Ref_Cd -N|halfwidth 6 -sigmax 2.0 -sigmin 0.3 "
00636               "-lens|lensno all -dbgtable dbg_wcalib -dbgspec dbg_gaussfit "
00637               "-degcont 3 -degcal 3 -W 0");
00638 
00639   init_snifs("$Name:  $");
00640   init_session(argv,argc,&arglabel,&argval);
00641 
00642   if (DEBUG) {
00643     print_msg("$Id: wcalib.c,v 1.16 2004/11/19 12:55:14 ycopin Exp $");
00644     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00645   }
00646 
00647   inname  = argval[0];
00648   refname = argval[1];
00649   get_argval(2,"%lf",&Nhw);                  /* Half-width line window */
00650 
00651   get_argval(3,"%lf",&sigmax);
00652   get_argval(4,"%lf",&sigmin);
00653 
00654   lensno = -1;
00655   if (DEBUG) {
00656     if (!strcmp(argval[5],"all")) lensno = -1;
00657     else get_argval(5,"%d",&lensno);
00658     dbgtable = argval[6];
00659     dbgspec  = argval[7];
00660   }
00661 
00662   get_argval(8,"%ld",&ncoefcont); ncoefcont++; /* Degree -> Nb of coeff. */
00663   get_argval(9,"%ld",&ncoefcal);  ncoefcal++;  /* Degree -> Nb of coeff. */
00664   if (ncoefcal > 9) {          /* There's not enough room in column name */ 
00665     print_error("Degree of calibration is too high! (<8)");
00666     exit_session(UNKNOWN);
00667   }
00668   npts_poly = ncoefcal * 2;
00669 
00670   get_argval(10,"%lf",&Wbound);              /* minimum distance (in sigma) 
00671                                                 of a line to the edges of the spectrum */
00672 
00673   /*======================================================================*/
00674   /* Ouverture de la table des raies d'emission de reference */
00675   /*======================================================================*/
00676 
00677   print_msg("Opening wavelength reference table %s", refname);
00678   if (open_table(&reftable,refname,"I")<0) {
00679     print_error("Unable to open reference table %s", refname);
00680     exit_session(ERR_OPEN);
00681   }
00682 
00683   /* ========== Catalogue ========== */
00684 
00685   CAT = (strstr(inname,".cat") != NULL);
00686   if (CAT) {
00687     print_msg("Calibration datacube catalog: %s", inname);
00688     RD_catalog(inname, cubename);
00689   } else {
00690     strcpy(cubename, inname);
00691   } 
00692 
00693   /* ========== Calcul sur tout le catalogue si necessaire ========== */
00694 
00695   do {
00696 
00697     NbTotGauss = 0;
00698 
00699     /*======================================================================*/
00700     /*                  Ouverture du cube TIGER                         */
00701     /*======================================================================*/
00702 
00703     print_msg("Opening datacube %s", cubename);
00704     if (open_tiger_frame(&cube,cubename,"I")<0){
00705       print_error("Unable to access the datacube");
00706       exit_session(ERR_OPEN);
00707     }
00708 
00709     switch ((class = read_file_class(&cube))) {
00710     case SUPER_CLASS:
00711     case RAW_CAL_CUBE: break;
00712     default:
00713       print_error("Wrong class type %d: %s should be "
00714                   "`Raw Calibration Datacube'",class,cube.name);
00715 
00716       print_msg("DEBUG: Pb with datacube FCLASS, keep going...");
00717       //      exit_session(ERR_BAD_TYPE);
00718     }
00719     nolens=(int *)malloc(sizeof(int)  *cube.nbspec);
00720     il =   (int *)malloc(sizeof(int)  *cube.nbspec);
00721     x0 = (float *)malloc(sizeof(float)*cube.nbspec);
00722     y0 = (float *)malloc(sizeof(float)*cube.nbspec);
00723 
00724     /* get_lenses_no(&cube,nolens); */
00725     get_lenses_coord(&cube,"X0","Y0",nolens,x0,y0,il);
00726 
00727     /*======================================================================*/
00728     /* Ouverture de la table associee au cube */
00729     /*======================================================================*/
00730 
00731     if (lensno == -1) {
00732       if (open_table(&tassoc,cube.table_name,"IO")<0){
00733         close_tiger_frame(&cube);
00734         print_error("Unable to access the table associated with the datacube");
00735         exit_session(ERR_OPEN);
00736       }
00737       switch ((class = read_file_class(&tassoc))) {
00738       case SUPER_CLASS:
00739       case TBL_FILE: break;
00740       case TBL_FILE_CAL:
00741         print_warning("Computating new calibration for table %s", tassoc.name);
00742         break;
00743       default:
00744         print_error("Wrong class type %d: %s should be a "
00745                     "`[Calibrated] table'",class,tassoc.name);
00746         exit_session(ERR_BAD_TYPE);
00747       }
00748       col_no = get_col_ref(&tassoc, LAB_COL_NO);
00749     } else {
00750       print_warning("DEBUG: table %s untouched", cube.table_name);
00751     }
00752 
00753     /*======================================================================*/
00754     /* retrieve filter bounds in the datacube descriptor */
00755     /*======================================================================*/
00756 
00757     if (RD_desc(&cube,LRANGE4,DOUBLE,4,dbl_desc) < 0) {
00758       print_error("Cannot read descriptor %s from datacube",LRANGE4);
00759 
00760       print_warning("Previous versions of extract_spec (bef. 10/11/04, v.1.31)\n"
00761                     "were writing an invalidly-long keyword LBDARANGn (n=1...4).\n"
00762                     "See if you can find it back and rename it to LRANGEn.\n");
00763       exit_session(ERR_NODESC);
00764     }
00765     else {
00766       lborne1 = dbl_desc[0];
00767       lborne2 = dbl_desc[1];
00768       delta_borne = lborne2 - lborne1;
00769       print_msg("Spectral domain: %.2f-%.2f A",lborne1,lborne2);
00770     }
00771 
00772     /*======================================================================*/
00773     /* convert sigma bounds to A */
00774     /*======================================================================*/
00775     sigmax *= cube.step;
00776     sigmin *= cube.step;
00777 
00778     /*======================================================================*/
00779     /* working on reference table */
00780     /*======================================================================*/
00781 
00782     print_msg("Extracting lines from reference table %s", reftable.name);
00783     LRefSize = Take_Ref(&reftable,&LRef,lborne1,lborne2);
00784     print_msg("Found %d lines between %.1fA and %.1fA",
00785               LRefSize, lborne1, lborne2);
00786    
00787     /*======================================================================*/
00788     /* Regroup peaks according to proximity */
00789     /*======================================================================*/
00790    
00791     RefPeak = Join_Ref_Peak(LRefSize,LRef,Nhw,cube.step);
00792 
00793     /*======================================================================*/
00794     /* allocation memoire */
00795     /*======================================================================*/
00796     X         = (double *)malloc(LRefSize * sizeof(double));
00797     dX        = (double *)malloc(LRefSize * sizeof(double));
00798     Weightcal = (double *)calloc(LRefSize,  sizeof(double));
00799     LAMBDA    = (double *)malloc(LRefSize * sizeof(double));
00800     I         = (double *)malloc(LRefSize * sizeof(double));
00801     SIGMA     = (double *)malloc(LRefSize * sizeof(double));
00802     FLAG      =    (int *)malloc(LRefSize * sizeof(int));
00803     LRefcal   = (double *)malloc(npts_poly * sizeof(double));
00804     dXcal     = (double *)malloc(npts_poly * sizeof(double));
00805     coefcont  = (double *)malloc(ncoefcont * sizeof(double));
00806     coefcal   = (double *)malloc(ncoefcal  * sizeof(double));
00807     errcal    = (double *)malloc(cube.nbspec * sizeof(double));
00808 
00809     /*======================================================================*/
00810     /* creation de la table  de debug et des colonnes */
00811     /*======================================================================*/
00812 
00813     if (DEBUG) {
00814       print_msg("Creating debug table %s", dbgtable);
00815       if (create_table(&table,dbgtable,-1,7,'Q',"Wcalib table")<0) {
00816         print_error("Unable to create the debug table...");
00817         exit_session(ERR_CREAT);
00818       }
00819       else {
00820         id_no  = create_col(&table,LAB_COL_NO,LONG, 'O',"I4",  "Lens Number");
00821         id_x0  = create_col(&table,"X0",      FLOAT,'O',"F9.6","px");
00822         id_y0  = create_col(&table,"Y0",      FLOAT,'O',"F9.6","px");
00823         id_lref= create_col(&table,"REF",     FLOAT,'O',"F9.6","A");
00824         id_x   = create_col(&table,"X",       FLOAT,'O',"F9.6","A");
00825         id_i   = create_col(&table,"I",       FLOAT,'O',"F9.6","e-");
00826         id_sig = create_col(&table,"SIG",     FLOAT,'O',"F9.6","A");
00827         id_dx  = create_col(&table,"DX",      FLOAT,'O',"F9.6","A");
00828         id_fit = create_col(&table,"XFIT",    FLOAT,'O',"F9.6","A");
00829         id_ifail     = create_col(&table,"IFAIL",LONG, 'O',"I2","");
00830         id_weightcal = create_col(&table,"WCAL", FLOAT,'O',"F3.1","number");
00831       }
00832     }
00833 
00834     k_min = -1;
00835 
00836     /*======================================================================*/
00837     /* On commence par une lentille */
00838     /*======================================================================*/
00839     if (lensno!=-1) { /* mode debug, une seule lentille demande */
00840 
00841       j=0;
00842       jstart=0;
00843       nbgoodcal = 0;
00844     
00845       print_msg("Starting computation for lens %d",lensno);
00846     
00847       get_tiger_spec(&cube,&spectre,NULL,lensno);
00848 
00849       /*======================================================================*/
00850       /* Ajustement du continu sous les raies: pour cela on met
00851          une rejection importante positive et un nombre d'iterations
00852          de 10 */
00853       /*======================================================================*/
00854       start_polcont = spectre.start;
00855       end_polcont   = spectre.end;
00856       fit_poly_rej_nag(&spectre, ncoefcont, siginf_cont, sigsup_cont, 
00857                        itermax_cont, coefcont,&niter,&error,&nrej,start_polcont,
00858                        end_polcont, &status);
00859 
00860       courant = RefPeak; k = 0; 
00861 
00862       print_msg("Scanning the calibration lines set");
00863       while (courant!=NULL) {
00864 
00865         /*====================================================================*/
00866         /* attention: Get_Spectrum_Data inclut la soustraction du continu
00867            qui a ete ajuste ci-dessus */
00868         /*====================================================================*/
00869 
00870         NbPixel = (long)Get_Spectrum_Data(&spectre, courant->lbda_inf, 
00871                                           courant->lbda_sup, coefcont, 
00872                                           ncoefcont, start_polcont,end_polcont);
00873 
00874         if (NbPixel < 5) {
00875           k += courant->nb;
00876           courant = courant->next;
00877           continue;
00878         }
00879 
00880         lbdastep = spectre.step;
00881         NbGauss  = courant->nb;
00882         NbTotGauss += NbGauss;
00883 
00884         /*====================================================================*/
00885         /* ajustement des gaussiennes */
00886         /*====================================================================*/
00887 
00888         if (k_min < 0) k_min = k;;
00889 
00890         Detection_peak(NbPixel, Wbound, &(X[k]), &(LRef[k]), &(SIGMA[k]),
00891                        &(I[k]), &(FLAG[k]), sigmin, sigmax);
00892 
00893         for (i=0;i < NbGauss;i++) {
00894           dX[k+i] = X[k+i] - LRef[k+i];
00895 
00896           /*==================================================================*/
00897           /* si l'ajustement de cette raie n'a pas marche on elimine */
00898           /*==================================================================*/
00899 
00900           if (FLAG[k+i] <= 1) /* BUG ???? <= 1 ou < 1 ??? */
00901             Weightcal[k+i] = 1.;
00902           else 
00903             Weightcal[k+i] = 1.e-99;
00904         }
00905      
00906         if (DEBUG) {
00907           for (i=0; i<NbGauss; i++) {
00908             WR_tbl(&table,j,id_no,&lensno);
00909             valfloat = LRef[k+i];  WR_tbl(&table,j,id_lref,&valfloat);
00910             valfloat = X[k+i];     WR_tbl(&table,j,id_x,  &valfloat);
00911             valfloat = I[k+i];     WR_tbl(&table,j,id_i,  &valfloat);
00912             valfloat = SIGMA[k+i]; WR_tbl(&table,j,id_sig,&valfloat);
00913             valfloat = dX[k+i];    WR_tbl(&table,j,id_dx, &valfloat);
00914             WR_tbl(&table,j,id_ifail,&(FLAG[k+i]));
00915             j++;
00916           }
00917         }
00918 
00919         k += courant->nb;
00920         courant=courant->next;
00921       }
00922 
00923       if (DEBUG) {
00924         print_msg("Saving fit result as debug spectrum file");
00925         disable_erase_flag();
00926         create_spec(&resultat,dbgspec,spectre.npts,spectre.start,
00927                     spectre.step,DOUBLE,"essai","A");
00928 
00929         SaveGauss(&resultat,X+k_min,SIGMA+k_min,I+k_min, 
00930                   NbTotGauss, coefcont, ncoefcont, 
00931                   start_polcont, end_polcont);
00932         write_file_class(&resultat, DEBUG_FILE);
00933         close_spec(&resultat);
00934       }
00935 
00936       /*======================================================================*/
00937       /* On determine les bornes a passer au polynome */
00938       /* ceci au cas ou le centre d'une raie ajustee tomberait hors spectre */
00939       /*======================================================================*/
00940       start_polcal = MIN(lborne1, LRef[0]);
00941       end_polcal = MAX(lborne2, LRef[LRefSize-1]);
00942 
00943       /*======================================================================*/
00944       /* Ajustement du polynome de calibration sur les raies ajustees */
00945       /*======================================================================*/
00946       if (NbTotGauss < ncoefcal) {
00947         flag = 3;
00948         print_error("Spectrum %d cannot be calibrated: lack of lines to fit",
00949                     lensno);
00950         exit_session(UNKNOWN);
00951       }
00952       else {
00953         fit_poly_rej_nag_tab(LRef+k_min,dX+k_min,Weightcal+k_min,
00954                              (long)NbTotGauss, 
00955                              ncoefcal,siginf_cal,sigsup_cal,itermax_cal,coefcal,
00956                              &niter, &error, &nrej, start_polcal, end_polcal, 
00957                              &status);
00958       }
00959 
00960       if (DEBUG) {
00961         for (i=0; i<NbTotGauss; i++) {
00962           if (Weightcal[i+k_min] < 1.e-20) valfloat = 0.;
00963           else                             valfloat = Weightcal[i+k_min];
00964           WR_tbl(&table,i,id_weightcal,&valfloat);
00965         }
00966       }
00967 
00968       /*======================================================================*/
00969       /*  Calcul du flag 0, 1 ou 2 du spectre courant */
00970       /*======================================================================*/
00971 
00972       switch(status) {
00973       case -1:
00974         print_error("Spectrum %d not calibrated: bounds error",lensno);
00975         exit_session(UNKNOWN);
00976       case -2:
00977         flag = 3;
00978         print_error("Spectrum %d not calibrated: lack of lines to fit",lensno);
00979         exit_session(UNKNOWN);
00980       default:
00981         nbgoodcal = nbgoodcal + 1;
00982         if (DEBUG) 
00983           for (i=0; i<NbTotGauss; i++) {
00984             valfloat = val_poly_nag(LRef[i+k_min], coefcal, ncoefcal, start_polcal,
00985                                     end_polcal, &val_status);
00986             WR_tbl(&table,i,id_fit,&valfloat);
00987           }
00988         break;
00989 
00990       }
00991 
00992       error_norm = error / lbdastep;         /* In pixels */
00993       perc_nrej = (double)nrej / (double)NbTotGauss;
00994       npts_cal = NbTotGauss - nrej;
00995       ndiff = npts_cal - ncoefcal;
00996     
00997       if (ndiff < 0) {
00998         flag = 3;
00999         print_error("Spectrum %d not calibrated: lack of lines to fit",lensno);
01000         exit_session(UNKNOWN);
01001       }
01002 
01003       flag = 0;
01004       if ((error_norm > CRITER_SIG_2) || 
01005           (perc_nrej > CRITER_NREJ_2) || 
01006           (ndiff < CRITER_DIFF_2)) {
01007         flag = 2;
01008       }
01009       else if ((error_norm > CRITER_SIG_1) || 
01010                (perc_nrej > CRITER_NREJ_1) || 
01011                (ndiff < CRITER_DIFF_1)) {
01012         flag = 1;
01013       }
01014 
01015       /*======================================================================*/
01016       /* Ecriture des resultats dans la table */
01017       /*======================================================================*/
01018       print_msg("%d lines in total: %d used, %d rejected, residual: %f",
01019                 NbTotGauss, npts_cal, nrej, error_norm);
01020       switch (flag) {
01021       case 0: print_msg("Fit quality: excellent"); break;
01022       case 1: print_msg("Fit quality: very good"); break;
01023       case 2: print_msg("Fit quality: good");      break;
01024       case 3: print_msg("Fit quality: poor");      break;
01025       }
01026 
01027       free_spec_mem(&spectre);
01028 
01029     }
01030     /*======================================================================*/
01031     /* maintenant pour toutes les lentilles */
01032     /*======================================================================*/
01033     else {
01034 
01035       nbligne = 0;
01036 
01037       /*======================================================================*/
01038       /* Creation des colonnes de sortie pour la table associee au cube */
01039       /*======================================================================*/
01040 
01041       id_xmin = create_col(&tassoc,"XMIN",DOUBLE,'R',"F9.6","A");
01042       id_xmax = create_col(&tassoc,"XMAX",DOUBLE,'R',"F9.6","A");
01043       id_deg  = create_col(&tassoc,"DEGREE",INT,'R',"I3","number");
01044 
01045       id_c      =   (int *)malloc(ncoefcal * sizeof(int));
01046       namecol_c = (char **)malloc(ncoefcal * sizeof(char *));
01047       for (i=0; i<ncoefcal; i++) {
01048         namecol_c[i] = (char *)malloc(3 * sizeof(char)); /* ncoefcal up to 9 */
01049         sprintf(&(namecol_c[i][0]),"C%d",i);
01050         id_c[i] = create_col(&tassoc,namecol_c[i],DOUBLE,'R',"F9.6","number");
01051       }
01052       id_nrej = create_col(&tassoc,"NREJ", INT,'R',"I4","number");
01053       id_error= create_col(&tassoc,"ERROR",DOUBLE,'R',"F9.6","px");
01054       id_npts = create_col(&tassoc,"NFIT",INT,'R',"I4","number");
01055       id_flag = create_col(&tassoc,"FLAG",INT,'R',"I4","number");
01056 
01057       print_msg("Starting computation for all lenses");
01058 
01059       reset_print_progress();
01060       for (i=0;i<cube.nbspec;i++) {        
01061         get_tiger_spec(&cube,&spectre,NULL,nolens[i]);
01062 
01063         print_progress("Calibration in progress",(i+1.)/cube.nbspec*100.,1.0);
01064               
01065         /*====================================================================*/
01066         /* Ajustement du continu sous les raies: pour cela on met
01067            une rejection importante positive et un nombre d'iterations
01068            de 10 */
01069         /*====================================================================*/
01070         start_polcont = spectre.start;
01071         end_polcont = spectre.end;
01072         fit_poly_rej_nag(&spectre, ncoefcont, siginf_cont, sigsup_cont, 
01073                          itermax_cont, coefcont, &niter, &error, &nrej, 
01074                          start_polcont, end_polcont, &status);
01075 
01076         courant = RefPeak; k = 0; NbTotGauss = 0;
01077 
01078         while (courant!=NULL) {
01079 
01080           /*==================================================================*/
01081           /* attention: Get_Spectrum_Data inclus la soustraction du continu
01082              qui a ete ajuste ci-dessus */
01083           /*==================================================================*/
01084           NbPixel = Get_Spectrum_Data(&spectre, courant->lbda_inf, 
01085                                       courant->lbda_sup,  coefcont, 
01086                                       ncoefcont, start_polcont, 
01087                                       end_polcont);
01088           if (NbPixel < 5) {
01089             k += courant->nb;
01090             courant = courant->next;
01091             continue;
01092           }
01093 
01094           lbdastep = spectre.step;
01095           NbGauss = courant->nb;
01096           NbTotGauss += NbGauss;
01097 
01098 
01099           /*==================================================================*/
01100           /* ajustement des gaussiennes */
01101           /*==================================================================*/
01102           
01103           if (k_min < 0) k_min = k;
01104 
01105           Detection_peak(NbPixel, Wbound, &(X[k]), &(LRef[k]), &(SIGMA[k]),
01106                          &(I[k]), &(FLAG[k]), sigmin, sigmax);
01107 
01108           for (j=0; j<NbGauss; j++) { 
01109             dX[k+j] = X[k+j] - LRef[k+j];
01110 
01111             /*================================================================*/
01112             /* si l'ajustement de cette raie n'a pas marche on elimine */
01113             /*================================================================*/
01114             if (FLAG[k+j] <= 1)   /* BUG ???? <= 1 ou < 1 ??? */
01115               Weightcal[k+j] = 1.;
01116             else
01117               Weightcal[k+j] = 1.e-99;
01118           }
01119 
01120           if (DEBUG) {
01121             for (j=0; j<NbGauss; j++) {
01122               WR_tbl(&table,nbligne,id_no,&(nolens[i]));
01123               WR_tbl(&table,nbligne,id_x0,&(x0[i]));
01124               WR_tbl(&table,nbligne,id_y0,&(y0[i]));
01125               WR_tbl(&table,nbligne,id_lref,&(LRef[k+j]));
01126               WR_tbl(&table,nbligne,id_x,&(X[k+j]));
01127               WR_tbl(&table,nbligne,id_i,&(I[k+j]));
01128               WR_tbl(&table,nbligne,id_sig,&(SIGMA[k+j]));
01129               WR_tbl(&table,nbligne,id_dx,&(dX[k+j]));
01130               WR_tbl(&table,nbligne,id_ifail,&(FLAG[k+i]));
01131               nbligne++;
01132             }
01133           }
01134           k += courant->nb;
01135 
01136           courant=courant->next;
01137         }    
01138        
01139         /*====================================================================*/
01140         /* On determine les bornes a passer au polynome */
01141         /*====================================================================*/
01142         start_polcal = MIN(lborne1, LRef[0]);
01143         end_polcal = MAX(lborne2, LRef[LRefSize-1]);
01144 
01145         /*====================================================================*/
01146         /* Ajustement du polynome sur les raies ajustees */
01147         /*====================================================================*/
01148         if (NbTotGauss < ncoefcal) {
01149           flag = 3;
01150           WR_tbl(&tassoc, il[i], id_flag, &flag);
01151           print_warning("Spectrum %d not calibrated: lack of lines to fit",
01152                         nolens[i]);
01153           fflush(stdout);
01154           continue;
01155         }
01156         else {
01157           fit_poly_rej_nag_tab(LRef+k_min, dX+k_min, Weightcal+k_min, 
01158                                (long)NbTotGauss, ncoefcal, 
01159                                siginf_cal, sigsup_cal, itermax_cal, coefcal, 
01160                                &niter, &error, &nrej, start_polcal, end_polcal,
01161                                &status);
01162         }
01163 
01164         /*====================================================================*/
01165         /*      Calcul du flag 0, 1 ou 2 du spectre courant */
01166         /*====================================================================*/
01167         switch(status) {
01168         case -1:
01169           flag = 3;
01170           WR_tbl(&tassoc, il[i], id_flag, &flag);
01171           print_warning("Spectrum %d not calibrated: error on bounds",
01172                         nolens[i]);
01173           fflush(stdout);
01174           continue;
01175         case -2:
01176           flag = 3;
01177           WR_tbl(&tassoc, il[i], id_flag, &flag);
01178           print_warning("Spectrum %d not calibrated: too many lines rejected",
01179                         nolens[i]);
01180           fflush(stdout);
01181           continue;
01182         default:
01183           nbgoodcal = nbgoodcal + 1;
01184           break;
01185 
01186         }
01187 
01188         error_norm = error / lbdastep;
01189         perc_nrej = (double)nrej / (double)NbTotGauss;
01190         npts_cal = NbTotGauss - nrej;
01191         ndiff = npts_cal - ncoefcal;
01192 
01193         if (ndiff < 0) {
01194           flag = 3;
01195           nflag3++;
01196           WR_tbl(&tassoc, il[i], id_flag, &flag);
01197           print_warning("Spectrum %d not calibrated: lack of lines to fit",
01198                         nolens[i]);
01199           continue;
01200         }
01201 
01202         flag = 0;
01203         if ((error_norm > CRITER_SIG_2) ||
01204             (perc_nrej  > CRITER_NREJ_2) ||
01205             (ndiff      < CRITER_DIFF_2)) {
01206           flag = 2;
01207           nflag2++;
01208         }
01209         else if ((error_norm > CRITER_SIG_1) || 
01210                  (perc_nrej  > CRITER_NREJ_1) ||
01211                  (ndiff      < CRITER_DIFF_1)) {
01212           flag = 1;
01213           nflag1++;
01214         }
01215         errcal[nb_error++] = error_norm;
01216 
01217         /*====================================================================*/
01218         /* Ecriture des resultats dans la table */
01219         /*====================================================================*/
01220         WR_tbl(&tassoc, il[i], id_xmin, &start_polcal);
01221         WR_tbl(&tassoc, il[i], id_xmax, &end_polcal);
01222         WR_tbl(&tassoc, il[i], id_nrej, &nrej);
01223         WR_tbl(&tassoc, il[i], id_npts, &npts_cal);
01224         WR_tbl(&tassoc, il[i], id_error,&error_norm);
01225         WR_tbl(&tassoc, il[i], id_flag, &flag);
01226         varlong = ncoefcal - 1;
01227         WR_tbl(&tassoc, il[i], id_deg, &varlong);
01228         for (j=0;j < ncoefcal;j++) {
01229           vardouble = coefcal[j];
01230           WR_tbl(&tassoc, il[i], id_c[j], &vardouble);
01231         }
01232         free_spec_mem(&spectre);
01233       }
01234       WR_desc(&tassoc, WLMAXDEGRE, INT, 1, &ncoefcal);
01235       /* write type of table (wavelength calibrated table) */
01236       write_file_class(&tassoc, TBL_FILE_CAL);
01237 
01238       /* statistic on error */
01239       err_mean = gsl_stats_mean(errcal,1,nb_error);
01240       err_std  = gsl_stats_sd_m(errcal,1,nb_error,err_mean);
01241       gsl_stats_minmax(&err_min,&err_max,errcal,1,nb_error);
01242 
01243       print_msg("Calibration Mean Error: %5.3f px +- %5.3f px [%5.3f - %5.3f]",
01244                 err_mean, err_std, err_min, err_max);
01245       rdesc = err_mean;
01246       WR_desc(&tassoc, "MEANERR", FLOAT, 1, &rdesc);
01247       rdesc = err_std;
01248       WR_desc(&tassoc, "STDERR", FLOAT, 1, &rdesc);
01249       rdesc = err_min;
01250       WR_desc(&tassoc, "MINERR", FLOAT, 1, &rdesc);
01251       rdesc = err_max;
01252       WR_desc(&tassoc, "MAXERR", FLOAT, 1, &rdesc);
01253     }
01254 
01255     /* ArcGroup_Delete(RefPeak); */
01256     free(LRef);
01257 
01258     close_tiger_frame(&cube);
01259     if (lensno == -1) close_table(&tassoc);
01260     if (DEBUG) {
01261       write_file_class(&table, DEBUG_FILE);
01262       close_table(&table);
01263     }
01264 
01265     /* Nombre total de spectres etalonnes avec succes, en mode non-DEBUG */
01266 
01267     if (!DEBUG) 
01268       print_msg("-> %d spectra successfully calibrated",nbgoodcal);
01269 
01270     /* Liberation de la memoire */
01271         
01272     free(nolens);
01273     free(X); free(dX);
01274     free(Weightcal); free(LRefcal); free(dXcal);
01275     free(LAMBDA); free(I); free(SIGMA); free(FLAG);
01276     free(coefcont); free(coefcal);
01277 
01278     /* Fin de la boucle sur le catalogue */
01279   } while (CAT && RD_catalog(inname, cubename));
01280 
01281   close_table(&reftable);
01282 
01283   exit_session(OK);
01284   return(0);
01285 }    
01286 
01287 /* 
01288 ;;; Local Variables: ***
01289 ;;; eval: (add-to-list 'c-font-lock-extra-types "ArcLine") ***
01290 ;;; eval: (add-to-list 'c-font-lock-extra-types "ArcGroup") ***
01291 ;;; End: ***
01292 */

Generated on Tue Nov 23 18:04:18 2004 for Snifs by doxygen 1.3.3