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

extract/source/create_mask.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00014 /* =========================================================== */
00015 
00016 #include <IFU_io.h>
00017 #include <IFU_math.h>
00018 #include <snifs.h>
00019 #include <extract.h>
00020 
00021 #include <gsl/gsl_statistics.h>
00022 
00023 //#define LOC_NCDX  3                   //!< Nb of coeff. in dx(lambda) (<9)
00024 //#define LOC_NCSIG 3                   //!< Nb of coeff. in sigma(lambda) (<9)
00025 //#define LOC_CLDX  3                   //!< Sigma-clipping for dx
00026 //#define LOC_CLSIG 3                   //!< Sigma-clipping for sigma
00027 //#define LOC_MXDX  3                   //!< Max iter for dx
00028 //#define LOC_MXSIG 3                   //!< Max iter for sigma
00029 static int const LOC_NCDX= 3;           
00030 static int const LOC_NCSIG=3;           
00031 static int const LOC_CLDX= 3;           
00032 static int const LOC_CLSIG=3;           
00033 static int const LOC_MXDX= 3;           
00034 static int const LOC_MXSIG=3;           
00035 
00036 #define MAXDX (0.5*SNIFS_INTERSP)
00037 
00038 SnifsOptics *gloptics;
00039 Maxima_Set *glmaxset;
00040 IMAGE2D *glframe;
00041 
00042 int glFITARC,glFITCONT, glnlbda;
00043 
00044 float *glxppup,*glyppup;
00045 double glpixsize,gllbdapup, *gllbda, glnormmax;  
00046 
00047 /* === Doxygen Comment ======================================= */
00068 /* =========================================================== */
00069 
00070 long nllsqfit_mask(long *mode, long *nlenses, long *npar, long *ldfj, 
00071                    double *par, double f[], double fjac[], long *nstate, 
00072                    long *iuser, double user[])
00073 {
00074   /* Nb of iterations and best chi2 */
00075   static int niter = 0;
00076   static double best_merit = MAXDOUBLE;
00077   static int maxstep = 0;
00078   static float lstep = 0;
00079 
00080   SnifsOptics *optics;
00081   Maxima_Set *maxset;
00082   IMAGE2D *frame;
00083   float *xppup, *yppup;
00084   double pixsize,lbdapup, *tlbda;  
00085 
00086   char text[132];
00087 
00088   int i,j, order;
00089   int nstep, iline,ipeak;
00090 
00091   double tpar[16];
00092   double farc,fcont;
00093   double merit,meritarc,meritcont, xmla,ymla, xccd,yccd, I,dx, lbda;
00094 
00095   /* ===== Initialization ============================== */
00096 
00097   /* Convert global variables to local */
00098   optics  = gloptics;
00099   maxset  = glmaxset;
00100   frame   = glframe;
00101   pixsize = glpixsize;
00102   lbdapup = gllbdapup;
00103   xppup   = glxppup;
00104   yppup   = glyppup;
00105   tlbda   = gllbda;
00106 
00107   /* Maxset step [px] and lambda step [AA = (AA/mm)*(mm/px)*px] */
00108   if (!niter) {
00109     maxstep = (maxset->line[maxset->nb_ycoords-1].ycoord - 
00110                maxset->line[0].ycoord) / (maxset->nb_ycoords - 1);
00111     /* If FITCONT = n, fit every n-th maxima line */
00112     lstep = optics->grism.w_disp*pixsize*maxstep * glFITCONT;
00113   }
00114 
00115   /* Unnormalize mask parameters and set optical parameters */
00116   memcpy((double *)tpar, (double *)par, 16*sizeof(double));
00117   unnorm_mask_param(optics->channel,tpar);
00118   mask_param2optics(tpar,optics); 
00119 
00120   niter++;
00121 
00122   /* ===== Merit functions ============================== */
00123 
00124   meritarc = meritcont = merit = 0;
00125 
00126   for (i=0; i<*nlenses; i++) { /* Loop over lenses */
00127 
00128     /* [CCD/px] -> [MLA/mm] */
00129     snifs_optics_CCD2MLA(optics, xppup[i],yppup[i], &xmla,&ymla, 
00130                          pixsize,lbdapup);
00131     
00132     /* ----- Arc frame ------------------------------ */
00133 
00134     if (glFITARC) {
00135 
00136       farc = 0;
00137 
00138       /* ..... 1st-order contribution .............................. */
00139 
00140       for (j=0; j<glnlbda; j++) { /* loop over arc lines */
00141 
00142         /* [MLA/mm] -> [CCD/px] */
00143         snifs_optics_MLA2CCD(optics, xmla,ymla, &xccd,&yccd, 
00144                              pixsize,tlbda[j],1);
00145 
00146         if (in_frame(frame,xccd,yccd)) {
00147           /* Interpolated intensity */
00148           I = MAX(0,RD_frame_interp(frame, xccd,yccd)) / glnormmax;
00149           
00150           /* Merit function */
00151           farc += exp(-I);
00152         }
00153           
00154       } /* j: loop over arc lines */
00155 
00156       /* ..... Higer orders contribution .............................. */
00157 
00158       /* Add other orders contribution (including blaze function) if needed */
00159 
00160     }
00161     else farc = 1;
00162     
00163     /* ----- Maxima ------------------------------ */
00164 
00165     if (glFITCONT) { 
00166 
00167       fcont = 0;
00168       
00169       nstep = 0;
00170 
00171       /* ..... 1st-order contribution .............................. */
00172 
00173       for (lbda  = optics->filter.inf_util; 
00174            lbda <= optics->filter.sup_util; 
00175            lbda += lstep) { /* loop over lbda */
00176 
00177         /* [MLA/mm] -> [CCD/px] */
00178         snifs_optics_MLA2CCD(optics, xmla,ymla, &xccd,&yccd, 
00179                              pixsize,lbda,1);
00180 
00181         if (in_frame(frame,xccd,yccd)) {
00182           /* Find nearest max position */
00183           if (find_nearest_max(maxset, xccd,yccd, &iline,&ipeak))
00184             dx = MAXDOUBLE;
00185           else 
00186             dx = ABS(xccd - maxset->line[iline].maxima[ipeak].xcoord);
00187           dx = MIN(dx,MAXDX);
00188 
00189           /* Merit function */
00190           if (dx < MAXDX) { /* Only on good points */
00191             fcont += SQ(dx);
00192             nstep ++;
00193           }
00194         }
00195     
00196       } /* lbda: loop over wavelength */
00197       
00198       /* ..... Higer orders contribution .............................. */
00199       
00200 //      for (order=0; order<=2; order+=2) { /* Oth & 2nd orders */
00201 //
00202 //        for (lbda  = optics->filter.inf_util; 
00203 //             lbda <= optics->filter.sup_util; 
00204 //             lbda += lstep) { /* loop over lbda */
00205 //
00206 //          /* [MLA/mm] -> [CCD/px] */
00207 //          snifs_optics_MLA2CCD(optics, xmla,ymla, &xccd,&yccd, 
00208 //                               pixsize,lbda,order);
00209 //
00210 //          if (in_frame(frame,xccd,yccd)) {
00211 //            /* Find nearest max position */
00212 //            if (find_nearest_max(maxset, xccd,yccd, &iline,&ipeak))
00213 //              dx = MAXDOUBLE;
00214 //            else 
00215 //              dx = ABS(xccd - maxset->line[iline].maxima[ipeak].xcoord);
00216 //            dx = MIN(dx,MAXDX);
00217 //
00218 //            /* Merit function */
00219 //            if (dx < MAXDX) { /* Only on good points */
00220 //              fcont += SQ(dx);
00221 //              nstep ++;
00222 //            }
00223 //          }
00224 //    
00225 //        } /* lbda: loop over wavelength */
00226 //      
00227 //      } /* orders */
00228 
00229       if (nstep) fcont /= nstep;
00230 
00231     }
00232     else fcont = 1.;
00233 
00234     /* ----- Global merit function ------------------------------ */
00235 
00236     f[i] = farc*fcont;
00237 
00238     meritarc  += farc;
00239     meritcont += fcont;
00240     merit     += f[i];
00241     
00242   } /* i: loop over lenses */
00243 
00244   if (*nstate || (!(*nstate) && merit < 0.99*best_merit)) {
00245     /* Printing results */
00246     sprintf(text,"Call %5d: Arc = %10.4e, Cont = %10.4e => Merit = %10.4e",
00247             niter,meritarc,meritcont,merit);
00248     print_label(text);
00249     if (DEBUG) print_mask_param(tpar);
00250     best_merit = merit;
00251   }
00252 
00253   if (!glFITARC && !glFITCONT) *mode = -1; /* Stop here if no fit */
00254   
00255   return(OK);
00256 }
00257 
00258 /* === Doxygen Comment ======================================= */
00310 /* =========================================================== */
00311 
00312 int main(int argc, char **argv)
00313 {
00314   SnifsConfig config;
00315   SnifsOptics optics;
00316 
00317   Maxima_Set maxset;
00318   IMAGE2D arcframe;   
00319 
00320   TABLE mask,reftab;
00321 
00322   char **argval, **arglabel;
00323   char *maskname, *maxname, *arcname, *refname;
00324   char tmpname[lg_name+1];
00325 
00326   int *nopup, *iflag;
00327   int i;
00328   int FITARC,FITCONT,SKIPFIT,DEFAULT,BLAZE,NOLOCAL;
00329   int npup, nlbda, npar, status;
00330   int colno,colxd,colyd,collbda,colmsk;
00331   
00332   float *xppup,*yppup, *tmplbda;
00333   float lbdapup;
00334 
00335   double *par,*epar,*parlb,*parub, *lbda;
00336   double tol, lbda_inf,lbda_sup, objf,gof, mean,sigma;
00337   
00338   SPECTRUM specS,specN;
00339   TIGERfile cubesig, cubedist;
00340 
00341   int j,k, nmax, ifail, niter, nrej;
00342   int maxstep;
00343 
00344   float fval;
00345 
00346   double *plbda, *pdx, *psigma;
00347   double varbef, varaft, dxerr, sigerr;
00348   double lstep;
00349 
00350   int ngpup   = 0;
00351   int nmaxtot = 0;
00352   double varbeftot  = 0;
00353   double varafttot  = 0;
00354 
00355   int dxnc  = LOC_NCDX;                      /* Nb of coeff. in dx(lbda) */
00356   int signc = LOC_NCSIG;                     /* Nb of coeff. in sig(lbda) */
00357   double dxcoeff[LOC_NCDX];
00358   double sigcoeff[LOC_NCSIG];
00359 
00360   int coldx[LOC_NCDX], colddx;               /* Column index */
00361   int colsig[LOC_NCSIG], coldsig; 
00362 
00363   /* ##### INTRODUCTION ############################## */
00364 
00365   set_purpose("Creates the extraction mask according to the optical model "
00366               "from max and arc frames");
00367   set_arglist("-mask mask -max max -cal arc -reftable arc_ref "
00368               "-fit|fitArcCont 1,1 -tol 1.e-2 -default -blaze -nolocal");
00369 
00370   init_snifs("$Name:  $");
00371   init_session(argv,argc,&arglabel,&argval);
00372 
00373   if (DEBUG) {
00374     print_msg("$Id: create_mask.c,v 1.23 2004/10/06 16:56:33 ycopin Exp $");
00375     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00376   }
00377 
00378   /* ===== Parameters ============================== */
00379 
00380   /* ----- Input ------------------------------ */
00381 
00382   maskname = argval[0]; /* Input mask */
00383   maxname  = argval[1]; /* Input max */
00384   arcname  = argval[2]; /* Input calibration frame */
00385   refname  = argval[3]; /* Input lambda ref. table */
00386 
00387   /* Select emission lamp and continuum functions to fit */
00388   sscanf(argval[4], "%d,%d", &FITARC, &FITCONT);
00389   SKIPFIT = !(FITARC>0 || FITCONT>0);        /* Skip fitting if no fit... */
00390 
00391   get_argval(5, "%lf", &tol);                /* Fitting tolerance */
00392   DEFAULT = is_true(argval[6]);              /* Configuration */
00393   BLAZE   = is_true(argval[7]);              /* Use blaze function */
00394   NOLOCAL = is_true(argval[8]);              /* Don't compute local adjust. */
00395 
00396   /* ----- Output ------------------------------ */
00397 
00398   print_msg("o Input mask: %s",maskname);
00399   print_msg("o Input max: %s",maxname);
00400   print_msg("o Input calibration frame: %s",arcname);
00401   print_msg("o Input wavelength reference table: %s",refname);
00402   if (SKIPFIT) print_msg("o No fit performed");
00403   else {
00404     if (FITARC>0)  print_msg("o Fitting calibration frame");
00405     if (FITCONT>0) print_msg("o Fitting maxima every %d line(s)",FITCONT);
00406   }
00407   print_msg("o Fitting tolerance: %.2f",tol);
00408   if (DEFAULT) print_msg("o Default configuration used");
00409   else         print_msg("o Reading configuration from mask");
00410   if (BLAZE) print_msg("o Using blaze function in fit "
00411                        "(i.e. fitting 0th- and 2nd-orders)");
00412   else       print_msg("o 1st-order fit");
00413   if (BLAZE) {
00414     print_error("BLAZE option not implemented yet");
00415     exit_session(ERR_NOIMPL);
00416   }
00417 
00418   /* ===== Mask ============================== */
00419 
00420   if ((status = open_mask(maskname,&mask,"IO"))) exit_session(status);
00421   handle_select_flag(&mask,'W',NULL);
00422 
00423   /* ===== Configuration ============================== */
00424 
00425   /* ----- Channel ------------------------------ */
00426 
00427   lbdapup = read_lbdaref(&mask);
00428 
00429   /* ----- Configuration ------------------------------ */
00430 
00431   if (DEFAULT) {
00432     Channel channel;
00433 
00434     if (!(channel = read_channel(&mask))) exit_session(channel);
00435     read_default_snifs_config(channel,&config);
00436   }
00437   else if (read_snifs_config_from_mask(&mask,&config))
00438     exit_session(ERR_NODESC);
00439 
00440   if (DEBUG) print_snifs_config(&config);
00441 
00442   optics = config.optics;
00443 
00444   /* ----- Variables ------------------------------ */
00445 
00446   lbda_inf = optics.filter.inf_util;         /* Useful wavelength coverage */
00447   lbda_sup = optics.filter.sup_util;
00448 
00449   /* ===== Mpup coordinates ============================== */
00450 
00451   npup = mask.row;
00452   print_msg("Reading %d mpup coordinates (%.1f A) from mask...",npup,lbdapup);
00453 
00454   /* Real distorted positions */
00455   colno = get_col_ref(&mask,LAB_COL_NO);
00456   colxd = get_col_ref(&mask,LAB_COL_XLD);
00457   colyd = get_col_ref(&mask,LAB_COL_YLD);
00458   if (colno<0 || colxd<0 || colyd<0) {
00459     print_error("Cannot read specific mask columns in %s",maskname);
00460     exit_session(ERR_NOCOL);
00461   }
00462 
00463   nopup =   (int *)malloc(npup*sizeof(int));    /* Mpup number */
00464   xppup = (float *)malloc(npup*sizeof(float));  /* Mpup position [px] */
00465   yppup = (float *)malloc(npup*sizeof(float)); 
00466 
00467   RD_col(&mask,colno,nopup);
00468   RD_col(&mask,colxd,xppup);
00469   RD_col(&mask,colyd,yppup);
00470 
00471   print_msg("   %d (1st-order) mpups",npup);
00472 
00473   /* ===== Maxima ============================== */
00474 
00475   if ((status = open_max(maxname, &maxset))) exit_session(status);
00476 
00477   if (DEBUG) print_warning("No way to check that the max channel is good");
00478   
00479   if (maxset.nb_ycoords <= 1) {
00480     print_error("Maxima frame contains a single line");
00481     exit_session(ERR_BAD_TYPE);
00482   }
00483 
00484   /* ===== Calibration frame ============================== */
00485 
00486   if ((status = open_image(&arcframe,arcname,"I",-1,"arc")))
00487     exit_session(status);
00488 
00489   frame_stats(&arcframe, &mean, &sigma, NULL, NULL);
00490   print_msg("   Mean: %g, Sigma: %g",mean,sigma);
00491 
00492   /* ===== Wavelength references ============================== */
00493 
00494   /* ----- Open ------------------------------ */
00495 
00496   print_msg("Reading wavelength reference table %s...", refname);
00497   if (open_table(&reftab,refname,"I") < 0) {
00498     print_error("Unable to open %s",refname);
00499     exit_session(ERR_OPEN);
00500   }
00501   handle_select_flag(&reftab,'W',NULL);
00502 
00503   /* open LAMBDA and MASK col */
00504   collbda = get_col_ref(&reftab,"LAMBDA");
00505   colmsk  = get_col_ref(&reftab,"MASK");
00506   if (collbda<0 || colmsk<0) {
00507     print_error("Cannot find columns LAMBDA and MASK in table %s",reftab.name);
00508     exit_session(ERR_NOCOL);
00509   }
00510 
00511   /* ----- Read ------------------------------ */
00512 
00513   /* Read *all* the lambda */ 
00514   tmplbda = (float *)malloc(reftab.row*sizeof(float));
00515   iflag   =   (int *)malloc(reftab.row*sizeof(int));
00516   RD_col(&reftab,collbda,tmplbda);
00517   RD_col(&reftab,colmsk,iflag);
00518 
00519   /* Select flagged lambda in the proper range */
00520   lbda = (double *)malloc(reftab.row*sizeof(double));
00521   for (nlbda=0, i=0; i<reftab.row; i++) {
00522     if (iflag[i] &&
00523         tmplbda[i] >= lbda_inf && tmplbda[i] <= lbda_sup) {
00524       lbda[nlbda++] = tmplbda[i];
00525     }
00526     else if (DEBUG && !iflag[i]) {
00527       print_msg("Line %d (F%d, %.1f A) flagged-out from [%.1f,%.1f]",
00528                 i+1,iflag[i],tmplbda[i],lbda_inf,lbda_sup);
00529     }
00530   }
00531 
00532   free(tmplbda); free(iflag);
00533   close_table(&reftab);
00534 
00535   print_msg("-> %d selected lines in range [%.1f-%.1f]",
00536             nlbda,lbda_inf,lbda_sup);
00537 
00538   if (nlbda<2) {
00539     print_error("Number of lines <2, aborting");
00540     exit_session(UNKNOWN);
00541   }
00542   if (DEBUG) 
00543     for (i=0; i<nlbda; i++) print_msg("   Line %d: %.1f A",i+1,lbda[i]);
00544   
00545   /* ##### FIT MASK ############################## */
00546 
00547   if (!SKIPFIT) {
00548 
00549     /* ===== Initialization ============================== */
00550 
00551     /* ----- Global variables ------------------------------ */
00552 
00553     glFITARC  = MAX(FITARC,0);
00554     glFITCONT = MAX(FITCONT,0); 
00555 
00556     gloptics = &optics;
00557     glmaxset = &maxset;
00558     glframe  = &arcframe;
00559   
00560     glxppup = xppup;
00561     glyppup = yppup;  
00562 
00563     glpixsize = config.ccd.pixsize*1e-3; /* Pixel size [mm] */
00564     gllbdapup = lbdapup;
00565 
00566     glnlbda = nlbda;
00567     gllbda  = lbda;
00568   
00569     glnormmax = mean;
00570     print_warning("*** glnormmax set to arc-frame mean value *** TO BE UPDATED!!!");
00571 
00572     //  /* estimate of normalization factor: take 0.1*max to get optimum contrast */
00573     //  glnormval = 0.1*search_max(nvars,xv); 
00574     //  if (DEBUG) {
00575     //    print_msg("Set Normalization value to %g", glnormval);
00576     //  }
00577 
00578     /* ----- Parameters ------------------------------ */
00579 
00580     /* PARAMETERS: 
00581        Collimator: f,E,lc[3]          (0,1,2..4)
00582        Camera:     f,E,lc[3]          (5,6,7..9)
00583        Grism:      A,tiltx,tilty,rotz (10,11,12,13)
00584        Center:     x,y                (14,15)
00585        => Total = 16
00586     */
00587 
00588     npar = 16;
00589   
00590     par   = (double *)malloc(npar*sizeof(double));
00591     epar  = (double *)malloc(npar*sizeof(double));
00592     parlb = (double *)malloc(npar*sizeof(double));
00593     parub = (double *)malloc(npar*sizeof(double));
00594 
00595     optics2mask_param(&optics,par);
00596     set_mask_param_limits(par,parlb,parub);
00597     print_mask_param_in(par,parlb,parub);
00598 
00599     /* Normalization */
00600 
00601     norm_mask_param(optics.channel,par);
00602     norm_mask_param(optics.channel,parlb);
00603     norm_mask_param(optics.channel,parub);
00604 
00605     /* ===== Fit ============================== */
00606 
00607     nllsqfit_noderiv();        /* No derivatives supplied */
00608 
00609     epar[0] = 1;               /* Eastern egg: Print-out: 1, Derivatives: 2 */
00610 
00611     status = nllsqfit_bnd(NULL, NULL, NULL, npup, par, parlb, parub, npar,
00612                           *nllsqfit_mask, NULL, NULL, epar, &objf,&gof);
00613 
00614     /* ===== Conclusion ============================== */
00615 
00616     print_msg("Fit result: status=%d, Obj=%g",status,objf);
00617 
00618     /* ##### MASK UPDATE ############################## */
00619 
00620     /* ===== Fit result ============================== */
00621 
00622     if (!status || 
00623         status == 6 /* Cannot be improved upon, still probably good */
00624         ) {
00625       /* Unnormalization */
00626       unnorm_mask_param(optics.channel,par);
00627       unnorm_mask_param(optics.channel,epar);
00628 
00629       mask_param2optics(par,&optics);
00630       print_mask_param_out(par,epar);
00631 
00632       config.optics = optics;
00633 
00634       sprintf(config.name,"%s (fit %d,%d)",config.name,FITARC,FITCONT);
00635       /* Store configuration, calibration and max files */
00636       if (write_snifs_config_to_mask(&mask,&config) ||
00637           WR_desc(&mask,MSK_ARC,CHAR,lg_name+1,arcname) ||
00638           WR_desc(&mask,MSK_MAX,CHAR,lg_name+1,maxname) ) 
00639         exit_session(ERR_BAD_DESC);
00640     }
00641     else {
00642       print_error("Aborting, mask configuration left untouched");
00643       exit_session(UNKNOWN);
00644     }
00645   } /* !SKIPFIT */
00646   else {
00647     print_warning("No fit performed, mask configuration left untouched");
00648   }
00649 
00650   /* ##### LOCAL ADJUSTMENT ############################## */
00651 
00652   if (!NOLOCAL && FITCONT) {
00653 
00654     maxstep = (maxset.line[maxset.nb_ycoords-1].ycoord - 
00655                maxset.line[0].ycoord) / (maxset.nb_ycoords - 1);
00656     lstep = optics.grism.w_disp*config.ccd.pixsize*1e-3*maxstep * ABS(FITCONT);
00657 
00658     print_msg("Fitting local adjustments to mask:\n"
00659               "   step=%.1f A, D(dx)=%d, D(sig)=%d, maxdx=%.2f px",
00660               lstep,dxnc-1,signc-1,MAXDX);
00661 
00662     /* ===== Prepare mask ============================== */
00663 
00664     WR_desc(&mask,MSK_DXNC, INT,1,&dxnc);    /* Keywords */
00665     WR_desc(&mask,MSK_SIGNC,INT,1,&signc);
00666     for (i=0; i<dxnc; i++) {                 /* DX0,DX1,... columns */
00667       sprintf(tmpname,"DX%d",i);
00668       coldx[i] = create_col(&mask,tmpname,FLOAT,'R',"F9.6","");
00669     }
00670     colddx = create_col(&mask,"DXRMS",FLOAT,'R',"F9.6","px");
00671     for (i=0; i<signc; i++) {                /* SIG0,SIG1,... columns */
00672       sprintf(tmpname,"SIG%d",i);
00673       colsig[i] = create_col(&mask,tmpname,FLOAT,'R',"F9.6","");
00674     }
00675     coldsig = create_col(&mask,"SIGRMS",FLOAT,'R',"F9.6","px");
00676 
00677     if (DEBUG) {                             /* Debug cubes */
00678       nlbda = (lbda_sup - lbda_inf)/lstep + 1;
00679 
00680       disable_erase_flag();
00681       print_msg("DEBUG: Creating datacubes 'dbg_mask_dx' and 'dbg_mask_sig' "
00682                 "(data in signal, model in noise)");
00683       create_tiger_frame(&cubedist,"dbg_mask_dx",-1,0,lstep,FLOAT,
00684                          maskname,"Local distance to Max",NULL);
00685       create_tiger_frame(&cubesig,"dbg_mask_sig",-1,0,lstep,FLOAT,
00686                          maskname,"Max cross-disp. sigma",NULL);
00687       restore_erase_flag();
00688     }
00689 
00690     /* ===== Lens & max association ============================== */
00691 
00692     for (i=0; i<npup; i++) {                 /* Loop over lenses */
00693 
00694       /* ----- Max detection ------------------------------ */
00695 
00696       /* Look for the nearest max along the spectrum */
00697       nmax = pup_get_maxdata(&maxset, xppup[i], yppup[i], lbdapup, 
00698                              &optics, config.ccd.pixsize*1e-3, lstep,
00699                              &plbda, &pdx, &psigma);
00700     
00701       if (nmax<3) {
00702         print_warning("Lens #%d (%.1fx%.1f): cannot find [m]any max "
00703                       "position close enough (%d), skipping",
00704                       nopup[i],xppup[i],yppup[i],nmax);
00705       }
00706       else {
00707         ngpup ++;                            /* Lenses w/ local adjustment */
00708 
00709         /* Local and total RMS */
00710         varbef = gsl_stats_variance_m(pdx,1,nmax,0);
00711         nmaxtot   += nmax;
00712         varbeftot += (nmax-1)*varbef;
00713 
00714         /* ----- Fit dx(lambda) ------------------------------ */
00715 
00716         status = fit_poly_rej_nag_tab(plbda, pdx, NULL, nmax, 
00717                                       dxnc, -LOC_CLDX,LOC_CLDX,LOC_MXDX, 
00718                                       dxcoeff, &niter, &dxerr, &nrej, 
00719                                       lbda_inf, lbda_sup, &ifail);
00720 
00721         if (!status) {
00722 
00723           /* Residuals variance: one has
00724              dxerr^2 = gsl_stats_variance_m(pdx[j] - val_poly_nag(j...),1,nmax,0) */
00725           varaft = SQ(dxerr);
00726           varafttot += (nmax-1)*varaft;
00727 
00728           if (DEBUG) {                         /* Storage */
00729             init_new_tiger_spec(&cubedist,&specS,nlbda,lbda_inf);
00730             init_new_tiger_spec(&cubedist,&specN,nlbda,lbda_inf);
00731             for (k=0,j=0; j<nlbda; j++) {
00732               WR_spec(&specS,j,                /* Data */
00733                       (k<nmax && (ABS(plbda[k])-coord_spec(&specS,j))<1e-3)? pdx[k++]:0);
00734               WR_spec(&specN,j,                /* Model */
00735                       val_poly_nag(coord_spec(&specS,j), dxcoeff, dxnc, 
00736                                    lbda_inf, lbda_sup, &ifail));
00737             }
00738             put_tiger_spec(&cubedist,&specS,&specN,nopup[i]);
00739           }
00740         }
00741 
00742         if (DEBUG && (status || nrej)) {
00743           print_warning("Lens #%d (%.1fx%.1f): "
00744                         "dx(lbda): status=%d, nrej=%d, niter=%d",
00745                         nopup[i],xppup[i],yppup[i],status,nrej,niter);
00746         }
00747 
00748         /* ----- Fit sigma(lambda) ------------------------------ */
00749 
00750         status = fit_poly_rej_nag_tab(plbda, psigma, NULL, nmax, 
00751                                       signc, -LOC_CLSIG,LOC_CLSIG,LOC_MXSIG, 
00752                                       sigcoeff, &niter, &sigerr, &nrej, 
00753                                       lbda_inf, lbda_sup, &ifail);
00754 
00755         if (!status && DEBUG) {              /* Storage */
00756           init_new_tiger_spec(&cubesig,&specS,nlbda,lbda_inf);
00757           init_new_tiger_spec(&cubesig,&specN,nlbda,lbda_inf);
00758           for (k=0, j=0; j<nlbda; j++) {
00759             WR_spec(&specS,j,                /* Data */
00760                     (k<nmax && (ABS(plbda[k])-coord_spec(&specS,j))<1e-3)? psigma[k++]:0);
00761             WR_spec(&specN,j,                /* Model */
00762                     val_poly_nag(coord_spec(&specS,j), sigcoeff, signc, 
00763                                  lbda_inf, lbda_sup, &ifail));
00764           }
00765           put_tiger_spec(&cubesig,&specS,&specN,nopup[i]);
00766         }
00767       
00768         if (DEBUG && (status || nrej)) {
00769           print_warning("Lens #%d (%.1fx%.1f): "
00770                         "sig(lbda): status=%d, nrej=%d, niter=%d",
00771                         nopup[i],xppup[i],yppup[i],status,nrej,niter);
00772         }
00773 
00774         if (DEBUG) {                         /* Summary print-out */
00775           print_msg("Lens #%3d (%4.0fx%4.0f): RMS(dx)=%.2f->%.3f px, "
00776                     "RMS(sig)=%.2f px",
00777                     nopup[i],xppup[i],yppup[i],
00778                     sqrt(varbef),sqrt(varaft),sigerr);
00779         }
00780 
00781         /* ===== Writing local adjustments in mask ========================= */
00782 
00783         for (j=0; j<dxnc; j++) {
00784           fval = dxcoeff[j]; 
00785           WR_tbl(&mask,i,coldx[j],&fval);
00786         }
00787         fval = dxerr; 
00788         WR_tbl(&mask,i,colddx,&fval);
00789         for (j=0; j<signc; j++) {
00790           fval = sigcoeff[j]; 
00791           WR_tbl(&mask,i,colsig[j],&fval);
00792         }
00793         fval = sigerr; 
00794         WR_tbl(&mask,i,coldsig,&fval);
00795 
00796       } /* Enough max to compute local adjustments */
00797 
00798       if (nmax) {                           /* Allocated in pup_get_maxdata */
00799         free(plbda); free(pdx); free(psigma); 
00800       }
00801 
00802     } /* Loop over lenses */
00803 
00804     if (DEBUG) {                             /* Close debug cubes */
00805       write_file_class(&cubesig,DEBUG_FILE);  close_tiger_frame(&cubesig);
00806       write_file_class(&cubedist,DEBUG_FILE); close_tiger_frame(&cubedist);
00807     }
00808 
00809     /* BEWARE: RMS(after) is wrong if some dx fits failed (nmaxtot is then
00810        updated, but not varafttot). */
00811     print_msg("=> %d/%d mpups w/ %d associated max, RMS(dx)=%.2f->%.4f px",
00812               ngpup,npup,nmaxtot,
00813               sqrt(varbeftot/(nmaxtot-1)),sqrt(varafttot/(nmaxtot-1)));
00814 
00815   } /* Local adjustments */
00816   
00817   /* ##### Conclusions ############################## */
00818 
00819   free(nopup); free(xppup); free(yppup); free(lbda);
00820   
00821   free(par); free(epar); free(parlb); free(parub);
00822 
00823   close_table(&mask);
00824   close_frame(&arcframe);
00825 
00826   exit_session(OK);
00827   return(OK);
00828 }
00829 
00830 /* 
00831 ;;; Local Variables: ***
00832 ;;; eval: (add-to-list 'c-font-lock-extra-types "SnifsConfig") ***
00833 ;;; eval: (add-to-list 'c-font-lock-extra-types "SnifsOptics") ***
00834 ;;; End: ***
00835 */

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