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

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