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

extract/source/new_extract_spec.c

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

Generated on Mon Sep 19 15:19:24 2005 for Snifs by doxygen 1.3.5