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

extract/source/extract_spec.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 //#include <stdbool.h> /* C99 Booleans (gcc 3.0) */
00021 
00022 #define DBG_PREFIX "dbg_ext"    
00023 
00024 static int const ARCWIDTH=128; 
00025 
00026 static int const ORDER=1;       
00027 
00028 static int const XOFY_DEGMAX=7;  
00029 static int const YOFL_DEGMAX=7;  
00030 static int const LOFY_DEGMAX=7;  
00031 static int const POLY_NPTS=50;   
00032 
00033 #define SMAX 2.0                
00034 #define SMIN 0.1                
00035 
00036 static int const DBGSIDEWIN=3;   
00037 static int const XTRALENGTH=30;  
00038 
00039 /* Window range to look for neighbors */
00040 #define NEIGHBRANGE_C 0.5  
00041 /* 1.3 (< 1.5) avoids to find a 3rd neighb.*/
00042 #define NEIGHBRANGE_F 1.3  
00043 /* No need to look for neighbors further than TOOFAR*interspec */
00044 #define TOOFAR 3           
00045 
00046 long *glnpeaks;
00047 double *glpeaks;
00048 double ***glScs, ***glIcs;
00049 
00050 /* Global variables describing the geometric PSF */
00051 int glnggeo;                                 /* Nb of gaussian components */
00052 double *glIgeo, *glXgeo, *glSgeo;
00053 
00054 /* Global variables describing the global PSF */
00055 int glngpsf;                                 /* Nb of gaussian components */
00056 
00057 /* Global variables describing the core PSF = geom. PSF (x) global PSF */
00058 double **glIpsf, **glSpsf;           /* Global variables set by set_corePSF */
00059 
00060 /* Global variables describing the local PSF */
00061 int glngloc;                                 /* Nb of gaussian components */
00062 
00063 /* === Doxygen Comment ======================================= */
00079 /* =========================================================== */
00080 
00081 long funct1(long *n, double *I, double *result) 
00082 { 
00083   int no, i, j, np;
00084   float dx;
00085   double model;
00086 
00087   *result = 0;
00088   np = *n;
00089 
00090   for (i=0; i<glnpeaks[0]; i++) {            /* Loop over the pixels */
00091     model = 0;
00092     for (j=0; j<np; j++) {                   /* Loop over the peaks */
00093       no = (int)glnpeaks[j+1];               /* Lens number */
00094       /* Offset between current position and peak position */
00095       dx = (float)glpeaks[i] - glpeaks[2*glnpeaks[0] + j];
00096       model += I[j]*Mpup_PSF(no,dx);
00097     }
00098 
00099     *result += SQ(glpeaks[glnpeaks[0]+i] - model);
00100   }
00101   return(OK);
00102 }
00103 
00104 /* === Doxygen Comment ======================================= */
00155 /* =========================================================== */
00156 
00157 int main(int argc, char **argv)
00158 {
00159   Channel channel;
00160 
00161   SnifsConfig config;
00162   SnifsOptics optics;
00163 
00164   TABLE mask, tab, dbg_tab1, dbg_tab2, dbg_tab3;
00165   IMAGE2D inframe, invarframe;
00166   TIGERfile outcube, dbgcube;
00167   SPECTRUM spectrum, noise, dbg_spec, dbg_noise;
00168   
00169   char **argval, **arglabel;
00170   char *maskname, *inname, *outname, *arcname, *arc2name;
00171   char refname[lg_name+1],dbgname[lg_name+1];
00172   char ident[lg_ident+1];
00173   
00174   int iq[3], *yofl_deg, *lofy_deg, *xofy_deg, *nol, *I, *J; 
00175   int i, l, k, n, p;
00176   int iima, jima, ilef, irig, jspec, l0, ln, lc, ldbg, ltab; /*di*/
00177   int ibegdbg, ienddbg, iinf, icen, isup;
00178   int status, max_npix, nlens, npix;
00179   int colno, coli, colj, colxdl, colydl, colxnd, colynd;
00180   int cola, cold, colao, coldo, fclass;
00181   int coly1dbg, colxydbg, colxfdbg, coldxdbg;
00182   int collldbg, colyldbg, colyfdbg, coldydbg;
00183   int coly2dbg, collydbg, collfdbg, coldldbg;
00184   int colnodbg, colxdbg, colydbg, colldbg;
00185   int signc, dxnc;
00186   int flag_out, ld_min, lg_min, no_debug, l_debug, dbgno;
00187   int lneighbors[2], rneighbors[2], nlneighbors, nrneighbors;
00188   int nggeo,ngpsf,ngloc;
00189   int NOISE, WEIGHT, OVERLAP, LENS, LOCAL, LINFIT, MEGA_VERBOSE, OFFSET, CAL2, VARIANCE;
00190   /* C99 Booleans (gcc 3.0) */
00191   //  bool NOISE, WEIGHT, OVERLAP, LENS, LOCAL, LINFIT, MEGA_VERBOSE, OFFSET, CAL2, VARIANCE;
00192   
00193   int ibeg, iend;
00194   long *iw, nvar, ibound, npixfit, liw, lw, ifail;
00195   double *bl, *bu, *w, *q, res;
00196   long e04jaf(), funct1();
00197   
00198   float **xofy_coeffs, **yofl_coeffs, **lofy_coeffs;
00199   float *lambda, *lbda_min, *lbda_max, *local_lbda_inf_total, *local_lbda_sup_total;
00200   float *xdl, *ydl;
00201   float *xlo, *ylo, *xup, *yup, *xlot, *ylot, *xupt, *yupt;
00202   float *signal, *variance;
00203   float fval[3], Dx[POLY_NPTS], xd[POLY_NPTS], yd[POLY_NPTS];
00204   float xlef, xrig, alpha, delta;
00205   float xinf, xsup, xref, yref, x, y, flef, frig;
00206   float xi, rms, lbda, interspec, toofar;
00207   float val, maxval, minval;
00208   float ref_airmass, cal_airmass=1, cal2_airmass=1, obj_airmass;
00209   float obj_jdate, cal_jdate, cal2_jdate; 
00210   float half_width, sum, varsum, lfrac, rfrac, fk, dx;
00211   float lbda_inf, lbda_sup, lstep, lbda_inf_total, lbda_sup_total;
00212   float lbdaref;
00213 
00214   double **sigcoeff, **dxcoeff;
00215   double *psf_param, *geo_param;
00216   double qloc[3], offset[2], dval[4], lbdasig[3], sigma[3];
00217   double xmla,ymla, xccd,yccd;
00218   double pinf, pcen, psup, vi, s_vi, s_pcen;
00219   double xc, yc, ftrans;
00220   double outstep, pixsize, RoN2;
00221   double lens_size_pix, lens_size_arcsec;
00222   double qinf=0., qcen=0., qsup=0.;
00223 
00224   /* ##### INTRODUCTION ############################## */
00225   
00226   set_purpose("Spectra extraction (CCD to datacube)");
00227   set_arglist("-mask mask -in none -out none -width 5 -cal none "
00228               "-weighted -lens -1 -overlap -shift|offset NULL "
00229               "-cal2 NULL -local -linfit -verbose -variance");
00230 
00231   init_snifs("$Name:  $");
00232   init_session(argv,argc,&arglabel,&argval);
00233   
00234   if (DEBUG) {
00235     print_msg("$Id: extract_spec.c,v 1.35 2004/11/18 23:03:32 ycopin Exp $");
00236     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00237   }
00238 
00239   /* ===== Parameters ============================== */
00240   
00241   /* ----- Input ------------------------------ */
00242   
00243   maskname = argval[0];                           /* Input mask */
00244   inname   = argval[1];                          /* Input frame */
00245   outname  = argval[2];                      /* Output datacube */
00246   get_argval(3,"%f",&half_width);    /* Extraction window width */
00247   arcname  = argval[4];                          /* Input frame */
00248   WEIGHT   = is_true(argval[5]);          /* Optimal extraction */
00249   get_argval(6,"%d",&no_debug);       /* Single lens extraction */
00250   LENS     = (no_debug>0);
00251   OVERLAP  = is_true(argval[7]);            /* Handle pollution */
00252   if ((OFFSET = is_set(argval[8])) &&    /* User-defined offset */
00253       sscanf(argval[8], "%lf,%lf",&offset[0],&offset[1]) != 2) { 
00254     print_error("Invalid option -offset dx,dy (has '%s')", argval[8]);
00255     exit_session(ERR_BAD_PARAM);
00256   }
00257   if ((CAL2 = is_set(argval[9])))
00258     arc2name = argval[9];              /* 2nd calibration frame */
00259   LOCAL  = is_true(argval[10]);             /* Local adjustment */
00260   LINFIT = is_true(argval[11]);    /* Fit of the X-disp profile */
00261   MEGA_VERBOSE = is_true(argval[12]);           /* Mega-verbose */
00262   VARIANCE = is_true(argval[13]);     /* Use variance extension */
00263 
00264   /* ----- Output ------------------------------ */
00265 
00266   print_msg("o Input mask: %s",maskname);
00267   print_msg("o Input frame: %s",inname);
00268   print_msg("o Output datacube: %s",outname);
00269   if (LENS) print_msg("o Single lens extraction: #%d",no_debug);
00270   print_msg("o Extraction window full-width: %.1f",half_width);
00271   print_msg("o Calibration frame: %s",arcname);
00272   if (CAL2) print_msg("o 2nd calibration frame: %s",arc2name);
00273   if (OFFSET) print_msg("o User-defined offset for calibration frame: "
00274                         "%.2f x %.2f px", offset[0],offset[1]);
00275   if (WEIGHT) print_msg("o Extraction using Weighted summation");
00276   if (OVERLAP) print_msg("o Extraction using Inter-spectrum depollution");
00277   if (LOCAL) print_msg("o Local adjustment to mask positions");
00278   if (LINFIT) print_msg("o Extraction using Linear fit");
00279 
00280   if (LINFIT) {
00281     print_warning("Option -linfit not fully tested yet, use with caution!");
00282   }
00283 
00284   if (VARIANCE) print_msg("o Using [signal] and [variance] extensions of input frame");
00285   print_msg("");
00286 
00287   /* ----- Misc. ------------------------------ */
00288 
00289   if (DEBUG) set_control_level(WARNING);
00290   else       set_control_level(FATAL);
00291 
00292   half_width /= 2;
00293 
00294   /* ===== Mask ============================== */
00295 
00296   if ((status = open_mask(maskname,&mask,"I"))) exit_session(status);
00297   handle_select_flag(&mask,'W',NULL);
00298  
00299   /* ----- Channel ------------------------------ */
00300 
00301   lbdaref = read_lbdaref(&mask);
00302 
00303   /* ----- Configuration ------------------------------ */
00304 
00305   if (read_snifs_config_from_mask(&mask,&config)) 
00306     exit_session(ERR_NODESC);
00307 
00308   if (DEBUG) print_snifs_config(&config);
00309 
00310   optics = config.optics;
00311 
00312   /* ----- Variables ------------------------------ */
00313 
00314   xc       = optics.center.xc;
00315   yc       = optics.center.yc;
00316   lbda_inf = optics.filter.inf_util;
00317   lbda_sup = optics.filter.sup_util;
00318   if (WEIGHT) {
00319     lbda_inf_total = optics.filter.inf_total;
00320     lbda_sup_total = optics.filter.sup_total;
00321   } 
00322 
00323   pixsize = config.ccd.pixsize*1e-3;         /* Pixel size [mm] */
00324 
00325   lens_size_pix    = optics.MLA.sizepix;
00326   lens_size_arcsec = optics.MLA.size / (config.telescope.scale*1e-3*optics.enlarger.gamma); 
00327   ftrans           = lens_size_arcsec / lens_size_pix;
00328 
00329   /* ===== PSF models ============================== */
00330 
00331   if (WEIGHT) {
00332 
00333     print_msg("Reading PSF models from mask descriptors");
00334 
00335     /* ----- Geometric PSF ------------------------------ */
00336 
00337     if ((nggeo = read_PSFgeo_param(&mask, &geo_param)) < 0)
00338         exit_session(ERR_NODESC);
00339 
00340     glnggeo = nggeo;
00341     glIgeo = (double *)malloc(nggeo*sizeof(double));
00342     glXgeo = (double *)malloc(nggeo*sizeof(double));
00343     glSgeo = (double *)malloc(nggeo*sizeof(double));  
00344     /* GEOPAR already include the scale parameter, and is properly normed */
00345     for (i=0; i<nggeo; i++) {
00346       glIgeo[i] = geo_param[3*i  ];
00347       glXgeo[i] = geo_param[3*i+1];
00348       glSgeo[i] = geo_param[3*i+2];
00349       if (DEBUG) print_msg("\tGaussian %d: I = %f, X = %f, S = %f",i+1,
00350                            geo_param[3*i],geo_param[3*i+1],geo_param[3*i+2]);
00351     }
00352     free(geo_param);
00353 
00354     /* ----- Global PSF ------------------------------ */
00355 
00356     if ((ngpsf = read_PSFglobal_param(&mask, &psf_param)) < 0)
00357       exit_session(ERR_NODESC);
00358 
00359     if (DEBUG) 
00360       for (i=0; i<ngpsf; i++) 
00361         print_msg("\tUnconvolved gaussian %d: I = %f, S = %f", 
00362                   i+1,psf_param[2*i+1],psf_param[2*i+2]);
00363 
00364     /* ----- Local PSF ------------------------------ */
00365 
00366     if ((ngloc = read_PSFlocal_ng(&mask)) == MAXINT)
00367       exit_session(ERR_NODESC);
00368 
00369     /* ----- Core PSF ------------------------------ */
00370 
00371     if (ngloc > 0) {
00372       print_msg("   Core PSF = Geometric PSF (x) Global PSF");
00373     }
00374     else {
00375       ngloc = -ngloc;
00376       ngpsf = 0;
00377       print_msg("   Core PSF = Geometric PSF (global PSF *discarded*)");
00378     }
00379 
00380     if (ngloc > 1) {
00381       print_error("Multi-gaussian local PSF model not implemented yet");
00382       exit_session(ERR_NOIMPL);
00383     }
00384 
00385     glngpsf = ngpsf;
00386     alloc2d(&glIpsf,nggeo,MAX(ngpsf,1),DOUBLE);
00387     alloc2d(&glSpsf,nggeo,MAX(ngpsf,1),DOUBLE);
00388 
00389     /* If ngpsf=0, psfpar is not needed, since core PSF = geom. PSF then */
00390     set_corePSF(psf_param);
00391 
00392     free(psf_param);
00393 
00394   } /* WEIGHT */
00395   
00396   /* ===== Mpup coordinates ============================== */
00397 
00398   print_msg("Reading mpup coordinates (%.1f A) from mask...",lbdaref);
00399 
00400   nlens = mask.row;
00401 
00402   nol  = (int *)  malloc(nlens*sizeof(int));
00403   /* DEF: xdl,ydl = real distorted positions */
00404   xdl  = (float *)malloc(nlens*sizeof(float));
00405   ydl  = (float *)malloc(nlens*sizeof(float));
00406   /* DEF: I,J = position indices */
00407   I    = (int *)  malloc(nlens*sizeof(int));
00408   J    = (int *)  malloc(nlens*sizeof(int));
00409 
00410   /* C99 Variable Length Arrays (unimplemented in gcc 3.4) */
00411   // int *J[nlens];                            /* Not reallocable */
00412 
00413   colno   = get_col_ref(&mask,LAB_COL_NO);
00414   colxdl  = get_col_ref(&mask,LAB_COL_XLD);
00415   colydl  = get_col_ref(&mask,LAB_COL_YLD);
00416   coli    = get_col_ref(&mask,LAB_COL_I);
00417   colj    = get_col_ref(&mask,LAB_COL_J);
00418   if ((colno<0) || (colxdl<0) || (colydl<0) || (coli<0) || (colj<0)) {
00419     print_error("Unable to read specific columns from mask");
00420     exit_session(ERR_BAD_COL);
00421   }
00422 
00423   RD_col(&mask,colno,   nol);
00424   RD_col(&mask,colxdl,  xdl);
00425   RD_col(&mask,colydl,  ydl);
00426   RD_col(&mask,coli,      I);
00427   RD_col(&mask,colj,      J);
00428 
00429   /* ----- Single lens ------------------------------ */
00430 
00431   if (LENS) {
00432     for (l=0; l<nlens && nol[l]!=no_debug; l++);
00433     if (l < nlens) {
00434       print_msg("Debug lens #%d (%7.2f x %7.2f) is at line %d", 
00435                 nol[l], xdl[l], ydl[l], l+1);
00436       l_debug = l;
00437     } 
00438     else {
00439       print_error("Debug lens #%d cannot be found", no_debug);
00440       exit_session(ERR_BAD_PARAM);
00441     }
00442   }
00443 
00444   /* ===== Weighted summation ============================== */
00445 
00446   if (WEIGHT) {
00447 
00448     /* Reading sigma(lambda) from max analysis in create_mask */
00449     if ((signc = read_local_model(&mask,2,&sigcoeff)) < 0)
00450       exit_session(ERR_NODESC);
00451 
00452     /* Computing cross-dispersion profile model */
00453     
00454     glScs = (double ***)malloc(nlens*sizeof(double **));
00455     glIcs = (double ***)malloc(nlens*sizeof(double **));
00456     for (l=0; l<nlens; l++) {                /* loop over the lenses */
00457       glScs[l] = (double **)malloc(glnggeo*sizeof(double *));
00458       glIcs[l] = (double **)malloc(glnggeo*sizeof(double *));
00459       for (i=0; i<glnggeo; i++) {      /* loop over the geometric gaussians */
00460         glScs[l][i] = (double *)malloc(glngpsf*sizeof(double));
00461         glIcs[l][i] = (double *)malloc(glngpsf*sizeof(double));
00462       }
00463     }
00464 
00465     /* Local spectral range */
00466 
00467     /* Use the model you want to set the local bandpass if needed. For Sauron,
00468        we use an ad-hoc model function of I,J, but could be stored in mask if
00469        you think it's better or more handy... */
00470 
00471     /* DEF: local_lbda_inf|sup_total: local (ie lens specific) total spectral range */
00472     local_lbda_inf_total = (float *)malloc(nlens*sizeof(float));
00473     local_lbda_sup_total = (float *)malloc(nlens*sizeof(float));    
00474     for (l=0; l<nlens; l++) {
00475       local_lbda_inf_total[l] = lbda_inf_total;
00476       local_lbda_sup_total[l] = lbda_sup_total;
00477     }
00478     
00479     if (LENS) {
00480       print_msg("Total spectral range for lens #%d: %7.2f -- %7.2f AA", 
00481                 nol[l_debug], 
00482                 local_lbda_inf_total[l_debug], local_lbda_sup_total[l_debug]);
00483     }
00484 
00485   } /* WEIGHT */
00486 
00487   /* ===== Local adjustment ============================== */
00488 
00489   if (LOCAL) {
00490 
00491     /* Reading dx(lambda) from max analysis in create_mask */
00492     if ((dxnc = read_local_model(&mask,1,&dxcoeff)) < 0)
00493       exit_session(ERR_NODESC);
00494 
00495   } /* LOCAL */
00496 
00497   /* ===== Calibration Frames ============================== */
00498 
00499   if (!OFFSET) {
00500     
00501     /* ----- Mask reference frame ------------------------------ */
00502 
00503     if (RD_desc(&mask,MSK_ARC,CHAR,lg_name+1, refname) <= 0) {
00504       print_error("Cannot read name of reference calibration frame in "
00505                   "descriptor %s of mask %s",MSK_ARC,maskname);
00506       exit_session(ERR_BAD_DESC);
00507     }
00508 
00509 //    /* Expanding the descriptor to full pathname: By default, the reference
00510 //       calibration frame is located along the mask */
00511 //    if (strcmp(mask.path,"") != 0 && 
00512 //        strcmp(mask.path,"./") != 0) 
00513 //      sprintf(refname, "%s/%s", mask.path, refname);
00514     print_warning("Reference calibration frame %s supposed to be in "
00515                   "current directory (TABLE.path)",refname);
00516 
00517     if ((status = comp_calibration_offsets(refname, &ref_airmass, 
00518                                            arcname, &cal_airmass, &cal_jdate,
00519                                            (CAL2? arc2name:NULL), 
00520                                            &cal2_airmass, &cal2_jdate,
00521                                            offset)))
00522       exit_session(status);
00523 
00524     if (CAL2)
00525       print_msg("Mean calibration offset: %.2f x %.2f px", offset[0],offset[1]);
00526     else 
00527       print_msg("Calibration offset: %.2f x %.2f px",offset[0],offset[1]);
00528 
00529   } /* !OFFSET */
00530   else {
00531     print_msg("Calibration offset (manual): %.2f x %.2f px",offset[0],offset[1]);
00532   }
00533   
00534   close_table(&mask);
00535 
00536   /* ===== Science frame ============================== */
00537 
00538   /* Open frame & leave file class check for latter */
00539   //  if ((status = open_image(&inframe,inname,"I",0,"input")))
00540   //    exit_session(status);
00541   if ((status = open_image_ext(&inframe,(VARIANCE? &invarframe:NULL),
00542                                inname,"I",0,"input")))
00543     exit_session(status);
00544 
00545   /* Check channel */
00546   if ((channel = read_channel(&inframe)) != optics.channel) {
00547     print_error("mask channel (%s) and input frame channel (%s) are different",
00548                 channel_names[optics.channel],channel_names[channel]);
00549     exit_session(ERR_BAD_PARAM);
00550   }
00551 
00552   /* Check file class */
00553   switch ((fclass = read_file_class(&inframe))) {
00554   case PRE_OBJ_FRAME:
00555   case PRE_SKY_FRAME:
00556   case PRE_DOM_FRAME: 
00557   case PRE_CON_FRAME:
00558   case FLT_OBJ_FRAME:
00559   case FLT_SKY_FRAME:
00560   case FLT_DOM_FRAME:
00561   case FLT_CON_FRAME: NOISE = 1; break;
00562   case SUPER_CLASS:
00563   case PRE_CAL_FRAME:
00564   case FLT_CAL_FRAME: NOISE = 0; break;
00565   default:
00566     print_error("Wrong file class %d: %s should be a "
00567                 "`Preproc./HFFF CCD exposure'",fclass,inframe.name);
00568     exit_session(ERR_BAD_PARAM);
00569   }
00570 
00571   /* ----- Check airmass and julian date ------------------------------ */
00572 
00573   obj_airmass = read_airmass(&inframe); /* Read AIRMASS descriptor */
00574 
00575   check_airmass(arcname, cal_airmass, inname, obj_airmass);
00576   if (CAL2) check_airmass(arc2name, cal2_airmass, inname, obj_airmass);
00577 
00578   obj_jdate = read_juliandate(&inframe); /* Read MJDATE descriptor */
00579 
00580   check_jdate(arcname, cal_jdate, inname, obj_jdate);
00581   if (CAL2) check_jdate(arc2name, cal2_jdate, inname, obj_jdate);
00582 
00583   /* ----- Other descriptors ------------------------------ */
00584  
00585   disable_user_warnings();
00586 
00587   /* CCD gain is not used anymore, as preprocessed frames are already in e-
00588      (and not anymore in ADU) */
00589 
00590   if (NOISE || WEIGHT) {             /* RoN is also needed for opt. extract */
00591     RoN2 = read_mean_RoN(&inframe, optics.channel);
00592     RoN2 *= RoN2;
00593   }
00594 
00595   if (RD_desc(&inframe,OBJECT,CHAR,lg_ident+1,ident) <= 0) {
00596     print_warning("Cannot read object name in descriptor %s of file %s",
00597                   OBJECT,inname);
00598 
00599     sprintf(ident,"Frame %s",inname);
00600     print_warning("%s descriptor defaulted to '%s'",OBJECT,ident);
00601   }
00602   if (DEBUG) print_msg("   Object: '%s'", ident);
00603 
00604   restore_user_warnings();
00605 
00606   /* ##### EXTRACTION PREPARATION ############################## */
00607 
00608   /* ===== Lambda step ============================== */
00609 
00610   /* DEF: xlo,ylo = lower (red) useful end (on CCD) of spectrum */
00611   xlo = (float *)malloc(nlens*sizeof(float));
00612   ylo = (float *)malloc(nlens*sizeof(float));
00613   /* DEF: xup,yup = upper (blue) useful end (on CCD) of spectrum */
00614   xup = (float *)malloc(nlens*sizeof(float));
00615   yup = (float *)malloc(nlens*sizeof(float));
00616 
00617   max_npix = 0; /* max number of points in the wavelength range for all lens */
00618   for (l=0; l<nlens; l++) {
00619 
00620     /* [CCD/px] -> [MLA/mm] */
00621     snifs_optics_CCD2MLA(&optics, xdl[l],ydl[l], &xmla,&ymla, 
00622                          pixsize,lbdaref);
00623 
00624     /* [MLA/mm] -> [CCD/px]: Useful spectral extension */
00625     snifs_optics_MLA2CCD(&optics, xmla,ymla, &xccd,&yccd, 
00626                          pixsize,lbda_sup,ORDER);
00627     xlo[l] = xccd; ylo[l] = yccd;
00628     
00629     snifs_optics_MLA2CCD(&optics, xmla,ymla, &xccd,&yccd,
00630                          pixsize,lbda_inf,ORDER);
00631     xup[l] = xccd; yup[l] = yccd;
00632 
00633     max_npix = MAX(max_npix,  /* spectral extension [px] (cf. def. of npix) */
00634                    rint(ABS(yup[l]-ylo[l])/inframe.stepy));
00635   }
00636 
00637   outstep = (lbda_sup - lbda_inf)/(max_npix - 1); /* step [AA/px] */
00638   outstep = floor(outstep*1e2)*1e-2;         /* 2-digits truncation */
00639 
00640   /* WARNING: npix is not directly comparable to max_npix because of ylo,yup
00641      which are updated several times in tiny ways. Take some precautions... */
00642   max_npix = (lbda_sup-lbda_inf)/outstep + 1 + 4; 
00643 
00644   /* ===== Ouptput datacube creation ============================== */
00645 
00646   print_msg("Creating output datacube '%s'...", outname);
00647   if (LENS) {
00648     print_warning("Single lens extraction will overwrite datacube '%s'",
00649                   outname);
00650   }
00651 
00652   create_tiger_frame(&outcube,outname,-1,lbda_inf,outstep,FLOAT,
00653                      outname,ident,"e-");
00654   print_msg("   Wavelength start: %.2f AA, step: %.2f AA",lbda_inf,outstep);
00655 
00656   print_msg("Creating associated lens table...");  /* Associated lens table */
00657   create_table(&tab,outname,-1,9,'Q',ident);
00658 
00659   colno = create_col(&tab,LAB_COL_NO,LONG,'N',"I4", "Lens Number");
00660   coli  = create_col(&tab,"I",     LONG, 'N',"I4",  "vector");
00661   colj  = create_col(&tab,"J",     LONG, 'N',"I4",  "vector");
00662   colxnd= create_col(&tab,"X0",    FLOAT,'N',"F6.2","px");
00663   colynd= create_col(&tab,"Y0",    FLOAT,'N',"F6.2","px");
00664   cola  = create_col(&tab,"A",     FLOAT,'N',"F6.3","arcsec");
00665   cold  = create_col(&tab,"D",     FLOAT,'N',"F6.3","arcsec");
00666   colao = create_col(&tab,"A_ORIG",FLOAT,'N',"F6.3","arcsec");
00667   coldo = create_col(&tab,"D_ORIG",FLOAT,'N',"F6.3","arcsec");
00668 
00669   /* DEF: xofy_deg = polynomial degree (<XOFY_DEGMAX) for xd   = f(yd) */
00670   /* DEF: yofl_deg = polynomial degree (<YOFL_DEGMAX) for yd   = f(lambda) */
00671   /* DEF: lofy_deg = polynomial degree (<LOFY_DEGMAX) for lbda = f(yd) */
00672   xofy_deg = (int *)malloc(nlens*sizeof(int));
00673   yofl_deg = (int *)malloc(nlens*sizeof(int));
00674   lofy_deg = (int *)malloc(nlens*sizeof(int));
00675 
00676   /* DEF: xofy_coeffs = polynomial coeffs of xd   = f(yd) */
00677   /* DEF: yofl_coeffs = polynomial coeffs of yd   = f(lambda) */
00678   /* DEF: lofy_coeffs = polynomial coeffs of lbda = f(yd) */
00679   alloc2d(&xofy_coeffs,nlens,3*(XOFY_DEGMAX+1),FLOAT);
00680   alloc2d(&yofl_coeffs,nlens,3*(YOFL_DEGMAX+1),FLOAT);
00681   alloc2d(&lofy_coeffs,nlens,3*(LOFY_DEGMAX+1),FLOAT);
00682 
00683   if (WEIGHT) {
00684     /* DEF: xlot,ylot,xupt,yupt = lower/upper total end (on CCD) of spectrum 
00685        now taking care of the local bandpass of the filter */
00686     xlot = malloc(nlens*sizeof(float));
00687     ylot = malloc(nlens*sizeof(float));
00688     xupt = malloc(nlens*sizeof(float));
00689     yupt = malloc(nlens*sizeof(float));
00690   }
00691 
00692   if (LENS) {
00693     print_msg("DEBUG: Creating tables %s_fit_xofy, _yofl and _lofy...",
00694               DBG_PREFIX);
00695 
00696     disable_erase_flag();
00697     sprintf(dbgname, "%s_fit_xofy", DBG_PREFIX);
00698     create_table(&dbg_tab1,dbgname,-1,4,'Q',"Polynomial fit: x=x(y)");
00699     coly1dbg = create_col(&dbg_tab1,"YD",FLOAT,'N',"F9.6","");
00700     colxydbg = create_col(&dbg_tab1,"XD",FLOAT,'N',"F9.6","");
00701     colxfdbg = create_col(&dbg_tab1,"XF",FLOAT,'N',"F9.6","");
00702     coldxdbg = create_col(&dbg_tab1,"DX",FLOAT,'N',"F9.6","");
00703 
00704     sprintf(dbgname, "%s_fit_yofl", DBG_PREFIX);
00705     create_table(&dbg_tab2,dbgname,-1,4,'Q',"Polynomial fit: y=y(lambda)");
00706     collldbg = create_col(&dbg_tab2,"LA",FLOAT,'N',"F9.6","");
00707     colyldbg = create_col(&dbg_tab2,"YD",FLOAT,'N',"F9.6","");
00708     colyfdbg = create_col(&dbg_tab2,"YF",FLOAT,'N',"F9.6","");
00709     coldydbg = create_col(&dbg_tab2,"DY",FLOAT,'N',"F9.6","");
00710 
00711     sprintf(dbgname, "%s_fit_lofy", DBG_PREFIX);
00712     create_table(&dbg_tab3,dbgname,-1,4,'Q',"Polynomial fit: lbda=lbda(y)");
00713     coly2dbg = create_col(&dbg_tab3,"YD",FLOAT,'N',"F9.6","");
00714     collydbg = create_col(&dbg_tab3,"LA",FLOAT,'N',"F9.6","");
00715     collfdbg = create_col(&dbg_tab3,"LF",FLOAT,'N',"F9.6","");
00716     coldldbg = create_col(&dbg_tab3,"DL",FLOAT,'N',"F9.6","");
00717 
00718     restore_erase_flag();
00719   }
00720 
00721   /* ===== Wavelength sampling ============================== */
00722 
00723   lstep = (lbda_sup - lbda_inf)/(POLY_NPTS-1);
00724 
00725   lambda = (float *)malloc(POLY_NPTS*sizeof(float));
00726   for (k=0; k<POLY_NPTS; k++) lambda[k] = lbda_inf + k*lstep;
00727 
00728   /* ===== Polynomial fits ============================== */
00729 
00730   /* DEF: lbda_min,lbda_max = wavelength bounds of spectrum */
00731   lbda_min  = (float *)malloc(nlens*sizeof(float));
00732   lbda_max  = (float *)malloc(nlens*sizeof(float));
00733 
00734   for (l=0; l<nlens; l++) { /* Loop over lenses */
00735 
00736     /* [CCD/px] -> [MLA/mm] */
00737     snifs_optics_CCD2MLA(&optics, xdl[l],ydl[l], &xmla,&ymla, 
00738                          pixsize,lbdaref);
00739 
00740     /* ----- Data ------------------------------ */
00741 
00742     /* DEF: xd,yd = position along current spectra */
00743     for (k = 0; k<POLY_NPTS; k++) {
00744 
00745       /* [MLA/mm] -> [CCD/px] */
00746       snifs_optics_MLA2CCD(&optics, xmla,ymla, &xccd,&yccd, 
00747                            pixsize,lambda[k],ORDER);
00748 
00749       /* Relative offset */
00750       xd[k] = xccd - offset[0]; yd[k] = yccd - offset[1];
00751 
00752       if (LOCAL)                 /* Apply local adjustment to mask position */
00753         xd[k] += val_poly_nag(lambda[k], dxcoeff[l], dxnc, 
00754                               lbda_inf, lbda_sup, &status);
00755       /* By construction, lbda_inf <= lambda[k] <= lbda_sup */
00756     }
00757 
00758     /* ----- Limits ------------------------------ */
00759 
00760     xlo[l] = xd[POLY_NPTS-1]; xup[l] = xd[0]; /* Blue is up */
00761     ylo[l] = yd[POLY_NPTS-1]; yup[l] = yd[0]; 
00762 
00763     /* Add a pixel at both ends */
00764     ylo[l] -= inframe.stepy; yup[l] += inframe.stepy;
00765 
00766     if (LENS && l==l_debug) {
00767       print_msg("Spectrum #%d on CCD: %.2fx%.2f-%.2fx%.2f",
00768                 nol[l],xlo[l],ylo[l],xup[l],yup[l]);
00769     }
00770 
00771     if (WEIGHT) {
00772       /* [MLA/mm] -> [CCD/px]: Total spectral extension (no local correction) */
00773       snifs_optics_MLA2CCD(&optics, xmla,ymla, &xccd,&yccd, 
00774                            pixsize,local_lbda_inf_total[l],ORDER);
00775       xupt[l] = xccd; yupt[l] = yccd;
00776       snifs_optics_MLA2CCD(&optics, xmla,ymla, &xccd,&yccd, 
00777                            pixsize,local_lbda_sup_total[l],ORDER);
00778       xlot[l] = xccd; ylot[l] = yccd;
00779 
00780       /* Relative offsets */
00781       xupt[l] -= offset[0]; yupt[l] -= offset[1];
00782       xlot[l] -= offset[0]; ylot[l] -= offset[1];
00783 
00784       /* Add a pixel at both ends */
00785       ylot[l] -= inframe.stepy; yupt[l] += inframe.stepy;
00786     }
00787 
00788     /* ----- Fit xd(yd) ------------------------------ */
00789 
00790     status = Fit_polynom2(yd,xd,POLY_NPTS,XOFY_DEGMAX,
00791                           &xofy_deg[l],xofy_coeffs[l],Dx,&rms); /* DEBUG DEBUG DEBUG */
00792 
00793     if (LENS && l==l_debug) {
00794       print_msg("x = f(y):    status=%d, degree=%d, RMS=%g",
00795                 status,xofy_deg[l],rms);
00796       for (ldbg=k=0; k<POLY_NPTS; k++,ldbg++) {
00797         WR_tbl(&dbg_tab1,ldbg,coly1dbg,&yd[k]);
00798         WR_tbl(&dbg_tab1,ldbg,colxydbg,&xd[k]);
00799 
00800         val = Polynom(xofy_deg[l],xofy_coeffs[l],yd[k]); 
00801         WR_tbl(&dbg_tab1,ldbg,colxfdbg,&val);
00802         if (MEGA_VERBOSE)
00803           print_msg("%d/%d: x = %7.2f, y = %7.2f, xf = %7.2f, dx = %12.7f",
00804                     k+1,POLY_NPTS,xd[k],yd[k],val,xd[k]-val);
00805 
00806         val = xd[k] - val;
00807         WR_tbl(&dbg_tab1,ldbg,coldxdbg,&val);
00808       }
00809     }
00810 
00811     /* ----- Fit yd(lambda) ------------------------------ */
00812 
00813     status = Fit_polynom2(lambda,yd,POLY_NPTS,YOFL_DEGMAX,
00814                           &yofl_deg[l],yofl_coeffs[l],Dx,&rms); /* DEBUG DEBUG DEBUG */
00815 
00816     if (LENS && l==l_debug) {
00817       print_msg("y = f(lbda): status=%d, degree=%d, RMS=%g",
00818                 status,yofl_deg[l],rms);
00819       for (k=0,ldbg=0; k<POLY_NPTS; k++,ldbg++) {
00820         WR_tbl(&dbg_tab2,ldbg,collldbg,&lambda[k]);
00821         WR_tbl(&dbg_tab2,ldbg,colyldbg,&yd[k]);
00822 
00823         val = Polynom(yofl_deg[l],yofl_coeffs[l],lambda[k]); 
00824         WR_tbl(&dbg_tab2,ldbg,colyfdbg,&val);
00825         val = yd[k] - val;
00826         WR_tbl(&dbg_tab2,ldbg,coldydbg,&val);
00827       }
00828     }
00829 
00830     /* ----- Fit lambda(yd) ------------------------------ */
00831 
00832     status = Fit_polynom2(yd,lambda,POLY_NPTS,LOFY_DEGMAX,
00833                           &lofy_deg[l],lofy_coeffs[l],Dx,&rms); /* DEBUG DEBUG DEBUG */
00834 
00835     if (LENS && l==l_debug) {
00836       print_msg("lbda = f(y): status=%d, degree=%d, RMS=%g",
00837                 status,lofy_deg[l],rms);
00838       for (k=0,ldbg=0; k<POLY_NPTS; k++,ldbg++) {
00839         WR_tbl(&dbg_tab3,ldbg,coly2dbg,&yd[k]);
00840         WR_tbl(&dbg_tab3,ldbg,collydbg,&lambda[k]);
00841 
00842         val = Polynom(lofy_deg[l],lofy_coeffs[l],yd[k]); 
00843         WR_tbl(&dbg_tab3,ldbg,collfdbg,&val);
00844         val = lambda[k] - val;
00845         WR_tbl(&dbg_tab3,ldbg,coldldbg,&val);
00846       }
00847     }
00848 
00849     /* Spectral range of current spectrum */
00850     /* YC: Actually, I don't see the need of this section, except for 
00851        truncating the spectrum of *one* pixel on each side (which could 
00852        be done in an easier way...) */
00853     lbda_min[l] = lbda_inf;
00854     lbda_max[l] = lbda_sup;
00855     /* Blue is up */
00856     i = 0;             while ((yd[i] > yup[l]) && (i<POLY_NPTS)) i++;
00857     lbda_min[l] = lambda[i];
00858     i = POLY_NPTS - 1; while ((yd[i] < ylo[l]) && (i>=0))        i--;
00859     lbda_max[l] = lambda[i];
00860 
00861   } /* Loop over lenses (polynomial fits) */
00862 
00863   if (LENS) {
00864     write_file_class(&dbg_tab1,DEBUG_FILE); close_table(&dbg_tab1);
00865     write_file_class(&dbg_tab2,DEBUG_FILE); close_table(&dbg_tab2);
00866     write_file_class(&dbg_tab3,DEBUG_FILE); close_table(&dbg_tab3);
00867   }
00868 
00869   if (WEIGHT) {
00870     interspec = lens_size_pix*ABS(sin(optics.MLA.angle));
00871     toofar    = TOOFAR*interspec;
00872 
00873     print_msg("Inter-spectra distance: %.2f pixels",interspec);
00874   }
00875 
00876   if (LENS) {
00877     sprintf(dbgname, "%s_xdisp", DBG_PREFIX);
00878     print_msg("DEBUG: Creating datacube '%s' with X-disp. profiles...",
00879               dbgname);
00880     disable_erase_flag();
00881     create_tiger_frame(&dbgcube,dbgname,-1,-1,inframe.stepx,FLOAT,dbgname,
00882                        "Extraction check","");
00883 
00884     create_table(&dbg_tab1,dbgname,-1,4,'Q',ident);
00885     colnodbg = create_col(&dbg_tab1,LAB_COL_NO,INT,'N',"I4", "Lens Number");
00886     colxdbg  = create_col(&dbg_tab1,"X",    FLOAT,'N',"F6.2","px");
00887     colydbg  = create_col(&dbg_tab1,"Y",    FLOAT,'N',"F6.2","px");
00888     colldbg  = create_col(&dbg_tab1,"LBDA", FLOAT,'N',"F6.2","px");
00889 
00890     restore_erase_flag();
00891   }
00892 
00893   /* ##### EXTRACTION ############################## */
00894 
00895   if (LENS) {
00896     l0 = l_debug; 
00897     ln = 1; 
00898     dbgno = 1;
00899   } 
00900   else {
00901     l0 = 0;      /* first lens index */
00902     ln = nlens;  /* number of lenses */
00903   }
00904 
00905   signal   = (float *)malloc(max_npix*sizeof(float)); /* Extracted signal */
00906   variance = (float *)malloc(max_npix*sizeof(float)); /* Associated variance */
00907 
00908   reset_print_progress();
00909   for (ltab=lc=0; lc<ln; lc++) { /* Loop over lenses */
00910 
00911     l = l0 + lc;
00912     print_progress("Extraction in progress", 100.*(lc+1.)/ln, 1.);
00913 
00914     /* ===== Selection ============================== */
00915 
00916     /* If the spectrum is fully outside the CCD, don't extract it! */
00917     if ( (!in_frame(&inframe,xlo[l],ylo[l])) && 
00918          (!in_frame(&inframe,xup[l],yup[l])) ) {
00919       if (DEBUG) 
00920         print_msg("Spectrum #%d (%.1fx%.1f-%.1fx%.1f) is "
00921                   "outside the CCD, skipping",
00922                   nol[l],xlo[l],ylo[l],xup[l],yup[l]);
00923       continue;
00924     }
00925 
00926     if (nol[l] == 0) {
00927       print_error("ARGH, there's a lens #0!!! (l=%d, %.1fx%.1f)",
00928                   l,xdl[l],ydl[l]);
00929       continue;
00930     }
00931 
00932     if (WEIGHT && (sigcoeff[l][0] < SMIN || sigcoeff[l][0] > SMAX)) {
00933       if (DEBUG) print_warning("Lens #%d discarded (sigma = %f)",
00934                                nol[l],sigcoeff[l][0]);
00935       continue;
00936     }
00937 
00938     /* ===== Truncation ============================== */
00939 
00940     /* Discard the spectra truncated by the left/right edges of the CCD
00941        but keep the spectra truncated by the upper/lower edges of the CCD
00942        and update the end positions */
00943     pixel_frame(&inframe,xup[l],yup[l],&iima,&jima);
00944     if ((iima < 0) || (iima > inframe.nx-1)) {
00945       print_msg("Spectrum #%d (Up:%.1fx%.1f) is "
00946                 "partially outside the CCD (left/right), skipping",
00947                 nol[l],xup[l],yup[l]);
00948       continue;
00949     }
00950     jima = MIN(MAX(0,jima),inframe.ny-1);
00951     coord_frame(&inframe,iima,jima,&xup[l],&yup[l]);
00952 
00953     pixel_frame(&inframe,xlo[l],ylo[l],&iima,&jima);
00954     if ((iima < 0) || (iima > inframe.nx-1)) {
00955       print_msg("Spectrum #%d (Low:%.1fx%.1f) is "
00956                 "partially outside the CCD (left/right), skipping",
00957                 nol[l],xlo[l],ylo[l]);
00958       continue;
00959     }
00960     jima = MIN(MAX(0,jima),inframe.ny-1);
00961     coord_frame(&inframe,iima,jima,&xlo[l],&ylo[l]);
00962 
00963     /* Update nb of pixel for current spectrum */
00964     npix = rint(ABS(yup[l]-ylo[l])/inframe.stepy) + 1;
00965 
00966     /* Compute CCD-position of the center (npix/2) of the spectrum */ 
00967     yref = ylo[l] + npix*inframe.stepy/2.;
00968     xref = Polynom(xofy_deg[l],xofy_coeffs[l],yref);
00969 
00970     if (LENS && l==l_debug) {
00971       print_msg("Spectrum #%d (%d px): %.2fx%.2f-%.2fx%.2f",
00972                 nol[l],npix,xlo[l],ylo[l],xup[l],yup[l]);
00973 
00974       /* set coordinates relative to central position of extracted
00975          spectrum, and not anymore in absolute CCD coordinates. */
00976       print_msg("Use position x=%.2f (y=%.2f) as reference in x "
00977                 "for %s", xref,yref,dbgname);
00978     }
00979 
00980     /* ===== Summation ============================== */
00981 
00982     /* Compute flux by summation for each yd */
00983     flag_out = 0;
00984     lg_min   =-1;
00985     ld_min   =-1;
00986 
00987     /* TO BE UPDATED:
00988        jima is coming from pixel_frame(&inframe,xlo[l],ylo[l],&iima,&jima); */
00989 
00990     for (jspec=0; jspec<npix; jspec++, jima++) {
00991 
00992       y = ylo[l] + jspec*inframe.stepy;
00993       x = Polynom(xofy_deg[l],xofy_coeffs[l],y);    /* get x = f(yd) */
00994 
00995       /* ----- Extraction window ------------------------------ */
00996 
00997       xlef = x - half_width;
00998       xrig = x + half_width;
00999 
01000       if (xlef < inframe.startx || xrig > inframe.endx) { 
01001         /* spectrum box is outside frame (at .5 px) */
01002         flag_out = 1;
01003         break;
01004       }
01005 
01006       /* Pixel X-coordinates of the extraction window (0-based) */
01007       flef = (xlef - inframe.startx)/inframe.stepx;
01008       frig = (xrig - inframe.startx)/inframe.stepx;
01009 
01010       /* Full pixels in <ilef,irig> */
01011       ilef = (int)(flef+0.5) + 1;
01012       irig = (int)(frig+0.5) - 1;
01013       /* Fraction of pixels for ilef-1 and irig+1 */
01014       lfrac = (ilef-0.5) - flef;
01015       rfrac = frig - (irig+0.5);
01016 
01017       signal[jspec]   = 0; 
01018       variance[jspec] = 0;
01019 
01020       if (!WEIGHT) { 
01021 
01022         /* ----- Raw extraction ------------------------------ */
01023 
01024         /* 1: sum over complete pixels */
01025 
01026         for (iima=ilef; iima<=irig; iima++) {
01027           signal[jspec] += RD_frame(&inframe,iima,jima); /* Input frame is in e- */
01028           /* noise estimate */
01029           if (NOISE) {
01030             if (VARIANCE) 
01031               variance[jspec] += RD_frame(&invarframe,iima,jima);
01032             else
01033               variance[jspec] += RD_frame(&inframe,iima,jima) + RoN2;
01034           }
01035           
01036         }
01037 
01038         /* 2: add boundaries */
01039 
01040         signal[jspec] += 
01041           lfrac*RD_frame(&inframe,ilef-1,jima) + 
01042           rfrac*RD_frame(&inframe,irig+1,jima);
01043         if (NOISE) {
01044           if (VARIANCE)
01045             variance[jspec] += 
01046               SQ(lfrac)*RD_frame(&invarframe,ilef-1,jima) + 
01047               SQ(rfrac)*RD_frame(&invarframe,irig+1,jima);
01048           else
01049             variance[jspec] += 
01050               SQ(lfrac)*( RD_frame(&inframe,ilef-1,jima) + RoN2 ) + 
01051               SQ(rfrac)*( RD_frame(&inframe,irig+1,jima) + RoN2 );
01052         }
01053 
01054         if (LENS) { /* LENS debug mode */
01055 
01056           lbda = Polynom(lofy_deg[l],lofy_coeffs[l],y);
01057           if (MEGA_VERBOSE) {
01058             print_msg("line %4d (x=%7.2f y=%7.2f) pixels i=%4d->%4d j=%4d: "
01059                       "lbda=%7.2f, I=%f",l+1,x,y,ilef,irig,jima,lbda,signal[jspec]);
01060           }
01061 
01062           WR_tbl(&dbg_tab1, dbgno-1, colnodbg, &dbgno);
01063           WR_tbl(&dbg_tab1, dbgno-1, colxdbg,  &x);  
01064           WR_tbl(&dbg_tab1, dbgno-1, colydbg,  &y);
01065           WR_tbl(&dbg_tab1, dbgno-1, colldbg,  &lbda);
01066 
01067           /* Debug window limits */
01068           ibegdbg = MAX((int)(ilef-DBGSIDEWIN*half_width-0.5),0);
01069           ienddbg = MIN((int)(irig+DBGSIDEWIN*half_width+0.5),inframe.nx-1);
01070 
01071           xi = inframe.startx + ibegdbg*inframe.stepx;
01072           /* Set coord. relative to central position of extracted spectrum */
01073           xi -= xref; 
01074           init_new_tiger_spec(&dbgcube,&dbg_spec, ienddbg-ibegdbg+1,xi);
01075           init_new_tiger_spec(&dbgcube,&dbg_noise,ienddbg-ibegdbg+1,xi);
01076 
01077           maxval = -MAXFLOAT;
01078           minval =  MAXFLOAT;
01079           for (iima=ibegdbg; iima<=ienddbg; iima++) {
01080             /* Cross-dispersion profile in signal spectrum */
01081             val = RD_frame(&inframe,iima,jima);
01082             WR_spec(&dbg_spec,iima-ibegdbg,val);
01083             maxval = MAX(maxval,val);
01084             minval = MIN(minval,val);
01085 
01086             /* Cross-dispersion extraction window in noise spectrum */
01087             if (iima>=ilef && iima<=irig) val = 1;
01088             else if (iima==ilef-1)        val = lfrac;
01089             else if (iima==irig+1)        val = rfrac;
01090             else                          val = 0;
01091             WR_spec(&dbg_noise,iima-ibegdbg,val);
01092           }
01093 
01094           /* Normalization */
01095           minval = MIN(maxval/1e2,minval); /* For log-display */
01096           if (maxval > 0) {
01097             for (iima=ibegdbg; iima<=ienddbg; iima++) {
01098               val = RD_spec(&dbg_noise,iima-ibegdbg);
01099               WR_spec(&dbg_noise,iima-ibegdbg,
01100                       (val!=0)? val*maxval:minval);
01101             }
01102           }
01103 
01104           put_tiger_spec(&dbgcube,&dbg_spec,&dbg_noise,dbgno);
01105           dbgno++;
01106 
01107         } /* LENS debug mode */
01108 
01109       } /* Raw extraction */
01110 
01111       else { 
01112 
01113         /* ----- Weighted extraction ------------------------------ */
01114 
01115         /* ..... Neighbors .............................. */
01116 
01117         /* Test that previous neighbours (if any) are undoubtfully good */
01118         if (lg_min >= 0)
01119           if ((y < ylo[lg_min]) || (y > yup[lg_min])) lg_min = -1;
01120         if (ld_min >= 0) 
01121           if ((y < ylo[ld_min]) || (y > yup[ld_min])) ld_min = -1;
01122 
01123         /* If one of the neighbour is missing, find the 2 neighbours again */
01124         if (lg_min < 0 || ld_min < 0) {
01125           nlneighbors = 0;
01126           nrneighbors = 0;
01127           for (i=0; i<nlens; i++) { /* Loop over lenses (neighbor) */
01128             /* Simple test on x to remove most of the lenses */
01129             if (x < (xlo[i] - toofar) || 
01130                 x > (xlo[i] + toofar)) continue; /* skip this lens */
01131             /* ylot/yupt are now the *real* local spectrum edges */
01132             if (y < (ylot[i] - XTRALENGTH) || 
01133                 y > (yupt[i] + XTRALENGTH)) continue; /* skip this lens */
01134             /* If (nol[i] == 0) continue; *//* take care of lenses #0 */
01135             /* Actually, keep lenses #0 as possible neighbors, as they
01136                can model spectra undetected in the mask
01137                creation/adjustment process but still there as source
01138                of pollution. Anyway, if there is indeed no real
01139                neighbor and you find a lens #0 as neighbor, you will
01140                remove nothing... */
01141 
01142             dx = x - Polynom(xofy_deg[i],xofy_coeffs[i],y);
01143             /* Look for neighbor spectra */
01144             /* CAUTION: We may find on each side 2 spectra, one up and one
01145                down, since the 2 total ranges may overlap in the "dark
01146                zone" => Find the good spectrum later on */
01147 
01148             /* <- on the left */
01149             if ((dx > NEIGHBRANGE_C*interspec) && 
01150                 (dx < NEIGHBRANGE_F*interspec)) {
01151               if (nlneighbors == 2) { /* never happen w/ good NEIGHBRANGE */
01152                 print_warning("Already found 2 left neighbors (%d and %d) "
01153                               "for lens #%d, skipping new one %d",
01154                               nol[lneighbors[0]],nol[lneighbors[1]],
01155                               nol[l],nol[i]);
01156               } 
01157               else { 
01158                 lneighbors[nlneighbors] = i;
01159                 nlneighbors++;
01160               }
01161             }
01162 
01163             /* -> on the right */
01164             if ((dx > -NEIGHBRANGE_F*interspec) && 
01165                 (dx < -NEIGHBRANGE_C*interspec)) {
01166               if (nrneighbors == 2) { /* never happen w/ good NEIGHBRANGE */
01167                 print_warning("Already found 2 right neighbors (%d and %d) "
01168                               "for lens #%d, skipping new one %d",
01169                               nol[rneighbors[0]],nol[rneighbors[1]],
01170                               nol[l],nol[i]);
01171               } 
01172               else { 
01173                 rneighbors[nrneighbors] = i;
01174                 nrneighbors++;
01175               }
01176             }
01177           } /* Loop over lenses (neighbor) */
01178 
01179           /* Find the good neighbors */
01180           switch (nlneighbors) {
01181           case 0: lg_min = -1;                                       break; /* None */
01182           case 1: lg_min = lneighbors[0];                            break; /* One */
01183           case 2: lg_min = good_neighbor(ylot, yupt, lneighbors, y); break; /* More */
01184           }
01185           switch (nrneighbors) {
01186           case 0: ld_min = -1;                                       break; /* None */
01187           case 1: ld_min = rneighbors[0];                            break; /* One */
01188           case 2: ld_min = good_neighbor(ylot, yupt, rneighbors, y); break; /* More */
01189           }
01190 
01191         } /* Looking for neighbors */
01192 
01193         /* ..... 1: Sum over complete pixels .............................. */
01194 
01195         /* Chromatic model of cross-dispersion sigma */
01196 
01197         for (i=0; i<3; i++) { /* Loop over neighbors */
01198           switch (i) {
01199           case 0: n = lg_min; break;         /* Left neighbor */
01200           case 1: n = l;      break;         /* Central lens */
01201           case 2: n = ld_min; break;         /* Right neighbor */
01202           } 
01203 
01204           if (n>=0) {                        /* Valid lens */
01205             /* get lbda = f(yd) */
01206             lbdasig[i] = Polynom(lofy_deg[n],lofy_coeffs[n],y);
01207 
01208             /* Estimate dispersion from chromatic linear model */
01209             if (sigcoeff[n][0] > 0 && sigcoeff[n][0] < SMAX &&
01210                 lbda_inf <= lbdasig[i] && lbdasig[i] <= lbda_sup) /* Too restrictive! */
01211               sigma[i] = MAX(SMIN,val_poly_nag(lbdasig[i], sigcoeff[n], signc, 
01212                                                lbda_inf, lbda_sup, &status));
01213             else /* in case of pb, take average dispersion from central lens */
01214               sigma[i] = sigcoeff[l][0];
01215 
01216             for (k=0; k<glnggeo; k++) {
01217               for (p=0; p<glngpsf; p++) {
01218                 glScs[n][k][p] = hypot(glSpsf[k][p],sigma[i]);
01219                 glIcs[n][k][p] = M_SQRT2PI*glIpsf[k][p]*glSpsf[k][p]*sigma[i]/
01220                   glScs[n][k][p];
01221               }
01222             }
01223           }
01224 
01225         } /* Loop over neighbors */
01226 
01227         if (LENS && MEGA_VERBOSE) {          /* verbose mode */
01228           if (lg_min >= 0) xinf = Polynom(xofy_deg[lg_min],xofy_coeffs[lg_min],y); 
01229           if (ld_min >= 0) xsup = Polynom(xofy_deg[ld_min],xofy_coeffs[ld_min],y); 
01230 
01231           if (lg_min >= 0 && ld_min >= 0) {
01232             print_msg("%3d (y = %.0f): #%4d @%7.2f (x = %7.2f) "
01233                       "<- #%4d @%7.2f (x = %7.2f) -> "
01234                       "#%4d @%7.2f (x = %7.2f)", 
01235                       dbgno,y,nol[lg_min],lbdasig[0],xinf,
01236                       nol[l],lbdasig[1],x,
01237                       nol[ld_min],lbdasig[2],xsup); 
01238           } else if (lg_min >= 0) {
01239             print_msg("%3d (y = %.0f): #%4d @%7.2f (x = %7.2f) "
01240                       "<- #%4d @%7.2f (x = %7.2f) -> XXXXX", 
01241                       dbgno,y,nol[lg_min],lbdasig[0],xinf,
01242                       nol[l],lbdasig[1],x); 
01243           } else if (ld_min >= 0) {
01244             print_msg("%3d (y = %.0f): XXXXX <- #%4d @%7.2f (x = %7.2f) -> "
01245                       "#%4d @%7.2f (x = %7.2f)", 
01246                       dbgno,y,nol[l],lbdasig[1],x,
01247                       nol[ld_min],lbdasig[2],xsup); 
01248           } 
01249           else {
01250             print_msg("%3d (y = %.0f): XXXXX <- #%4d @%7.2f (x = %7.2f) "
01251                       "-> XXXXX", dbgno,y,nol[l],lbdasig[1],x); 
01252           }
01253 
01254           /* 
01255           if (lg_min >= 0 && ld_min >= 0) { 
01256             print_msg("Sigma:  inf = %9.4f, cen = %9.4f, sup = %9.4f", 
01257                     sigma[0], sigma[1], sigma[2]);
01258           } else if (lg_min >= 0) {  
01259             print_msg("Sigma:  inf = %9.4f, cen = %9.4f, sup = ---", 
01260                     sigma[0], sigma[1]);
01261           } else if (ld_min >= 0) { 
01262             print_msg("Sigma:  inf = ---, cen = %9.4f, sup = %9.4f", 
01263                     sigma[1], sigma[2]);
01264           } 
01265           else { 
01266             print_msg("Sigma:  inf = ---, cen = %9.4f, sup = ---", sigma[1]);
01267           }
01268           */
01269         }
01270 
01271         /* if no LINFIT: use median (or mean) estimate for intensity,
01272            ie, normalize the model profile Mpup_PSF by the median (or mean) of the 
01273            3 central pixel ratio q = frame/model */
01274         /* if LINFIT: if no 1st guess for the linear fit (q=0), 
01275            take the max estimate as first guess */
01276 
01277         if (!LINFIT || 
01278             (LINFIT && (qinf == 0 || qcen == 0 || qsup == 0))) {
01279 
01280           /* === Median/Max estimate === */
01281 
01282           qinf = 0; qcen = 0; qsup = 0;
01283           xinf = 0; xsup = 0;
01284 
01285           /* Central contribution: rough estimate (including neighbor pollution) */
01286           icen = (int)((x - inframe.startx)/inframe.stepx + 0.5);
01287           /* Normalization */
01288           for (n=-1; n<=1; n++) {            /* loop over 3 central pixels */
01289             xi = inframe.startx + (icen+n)*inframe.stepx;
01290             qloc[n+1] = RD_frame(&inframe,icen+n,jima) / 
01291               Mpup_PSF(l,(xi - x)); /* ratio real/model */
01292           }
01293           // qcen = (qloc[0]+qloc[1]+qloc[2])/3.; /* normalize at mean */
01294           indexx(3,qloc,iq); qcen = qloc[iq[1]]; /* normalize at median */
01295           qcen = MAX(0,qcen);                /* DEBUG DEBUG DEBUG */
01296 
01297           /* Left contribution: central contribution removed */
01298           ibegdbg = MAX((int)(ilef-DBGSIDEWIN*half_width-0.5),0);
01299           if (lg_min >= 0) {
01300             xinf = Polynom(xofy_deg[lg_min],xofy_coeffs[lg_min],y); 
01301             iinf = (int)((xinf - inframe.startx) / inframe.stepx + 0.5);
01302 
01303             if (LINFIT) { /* LINFIT: start of the fitting window */
01304               ibeg = (int)((xinf - inframe.startx) / inframe.stepx - interspec/2 + 0.5);
01305               ibeg = MAX(ibeg,0); 
01306             }
01307             /* debug window left limit */
01308             ibegdbg = MAX((int)(iinf-DBGSIDEWIN*half_width-0.5),0);
01309 
01310             if (OVERLAP && iinf-1>=0) {
01311               /* left neighbour median normalization */
01312               for (n=-1; n<=1; n++) {        /* loop over 3 central pixels */
01313                 xi = inframe.startx + (iinf+n)*inframe.stepx;
01314                 qloc[n+1] = ( RD_frame(&inframe,iinf+n,jima) - qcen*Mpup_PSF(l,(xi - x)) ) / 
01315                   Mpup_PSF(lg_min,(xi - xinf)); /* ratio real/model */
01316               }
01317               // qinf = (qloc[0]+qloc[1]+qloc[2])/3.; /* normalize at mean */
01318               // qinf = RD_frame(&inframe,iinf,jima)/Mpup_PSF(lg_min,(inframe.startx + iinf*inframe.stepx - xinf)); /* normalize at max */
01319               indexx(3,qloc,iq); qinf = qloc[iq[1]]; /* normalize at median */
01320               qinf = MAX(0,qinf);                /* DEBUG DEBUG DEBUG */
01321             }
01322           }
01323 
01324           /* Right contribution: central contribution removed  */
01325           ienddbg = MIN((int)(irig+DBGSIDEWIN*half_width+0.5),inframe.nx-1);
01326           if (ld_min >= 0) {
01327             xsup = Polynom(xofy_deg[ld_min],xofy_coeffs[ld_min],y); 
01328             isup = (int)((xsup - inframe.startx)/inframe.stepx + 0.5);
01329 
01330             if (LINFIT) { /* LINFIT: end of the fitting window */
01331               iend = (int)((xsup - inframe.startx) / inframe.stepx + interspec/2 + 0.5);
01332               iend = MIN(iend,inframe.nx-1);
01333             }
01334             /* debug window righ limit */
01335             ienddbg = MIN((int)(isup+DBGSIDEWIN*half_width+0.5),inframe.nx-1);
01336 
01337             if (OVERLAP && isup+1 < inframe.nx) {
01338               /* right neighbour median normalization */
01339               for (n=-1; n<=1; n++) {        /* loop over 3 central pixels */
01340                 xi = inframe.startx + (isup+n)*inframe.stepx;
01341                 qloc[n+1] = ( RD_frame(&inframe,isup+n,jima) - qcen*Mpup_PSF(l,(xi - x)) ) / 
01342                   Mpup_PSF(ld_min,(xi - xsup)); /* ratio real/model */
01343               }
01344               // qsup = (qloc[0]+qloc[1]+qloc[2])/3.; /* normalize at mean */ 
01345               indexx(3,qloc,iq); qsup = qloc[iq[1]]; /* normalize at median */
01346               qsup = MAX(0,qsup);                /* DEBUG DEBUG DEBUG */
01347             }
01348           }
01349 
01350           /* Update central contribution: neighbour contributions removed */
01351           for (n=-1; n<=1; n++) {
01352             xi = inframe.startx + (icen+n)*inframe.stepx;
01353             qloc[n+1] = ( RD_frame(&inframe,icen+n,jima) /* ratio (real-neighbors)/model */
01354                           - qinf*Mpup_PSF(lg_min,(xi-xinf))
01355                           - qsup*Mpup_PSF(ld_min,(xi-xsup)) ) / Mpup_PSF(l,(xi - x));
01356           }
01357           // qcen = (qloc[0]+qloc[1]+qloc[2])/3.; /* normalize at mean */
01358           indexx(3,qloc,iq); qcen = qloc[iq[1]]; /* normalize at median */
01359           qcen = MAX(0,qcen);                /* DEBUG DEBUG DEBUG */
01360 
01361           /* if (LENS && MEGA_VERBOSE) {
01362              print_msg("Median estimate: qinf = %f, qcen = %f, qsup = %f",
01363              qinf,qcen,qsup);
01364              } */
01365 
01366         } /* end of Median estimate */
01367 
01368         if (LINFIT && qinf && qcen && qsup) {
01369 
01370           /* === Linear fitting estimate === */
01371 
01372           /* LINFIT: Fit the intensity to take into account the pollutions */
01373           /* HYP: the current spectrum has 2 valid neighbours */
01374   
01375           /* Initialization */
01376   
01377           nvar    = 3;                       /* fit 3 peak intensities */
01378           ibound  = 2;                       /* all intensities are > 0 */
01379           npixfit = iend - ibeg + 1;         /* nb of pixels in fit */
01380 
01381           liw = nvar + 2;
01382           lw  = MAX(nvar*(nvar-1)/2+12*nvar,13);
01383 
01384           iw     =   (long *)malloc(liw     *sizeof(long));
01385           glnpeaks = (long *)malloc((nvar+1)*sizeof(long));
01386           q      = (double *)malloc(nvar            *sizeof(double));
01387           bl     = (double *)malloc(nvar            *sizeof(double));
01388           bu     = (double *)malloc(nvar            *sizeof(double));
01389           w      = (double *)malloc(lw              *sizeof(double));
01390           glpeaks= (double *)malloc((2*npixfit+nvar)*sizeof(double));
01391           glnpeaks[0] = npixfit;
01392 
01393           for (n=0; n<npixfit; n++) {
01394             glpeaks[n] = inframe.startx + (ibeg+n)*inframe.stepx;
01395             glpeaks[npixfit + n] = RD_frame(&inframe,(ibeg+n),jima);
01396           }
01397 
01398           /* Lens numbers and positions */
01399           glnpeaks[1] = lg_min;
01400           glnpeaks[2] = l;
01401           glnpeaks[3] = ld_min;
01402 
01403           glpeaks[2*npixfit  ] = xinf;
01404           glpeaks[2*npixfit+1] = x;
01405           glpeaks[2*npixfit+2] = xsup;
01406 
01407           /* Initial estimates (from previous fit) */
01408           q[0] = qinf; q[1] = qcen; q[2] = qsup;
01409 
01410           ifail = -1;
01411           e04jaf(&nvar,&ibound,bl,bu,q,&res,iw,&liw,w,&lw,&ifail);
01412 
01413           if (ifail) {
01414             print_warning("Fit: ifail=%d in intensity fit at pixel %d, "
01415                           "keeping initial guess from previous fit", 
01416                           (int)ifail,dbgno);
01417             /* keep initial guess as solution */
01418           } 
01419           else {
01420             /* Final solution */
01421             qinf = q[0];
01422             qcen = q[1];
01423             qsup = q[2];
01424           }
01425 
01426           free(iw); free(glnpeaks); free(q);
01427           free(bl); free(bu); free(w); free(glpeaks);
01428 
01429           if (LENS && MEGA_VERBOSE) {
01430             print_msg("Linear fit: qinf = %f, qcen = %f, qsup = %f",
01431                       qinf,qcen,qsup);
01432           }
01433 
01434         } /* LINFIT */
01435 
01436         if (LENS) {                          /* LENS debug mode */
01437           WR_tbl(&dbg_tab1, dbgno-1, colnodbg, &dbgno);
01438           WR_tbl(&dbg_tab1, dbgno-1, colxdbg,  &x);  
01439           WR_tbl(&dbg_tab1, dbgno-1, colydbg,  &y);
01440           fval[0] = lbdasig[1];
01441           WR_tbl(&dbg_tab1, dbgno-1, colldbg,  fval);
01442 
01443           xi = inframe.startx + ibegdbg*inframe.stepx;
01444           /* set coord. relative to central position of extracted spectrum */
01445           xi -= xref; 
01446           init_new_tiger_spec(&dbgcube,&dbg_spec,(ienddbg-ibegdbg+1),xi);
01447           init_new_tiger_spec(&dbgcube,&dbg_noise,(ienddbg-ibegdbg+1),xi);
01448 
01449           for (iima=ibegdbg; iima<=ienddbg; iima++) {
01450             /* Cross-dispersion profile in signal spectrum */
01451             WR_spec(&dbg_spec,iima-ibegdbg,RD_frame(&inframe,iima,jima));
01452 
01453             /* Cross-dispersion model in noise spectrum */
01454             xi = inframe.startx + iima*inframe.stepx;
01455             val = qinf * Mpup_PSF(lg_min, xi-xinf)
01456               +   qcen * Mpup_PSF(l,      xi-x   ) 
01457               +   qsup * Mpup_PSF(ld_min, xi-xsup);
01458             WR_spec(&dbg_noise,iima-ibegdbg,val);
01459           }
01460 
01461           put_tiger_spec(&dbgcube,&dbg_spec,&dbg_noise,dbgno);
01462           dbgno++;
01463         } /* LENS debug mode */
01464 
01465         s_pcen = 0;   /* Integral of the unnormalized central profile model */
01466         for (iima=ilef; iima<=irig; iima++) {
01467           xi = inframe.startx + iima*inframe.stepx;
01468           s_pcen += Mpup_PSF(l,(xi - x));
01469         }
01470 
01471         /* Eq. numbers refer to Robertson 1986 */
01472         s_vi = 0;
01473         for (iima=ilef; iima<=irig; iima++) {
01474           xi = inframe.startx + iima*inframe.stepx;
01475           pinf = Mpup_PSF(lg_min, xi-xinf);
01476           pcen = Mpup_PSF(l,      xi-x   );
01477           psup = Mpup_PSF(ld_min, xi-xsup);
01478           val = qinf*pinf + qcen*pcen + qsup*psup; /* X-disp. model (eq. II.8) */
01479           vi  = RoN2 + val;                  /* Denominator of eq. II.7 */
01480 
01481           pcen /= s_pcen;    /* Normalized central profile model (eq. II.9) */
01482           /* Unpolluted signal (Peff, eq. II.6) */
01483           val = RD_frame(&inframe,iima,jima) - qinf*pinf - qsup*psup;
01484           signal[jspec] += val * (pcen / vi); /* ??? */
01485 
01486           s_vi += SQ(pcen)/vi;               /* ??? */
01487         }
01488 
01489         /* ..... 2: Add boundaries .............................. */
01490 
01491         /* Left */
01492         xi   = inframe.startx + (ilef-1)*inframe.stepx;
01493         pinf = Mpup_PSF(lg_min,(xi-xinf));
01494         pcen = Mpup_PSF(l,     (xi-x));
01495         psup = Mpup_PSF(ld_min,(xi-xsup));
01496         vi   = RoN2 + (qinf*pinf + qcen*pcen + qsup*psup);
01497 
01498         pcen /= s_pcen;
01499         signal[jspec] += lfrac * 
01500           ( RD_frame(&inframe,ilef-1,jima) - qinf*pinf - qsup*psup )*pcen/vi;
01501 
01502         s_vi += SQ(pcen)/vi;
01503 
01504         /* Right */
01505         xi   = inframe.startx + (irig+1)*inframe.stepx;
01506         pinf = Mpup_PSF(lg_min,(xi-xinf));
01507         pcen = Mpup_PSF(l,     (xi-x));
01508         psup = Mpup_PSF(ld_min,(xi-xsup));
01509         vi   = RoN2 + (qinf*pinf + qcen*pcen + qsup*psup);
01510 
01511         pcen /= s_pcen;
01512         signal[jspec] += rfrac * 
01513           ( RD_frame(&inframe,irig+1,jima) - qinf*pinf - qsup*psup )*pcen/vi;
01514 
01515         s_vi += SQ(pcen)/vi;
01516 
01517         if (s_vi > 0) { /* Patch to prevent a negative variance ?!?... */
01518           signal[jspec] /= s_vi;
01519           if (NOISE) variance[jspec] = 1/s_vi;
01520         } 
01521         else { 
01522           print_error("Argh, negative variance at px %d of lens #%d",jspec,nol[l]);
01523           signal[jspec] = 0;
01524           if (NOISE) variance[jspec] = 0;
01525         }
01526 
01527       } /* Weighted summation */
01528 
01529     } /* Loop over spectrum pixels */
01530 
01531     /* ===== Interpolation & Storage ============================== */
01532 
01533     if (flag_out) continue; /* Skip if there is a problem */
01534 
01535     npix = (int)(lbda_max[l] - lbda_min[l])/outstep + 1;
01536 
01537     init_new_tiger_spec(&outcube,&spectrum,npix,lbda_min[l]);
01538     if (NOISE) init_new_tiger_spec(&outcube,&noise,npix,lbda_min[l]);
01539 
01540     /* Interpolate flux at the given wavelength, using flux(yd) */
01541     for (i=0; i<spectrum.npts; i++) {
01542       lbda = lbda_min[l] + i*outstep;
01543       y = Polynom(yofl_deg[l],yofl_coeffs[l],lbda); /* get y=f(lbda) */
01544 
01545       fk = (y - ylo[l])/inframe.stepy;
01546       k = (int)fk;
01547       lfrac = fk - k;
01548       rfrac = 1 - lfrac;
01549       sum = rfrac*signal[k] + lfrac*signal[k+1]; /* Linear interpolation */
01550 
01551       WR_spec(&spectrum,i,sum);
01552 
01553       if (NOISE) {
01554         lfrac *= lfrac;                      /* Quadractic interpolation */
01555         rfrac *= rfrac;
01556         varsum = (lfrac*variance[k+1] + rfrac*variance[k])/(lfrac + rfrac); 
01557         WR_spec(&noise,i,varsum);
01558       }
01559     }
01560 
01561     if (NOISE) put_tiger_spec(&outcube,&spectrum,&noise,nol[l]);
01562     else       put_tiger_spec(&outcube,&spectrum,NULL,  nol[l]);
01563 
01564     WR_tbl(&tab, ltab, colno, &nol[l]);
01565     WR_tbl(&tab, ltab, coli,  &I[l]);        /* position indices */
01566     WR_tbl(&tab, ltab, colj,  &J[l]);
01567     WR_tbl(&tab, ltab, colxnd, &xref);       /* real mean CCD-position */
01568     WR_tbl(&tab, ltab, colynd, &yref);
01569 
01570     /* Spatial coordinates: from undistorted positions */
01571     /* There is no "undistorted positions", so just use I,J and SNIFS_SAMP */
01572     alpha = I[l]*SNIFS_SAMPL;
01573     WR_tbl(&tab, ltab, cola, &alpha);
01574     WR_tbl(&tab, ltab, colao,&alpha);
01575 
01576     delta = J[l]*SNIFS_SAMPL;
01577     WR_tbl(&tab, ltab, cold, &delta);
01578     WR_tbl(&tab, ltab, coldo,&delta);
01579 
01580     ltab++;
01581 
01582   } /* Loop over lenses (extraction) */
01583 
01584   /* ##### CONCLUSIONS ############################## */
01585 
01586   if (DEBUG) print_msg("Updating and closing datacube %s and "
01587                        "associated table %s", outcube.name,tab.name);
01588 
01589   if (LENS) {
01590     write_file_class(&dbgcube,DEBUG_FILE);
01591     close_tiger_frame(&dbgcube);
01592     write_file_class(&dbg_tab1,DEBUG_FILE);
01593     close_table(&dbg_tab1);
01594   }
01595 
01596   write_file_class(&tab, TBL_FILE); /* write file class */
01597   close_table(&tab);
01598 
01599   /* Lens size */
01600   fval[0] = lens_size_pix;                   /* Used by remcosmic */
01601   fval[1] = lens_size_arcsec;
01602   fval[2] = 0.;                              /* Don't know what it was... */
01603   WR_desc(&outcube, LENSIZE3, FLOAT, 3, fval);
01604 //  WR_desc(&outcube, LENSIZE3, FLOAT, 3,  /* C99 compound literals (gcc 3.1) */
01605 //          (double[]){ lens_size_pix,         /* Used by remcosmic */
01606 //                      lens_size_arcsec, 
01607 //                      0. });                 /* Don't know what it was... */
01608 
01609   /* Wavelength range (useful & common) */
01610   dval[0] = optics.filter.inf_util;
01611   dval[1] = optics.filter.sup_util;
01612   dval[2] = optics.filter.inf_util; /* Will contain common wavelength domain */
01613   dval[3] = optics.filter.sup_util;
01614   WR_desc(&outcube, LRANGE4, DOUBLE, 4, dval);
01615 //  WR_desc(&outcube, LRANGE4, DOUBLE, 4,  /* C99 compound literals (gcc 3.1) */
01616 //          (double[]) { optics.filter.inf_util,
01617 //                       optics.filter.sup_util,
01618 //                       optics.filter.inf_util,
01619 //                       optics.filter.sup_util });
01620 
01621   /* Extraction-related quantities */
01622   WR_desc(&outcube,CALARC,CHAR,80,arcname);
01623   if (CAL2) WR_desc(&outcube,CALARC2,CHAR,80,arc2name);
01624   WR_desc(&outcube, EXTMASK, CHAR, 80,maskname);
01625   WR_desc(&outcube, OFFMASK, DOUBLE,2,offset);
01626 
01627   CP_non_std_desc(&inframe, &outcube);
01628 
01629   /* Write type of cube (raw datacube) */
01630   switch (fclass) {
01631   case PRE_OBJ_FRAME:
01632   case FLT_OBJ_FRAME: write_file_class(&outcube, RAW_OBJ_CUBE); break;
01633   case PRE_SKY_FRAME:
01634   case FLT_SKY_FRAME: write_file_class(&outcube, RAW_SKY_CUBE); break;
01635   case PRE_DOM_FRAME: 
01636   case FLT_DOM_FRAME: write_file_class(&outcube, RAW_DOM_CUBE); break;
01637   case PRE_CON_FRAME: 
01638   case FLT_CON_FRAME: write_file_class(&outcube, RAW_CON_CUBE); break;
01639   case PRE_CAL_FRAME: 
01640   case FLT_CAL_FRAME: write_file_class(&outcube, RAW_CAL_CUBE); break;
01641   case SUPER_CLASS:   write_file_class(&outcube, SUPER_CLASS);  break;
01642   }
01643 
01644   close_tiger_frame(&outcube);
01645   close_frame(&inframe);
01646 
01647   free2d(&xofy_coeffs,FLOAT);
01648   free2d(&yofl_coeffs,FLOAT);
01649   free2d(&lofy_coeffs,FLOAT);
01650 
01651   free(nol); free(I); free(J);
01652   free(xdl); free(ydl); 
01653   free(xlo); free(ylo); free(xup); free(yup); 
01654   free(signal); free(variance);
01655   free(lambda); free(lbda_min); free(lbda_max);
01656   free(xofy_deg); free(yofl_deg); free(lofy_deg);
01657   if (WEIGHT) {
01658     for (l=0; l<nlens; l++) {
01659       for (i=0; i<glnggeo; i++) { free(glScs[l][i]); free(glIcs[l][i]); }
01660       free(glScs[l]); free(glIcs[l]);
01661     }
01662     free(glScs); free(glIcs);
01663 
01664     free(glIgeo); free(glXgeo); free(glSgeo);
01665     free2d(&glIpsf,DOUBLE); free2d(&glSpsf,DOUBLE);
01666 
01667     free2d(&sigcoeff,DOUBLE);
01668     free(xlot); free(ylot); free(xupt); free(yupt); 
01669     free(local_lbda_inf_total); free(local_lbda_sup_total); 
01670   }
01671   if (LOCAL) {
01672     free2d(&dxcoeff,DOUBLE);
01673   }
01674 
01675   exit_session(OK);
01676   return(OK);
01677 }
01678 
01679 /* 
01680 ;;; Local Variables: ***
01681 ;;; eval: (add-to-list 'c-font-lock-extra-types "SnifsConfig") ***
01682 ;;; eval: (add-to-list 'c-font-lock-extra-types "SnifsOptics") ***
01683 ;;; End: ***
01684 */

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