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

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