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

calib/source/comp_lfff.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00015 /* =========================================================== */
00016 
00017 #include <IFU_io.h>
00018 #include <IFU_math.h>
00019 #include <snifs.h>
00020 #include <calib.h>
00021 
00023 #define NORM_CONT 1
00024 
00025 /* === Doxygen Comment ======================================= */
00034 /* =========================================================== */
00035 
00036 int find_central_indices(TIGERfile *cube, int nmax, int idx[])
00037 {
00038   int *no, *il;
00039   int i,n;
00040   float *a, *d;
00041   double *r;
00042   
00043   no =    (int *)malloc(cube->nbspec*sizeof(int));
00044   il =    (int *)malloc(cube->nbspec*sizeof(int));
00045   a  =  (float *)malloc(cube->nbspec*sizeof(float));
00046   d  =  (float *)malloc(cube->nbspec*sizeof(float));
00047   r  = (double *)malloc(cube->nbspec*sizeof(double)); /* Needed by indexx */
00048   
00049   /* Get lens coords */
00050   get_lenses_coord(cube,"A","D",no,a,d,il);
00051 
00052   /* Compute radii */
00053   for (i=0; i<cube->nbspec; i++) r[i] = hypot(a[i],d[i]);
00054 
00055   /* Sort radii increasingly */
00056   indexx(cube->nbspec,r,il);
00057 
00058   n = MIN(nmax,cube->nbspec);
00059   print_msg("%d central lenses found (%d requested)",n,nmax);
00060   for (i=0; i<n; i++) {
00061     idx[i] = il[i];
00062 
00063     if (DEBUG) 
00064       print_msg("   %3d Lens #%d [%d]:\t%.2f x %.2f (r=%.2f)",
00065                 i+1,cube->signal[idx[i]].specId,idx[i],
00066                 a[idx[i]],d[idx[i]],r[idx[i]]);
00067   }
00068 
00069   free(no); free(il); free(a); free(d); free(r);
00070 
00071   return(n);
00072 }
00073 
00074 /* === Doxygen Comment ======================================= */
00082 /* =========================================================== */
00083 
00084 int find_lens(TIGERfile *frame, int no)
00085 {
00086   int i;
00087 
00088   for (i=0; i<frame->nbspec; i++)
00089     if (frame->signal[i].specId == no) return(i);
00090 
00091   return(-1);
00092 }
00093 
00094 /* === Doxygen Comment ======================================= */
00103 /* =========================================================== */
00104 
00105 float compute_mean_spec(SPECTRUM *spec, float frac, int radius)
00106 {
00107   int i, m, i1, l;
00108   float *work, *fwork, mean;
00109 
00110   /* median filtering of the central part */
00111   m = frac*spec->npts;
00112   work  = (float *)malloc(m*sizeof(float));
00113   fwork = (float *)malloc(m*sizeof(float));
00114   i1 = (1-frac)/2*spec->npts;
00115   for (i=i1, l=0; l<m; l++, i++) work[l] = RD_spec(spec, i);
00116   filter_med_array(work, m, radius, fwork);
00117 
00118   /* compute mean */
00119   for (mean=l=0; l<m; l++) mean += fwork[l];
00120   mean /= m;
00121 
00122   free(work); free(fwork);
00123 
00124   return(mean);
00125 }
00126 
00127 /* === Doxygen Comment ======================================= */
00137 /* =========================================================== */
00138 
00139 void compute_mean_cube(TIGERfile *cube, float frac, int radius, 
00140                        float cont[], int dbglens)
00141 {
00142   SPECTRUM spec;
00143   int k;
00144   
00145   reset_print_progress();
00146   for (k=0; k<cube->nbspec; k++) {
00147     print_progress("Computing mean level", 100.*(k+1)/cube->nbspec,5.);
00148 
00149     /* ------------------------------------------------------- */
00150     /* median window filtering on the central part             */
00151     /* ------------------------------------------------------- */
00152     get_tiger_spec(cube, &spec, NULL, cube->signal[k].specId);
00153     cont[k] = compute_mean_spec(&spec, frac, radius);
00154     free_spec_mem(&spec);
00155 
00156     /* ------------------------------------------------------- */
00157     /* SINGLE LENS DEBUG MODE : print the mean continuum       */
00158     /* value for the debug lens                                */
00159     /* ------------------------------------------------------- */
00160     if (dbglens && cube->signal[k].specId == dbglens) 
00161       print_msg("   Mean level for lens #%d: %f",dbglens, cont[k]);
00162   }
00163   print_msg("   Mean level computed for %d spectra",cube->nbspec); 
00164 }
00165 
00166 /* === Doxygen Comment ======================================= */
00185 /* =========================================================== */
00186 
00187 int main(int argc,char **argv)
00188 {
00189 #if NORM_CONT
00190   SPECTRUM refcontspec;
00191 #endif
00192   TIGERfile *p_refcube;
00193   TIGERfile contcube, flatcube, skycube;
00194   SPECTRUM inspec, noisespec, outspec, skyspec, refskyspec, skynoise;
00195   SPECTRUM tmpspec, tmpnoise;
00196   SPECTRUM dbg_flat, dbg_noiseflat, dbg_nflat, dbg_wflat;
00197 
00198   char **label,**param;
00199   char *contname, *flatname, *skyname;
00200 
00201   int *centrallensk;
00202   int SKY, NOISE, TWILIGHT, SPLINE, i, k, l, m;
00203   int nbref, status;
00204   int Ncentral, Nfound, nolens, reflens;
00205   int singlelens, nbspec, indin, indsky, radius;
00206   int niter, nrej, npts, deg_sky;
00207 
00208   float *mean_sky, *mean_cont;
00209   float frac, fval, AveCenter, norm;
00210 
00211   double *refskycont;
00212   double fspec, fnoise, smooth_param, lbda, error;
00213 
00214 #if !NORM_CONT
00215   int    const ncoef  =15;
00216   int    const itermax=20;
00217   double const siginf =-1.;
00218   double const sigsup = 3.;
00219   double coef[ncoef], lambda, end, start;
00220 #endif
00221   /* SPLINE */
00222   long nknots;
00223   double *coefspline, *xknots, *weight;
00224   double threshold, khi2;
00225   /* !SPLINE */
00226   int    const ncoef_cont  =30;
00227   int    const itermax_cont=30;
00228   double const siginf_cont =-3.;
00229   double const sigsup_cont = 3.;
00230   double coef_cont[ncoef_cont], start_cont, end_cont;
00231 
00232   /* ##### INTRODUCTION ############################## */
00233 
00234   set_purpose("Low-frequency spectro-spatial flat-field (computation)");
00235   set_arglist("-in none -out none -N 10 -smooth 1 -lens null -sky null "
00236               "-frac 0.8 -radius 10 -twilight -reflens|lens_ref 113 -deg_sky 15");
00237 
00238   init_snifs("$Name:  $");
00239   init_session(argv,argc,&label,&param);
00240 
00241   if (DEBUG) {
00242     print_msg("$Id: comp_lfff.c,v 1.9 2004/11/09 17:07:45 ycopin Exp $");
00243     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00244   }
00245 
00246   /* ===== Parameters ============================== */
00247 
00248   /* ----- Input ------------------------------ */
00249 
00250   contname = param[0];
00251   flatname = param[1];
00252   get_argval(2,"%d",&Ncentral);
00253   get_argval(3,"%lf",&smooth_param);
00254   SPLINE = (smooth_param > 0);
00255 
00256   if ((singlelens = is_set(param[4]))) get_argval(4,"%d",&singlelens);
00257   if ((SKY = is_set(param[5]))) skyname = param[5]; /* name of the twilight cube */
00258 
00259   get_argval(6, "%f", &frac);
00260   get_argval(7, "%d", &radius);
00261 
00262   TWILIGHT = is_true(param[8]);
00263   get_argval(9, "%d", &reflens);
00264   get_argval(10, "%d", &deg_sky);
00265 
00266   /* ----- Verifications ------------------------------ */
00267 
00268   if (Ncentral <= 0) {
00269     print_error("N (%d) must be strictly positive.", Ncentral);
00270     exit_session(ERR_BAD_PARAM);
00271   }
00272   if (deg_sky <= 0) {                        /* Check if deg_cal is ok */
00273     print_error("deg_sky (%d) must be strictly positive", deg_sky);
00274     exit_session(ERR_BAD_PARAM);
00275   }
00276   if (TWILIGHT && !SKY) {                    /* -twilight requires -sky */
00277     print_error("When using the -twilight option, "
00278                 "the flat sky expo should be specified");
00279     exit_session(ERR_BAD_PARAM);
00280   }
00281 
00282   /* ##### Let's go ############################## */
00283 
00284   if (SKY) {                         /* Use twilight cube for normalization */
00285 
00286     /* ======================================================== */
00287     /* Opening flat sky datacube, checking if everything is OK  */
00288     /* ======================================================== */
00289 
00290     print_msg("SKY: Opening twilight datacube %s", skyname);
00291     if (open_tiger_frame(&skycube, skyname, "I") < 0) {
00292       print_error("Cannot open frame %s", skyname);
00293       exit_session(ERR_OPEN);
00294     }
00295 
00296     switch (read_file_class(&skycube)) {
00297     case SUPER_CLASS:
00298     case WAV_OBJ_CUBE:
00299     case COS_OBJ_CUBE:
00300     case WAV_DOM_CUBE:
00301     case COS_DOM_CUBE:
00302     case WAV_SKY_CUBE:
00303     case COS_SKY_CUBE: break;
00304     default:
00305       print_error("Wrong file class "
00306                   "(should be lambda-calibrated twilight/dome datacube)");
00307       exit_session(ERR_BAD_TYPE);
00308     }
00309 
00310   } /* SKY */
00311 
00312   if (TWILIGHT) {
00313 
00314     /* ========================================================== */
00315     /* Doing the calculation for the twilight option It will set  */
00316     /* contcube as the normalised datacube coming from skyflat exp. */
00317     /* ========================================================== */
00318 
00319     NOISE = exist_noise(&skycube);
00320 
00321     /* -------------------------------------------------------- */
00322     /* Get the reference sky spectrum                           */
00323     /* -------------------------------------------------------- */
00324     print_msg("Look for ref. spectrum #%d in twilight cube %s",reflens,skycube.name);
00325     if (get_tiger_spec(&skycube, &refskyspec, NULL, reflens) < 0) {
00326       print_error("Reference lens #%d does not exist "
00327                   "in twilight datacube (option -reflens)", reflens);
00328       exit_session(ERR_OPEN);
00329     }
00330 
00331 #if NORM_CONT
00332 
00333     print_msg("TWILIGHT: Using continuum spectrum for renormalization");
00334 
00335     /* -------------------------------------------------------- */
00336     /* Opening continuum data cube and check fclass             */
00337     /* -------------------------------------------------------- */
00338 
00339     print_msg("Opening continuum datacube %s", contname);
00340     if (open_tiger_frame(&contcube,contname,"I") < 0) {
00341       print_error("Cannot open frame %s", contname);
00342       exit_session(ERR_OPEN);
00343     }
00344 
00345     switch (read_file_class(&contcube)) {
00346     case SUPER_CLASS:
00347     case WAV_CON_CUBE:
00348     case COS_CON_CUBE: 
00349     case WAV_DOM_CUBE:
00350     case COS_DOM_CUBE: break;
00351     default:
00352       print_error("Wrong file class for input datacube"
00353                   "(should be lambda-calibrated continuum/dome datacube)");
00354       exit_session(ERR_BAD_TYPE);
00355     }
00356 
00357     /* -------------------------------------------------------- */
00358     /* get the reference continuum spectrum */
00359     /* -------------------------------------------------------- */
00360     print_msg("Extract spectrum #%d from cube %s",reflens,contcube.name);
00361     if (get_tiger_spec(&contcube, &refcontspec, NULL, reflens) < 0) {
00362       print_error("Reference lens #%d does not exist "
00363                   "in continuum datacube (option -reflens)", reflens);
00364       exit_session(ERR_OPEN);
00365     }
00366 
00367     /* -------------------------------------------------------- */
00368     /* Closing the continuum datacube (not needed anymore)      */
00369     /* -------------------------------------------------------- */
00370     close_tiger_frame(&contcube);
00371 
00372     /* -------------------------------------------------------- */
00373     /* Get the continuum spectrum over the twilight range       */
00374     /* -------------------------------------------------------- */
00375     inter_spec(&refskyspec, &refcontspec);   /* Compute common domain */
00376     if ((refskyspec.iwend - refskyspec.iwstart + 1) != refskyspec.npts) {
00377       print_error("Cont. ref. spectrum is not larger than twilight ref. spectrum");
00378       exit_session(ERR_OPEN);
00379     }
00380 
00381     print_msg("Extracting cont. over twilight range [%d-%d/%d]",
00382               refcontspec.iwstart+1,
00383               refcontspec.iwstart+refskyspec.npts,refcontspec.npts);
00384     refskycont = (double *)calloc(refskyspec.npts,sizeof(double));
00385     for (k=refcontspec.iwstart, i=0; i<refskyspec.npts; i++, k++) 
00386       refskycont[i] = RD_spec(&refcontspec, k);
00387 
00388     free_spec_mem(&refcontspec);
00389 
00390 #else
00391 
00392     print_msg("TWILIGHT: Using low-freq. twilight spectrum for renormalization");
00393 
00394     /* -------------------------------------------------------- */
00395     /* Low-freq. polynomial fit                                 */
00396     /* -------------------------------------------------------- */
00397 
00398     start = refskyspec.start * 0.95;
00399     end   = refskyspec.end   * 1.05;
00400     
00401     fit_poly_rej_nag(&refskyspec, ncoef, siginf, sigsup, itermax, 
00402                      coef, &niter, &error, &nrej, start, end, &status);
00403     
00404     refskycont = (double *)malloc(refskyspec.npts * sizeof(double));
00405     for (lambda=refskyspec.start, i=0; i<refskyspec.npts; 
00406          i++, lambda += refskyspec.step) 
00407       refskycont[i] = val_poly_nag(lambda, coef, ncoef, start, end, &status);
00408 
00409 #endif
00410 
00411     /* =========================================================== */
00412     /* Create the temporary datacube which will serve as input     */
00413     /* for the calculation of the correction flat-field            */
00414     /* Here we will use the resulting calculated datacube          */
00415     /* =========================================================== */
00416 
00417     print_msg("TWILIGHT: creating temp. datacube with renormalized twilight");
00418 
00419     if (create_tiger_frame(&contcube, "dbg_comp_lfff", -1, 0., skycube.step, 
00420                            skycube.data_type, skycube.table_name, 
00421                            "Input for comp_lfff", NULL) < 0) {
00422       print_error("Cannot create temporary datacube dbg_comp_lfff");
00423       exit_session(ERR_CREAT);
00424     }
00425 
00426     /* -------------------------------------------------------- */
00427     /* We transfer everything into the temp datacube            */
00428     /* -------------------------------------------------------- */
00429     CP_non_std_desc(&skycube, &contcube);
00430     for (nrej=i=0; i<skycube.nbspec; i++) {
00431       if (NOISE) get_tiger_spec(&skycube,&skyspec,&skynoise,skycube.signal[i].specId);
00432       else       get_tiger_spec(&skycube,&skyspec,NULL,     skycube.signal[i].specId);
00433 
00434       /* -------------------------------------------------------- */
00435       /* For each point, we normalise by the ref spectrum itself  */
00436       /* continuum divided                                        */
00437       /* -------------------------------------------------------- */
00438 
00439       inter_spec(&skyspec, &refskyspec);
00440       npts = skyspec.iwend - skyspec.iwstart + 1;
00441       if (npts < 0) {
00442         print_warning("   Lens #%d discarded (no common wavelength range)",
00443                       skycube.signal[i].specId);
00444         nrej++;
00445         free_spec_mem(&skyspec);
00446         if (NOISE) free_spec_mem(&skynoise);
00447         continue;
00448       }
00449 
00450       init_new_tiger_spec(&contcube, &tmpspec, npts, skyspec.wstart);
00451       if (NOISE) init_new_tiger_spec(&contcube, &tmpnoise, npts, skyspec.wstart);
00452       
00453       for (k=0, l=skyspec.iwstart, m=refskyspec.iwstart;k < npts;k++, l++, m++) {
00454         WR_spec(&tmpspec, k, 
00455                 (double)refskycont[m] * RD_spec(&skyspec,l)/RD_spec(&refskyspec,m));
00456         if (NOISE) WR_spec(&tmpnoise, k, 1.);
00457       }
00458       if (NOISE) put_tiger_spec(&contcube,&tmpspec,&tmpnoise,skycube.signal[i].specId);
00459       else       put_tiger_spec(&contcube,&tmpspec,NULL,     skycube.signal[i].specId);
00460 
00461       free_spec_mem(&skyspec);
00462       if (NOISE) free_spec_mem(&skynoise);
00463     }
00464     print_msg("   -> %d spectra rejected from twilight cube", nrej);
00465 
00466     free(refskycont);
00467     free_spec_mem(&refskyspec);
00468 
00469   } /* TWILIGHT */
00470 
00471   else {
00472 
00473     /* ========================================================== */
00474     /* Opening continuum datacube, checking if everything is OK   */
00475     /* ========================================================== */
00476 
00477     print_msg("Opening continuum datacube %s", contname);
00478     if (open_tiger_frame(&contcube,contname,"I") < 0) {
00479       print_error("Cannot open cube %s", contname);
00480       exit_session(ERR_OPEN);
00481     }
00482 
00483     switch (read_file_class(&contcube)) {
00484     case SUPER_CLASS:
00485     case WAV_DOM_CUBE:
00486     case COS_DOM_CUBE:
00487     case WAV_CON_CUBE:
00488     case COS_CON_CUBE: break;
00489     default:
00490       print_error("Wrong file class "
00491                   "(should be lambda-calibrated continuum/dome datacube)");
00492       exit_session(ERR_BAD_TYPE);
00493     }
00494 
00495     NOISE = exist_noise(&contcube);
00496 
00497     if (SKY && skycube.nbspec != contcube.nbspec) 
00498       print_warning("Sky and continuum cubes have different "
00499                     "number of spectra: %d vs. %d",skycube.nbspec,contcube.nbspec);
00500 
00501   } /* !TWILIGHT */
00502 
00503   /* =========================================================== */
00504   /* Mode debug sur un seul spectre                              */
00505   /* Checking that the lens and the corresponding spectrum exist */
00506   /* =========================================================== */
00507 
00508   if (singlelens) {
00509     /* ------------------------------------------------------- */
00510     /* Check for the existence of this lens in [each] datacube */
00511     /* ------------------------------------------------------- */
00512     if (exist_lens(&contcube, singlelens) < 0) {
00513       print_error("Lens %d not found in datacube %s", singlelens, contcube.name);
00514       exit_session(ERR_BAD_PARAM);
00515     }
00516     if (SKY && exist_lens(&skycube, singlelens) < 0) {
00517       print_error("Lens %d not found in sky datacube %s", singlelens, skycube.name);
00518       exit_session(ERR_BAD_PARAM);
00519     }
00520 
00521     nbspec = 1;
00522     nolens = singlelens;
00523   } 
00524   else nbspec = contcube.nbspec;
00525 
00526   /* ======================================================== */
00527   /* Compute the mean of each continuum spectrum              */
00528   /* and storing the computed value in the mean_cont[] array  */
00529   /* ======================================================== */
00530 
00531   print_msg("Computing mean continuum level");
00532   mean_cont = (float *)malloc(contcube.nbspec*sizeof(float));
00533   compute_mean_cube(&contcube, frac, radius, mean_cont, singlelens);
00534 
00535   if (SKY) {
00536     /* ======================================================== */
00537     /* Storing the mean value of each sky spectrum into the     */
00538     /* mean_sky[] array.                                        */
00539     /* ======================================================== */
00540 
00541     print_msg("Computing mean twilight level");
00542     mean_sky = (float *)malloc(skycube.nbspec*sizeof(float));
00543     compute_mean_cube(&skycube, frac, radius, mean_sky, singlelens);
00544   } 
00545 
00546   /* ============================================================ */
00547   /* Looking for Ncentral lenses and storing their index into     */
00548   /* centrallensk[]. Using either the SKY or CONTINUUM datacubes. */
00549   /* ============================================================ */
00550 
00551   centrallensk = (int *)malloc(Ncentral * sizeof(int));
00552 
00553   /* Reference datacube */
00554   if (SKY) {
00555     print_msg("Sky datacube used as reference datacube");
00556     p_refcube = &skycube;
00557   } 
00558   else {
00559     print_msg("Continuum datacube used as reference datacube");
00560     p_refcube = &contcube;
00561   }
00562   nbref = p_refcube->nbspec;
00563 
00564   /* Searching for central lens indices */
00565 
00566   if (Ncentral > nbref) {
00567     print_error("N (%d) is larger than the total number of"
00568                 " spectra in the reference datacube (%d)", Ncentral,nbref);
00569     exit_session(ERR_BAD_PARAM);
00570   }
00571 
00572   Nfound = find_central_indices(p_refcube, Ncentral, centrallensk);
00573 
00574   /* ================================================================== */
00575   /* Compute the normalization factor using the selected central lenses */
00576   /* ================================================================== */
00577 
00578   AveCenter = 0;
00579   if (SKY) for (i=0; i<Nfound; i++) AveCenter += mean_sky[centrallensk[i]];
00580   else     for (i=0; i<Nfound; i++) AveCenter += mean_cont[centrallensk[i]];
00581   AveCenter /= Nfound;
00582   print_msg("-> Mean reference value from %d central spectra: %f", 
00583             Nfound, AveCenter);
00584 
00585   free(centrallensk);
00586 
00587   /* ============================================================== */
00588   /* Open the output data cube (if we are not in the single lens    */
00589   /* debug mode)                                                    */
00590   /* ============================================================== */
00591 
00592   if (!singlelens) {
00593     print_msg("Creating output datacube %s", flatname);
00594     if (create_tiger_frame(&flatcube, flatname, -1, 0., contcube.step, contcube.data_type, 
00595                            contcube.table_name, contcube.ident, contcube.cunit) < 0) {
00596       print_error("Cannot create output datacube");
00597       exit_session(ERR_CREAT);
00598     }
00599   }
00600 
00601   /* ============================================================== */
00602   /* Starting the serious stuff of computing the low-frequency flat */
00603   /* spectra                                                        */
00604   /* ============================================================== */
00605 
00606   norm = AveCenter; /* will be changed if SKY */
00607 
00608   /* --------------------------------------------------------------- */
00609   /* Starting the loop on indices of lenses (contcube.nbspec or 1)     */
00610   /* --------------------------------------------------------------- */
00611 
00612   reset_print_progress();
00613   for (indin=0; indin<nbspec; indin++) {
00614     /* --------------------------------------------------------------------- */
00615     /* Get nolens from indin in [each] datacube */
00616     /* (or indin from nolens if singlelens) */
00617     /* --------------------------------------------------------------------- */
00618     if (singlelens) indin = find_lens(&contcube,nolens);
00619     else {
00620       print_progress("Computing low frequency flat-field", 
00621                      100.*(indin+1)/nbspec, 1.);
00622       nolens = contcube.signal[indin].specId;
00623     }
00624 
00625     if (SKY) {
00626       indsky = find_lens(&skycube,nolens);
00627 
00628       if (indsky < 0) {
00629         print_warning("Lens #%d [%d] in continuum datacube does not exist"
00630                       " in sky datacube, skipping...", nolens, indin);
00631         continue; /* cannot happen if singlelens, because already tested */
00632       } 
00633       else if (mean_sky[indsky] <= 0) {
00634         print_warning("Negative or null mean value of the sky spectrum"
00635                       "for lens # %d [%d], skipping...", nolens, indsky);
00636         if (singlelens) {
00637           print_error("Results for one lens not applicable...");
00638           norm = 1.;
00639         } 
00640         else continue;
00641       } 
00642       else {
00643         /* ------------------------------------------------------------ */
00644         /* Compute the normalisation coefficient                        */
00645         /* AveCenter(continuum)  if no sky datacube is available        */
00646         /* AveCenter(sky) * mean_cont / mean_sky otherwise              */ 
00647         /* ------------------------------------------------------------ */
00648         norm = AveCenter * mean_cont[indin] / mean_sky[indsky];
00649 
00650         if (norm <= 0) {
00651           print_error("ARGH!!! Norm = %f <= 0 for lens #%d", norm, nolens);
00652           norm = 1;
00653         }
00654       }
00655     }
00656 
00657     /* ------------------------------------------------------------ */
00658     /* Open the continuum spectrum (and its corresponding noise     */
00659     /* spectrum if present)                                         */
00660     /* ------------------------------------------------------------ */
00661 
00662     if (NOISE) get_tiger_spec(&contcube, &inspec, &noisespec, nolens);
00663     else       get_tiger_spec(&contcube, &inspec, NULL,       nolens);
00664 
00665     /* ------------------------------------------------------------ */
00666     /* Open the output spectrum or the debug spectra in the case    */
00667     /* of the single lense mode                                     */
00668     /* ------------------------------------------------------------ */
00669 
00670     if (singlelens) {
00671       print_msg("Normalisation coefficient for lens #%d [%d]: %f",
00672                 nolens, indin, norm);
00673 
00674       print_msg("DEBUG - Creating spectra:");
00675       disable_erase_flag();
00676       print_msg("-> dbg_flat: continuum spectrum");
00677       create_spec(&dbg_flat, "dbg_flat", inspec.npts, inspec.start, inspec.step,
00678                   inspec.data_type, inspec.ident, inspec.cunit);
00679       print_msg("-> dbg_nflat: normalized continuum spectrum");
00680       create_spec(&dbg_nflat, "dbg_nflat", inspec.npts, inspec.start, inspec.step,
00681                   inspec.data_type, inspec.ident, inspec.cunit);
00682       print_msg("-> dbg_fitflat: fit to the normalized continuum spectrum");
00683       create_spec(&outspec, "dbg_fitflat", inspec.npts, inspec.start, inspec.step,
00684                   inspec.data_type, inspec.ident, inspec.cunit);
00685       if (SPLINE && NOISE) {
00686         print_msg("-> dbg_noiseflat: continuum noise spectrum");
00687         create_spec(&dbg_noiseflat, "dbg_noiseflat", inspec.npts, inspec.start, 
00688                     inspec.step,inspec.data_type, inspec.ident, inspec.cunit);
00689         print_msg("-> dbg_wflat: weights used in the fit");
00690         create_spec(&dbg_wflat, "dbg_wflat", inspec.npts, inspec.start, inspec.step,
00691                     inspec.data_type, inspec.ident, inspec.cunit);
00692       }
00693       restore_erase_flag();
00694     }
00695     else init_new_tiger_spec(&flatcube, &outspec, inspec.npts, inspec.start);
00696 
00697     /* ------------------------------------------------------------ */
00698     /* Temporarly (it will be overwritten later) store the          */
00699     /* normalized continuum spectrum in outspec. Computes the       */
00700     /* weights using the noise spectrum                             */
00701     /* ------------------------------------------------------------ */
00702 
00703     for (i=0; i<inspec.npts; i++) {
00704       fspec = RD_spec(&inspec, i);
00705       WR_spec(&outspec, i, fspec/norm);
00706 
00707       if (singlelens) {
00708         WR_spec(&dbg_flat, i, fspec);
00709         WR_spec(&dbg_nflat, i, fspec/norm);
00710       }
00711     }
00712 
00713     /* ------------------------------------------------------------ */
00714     /* Fit the continuum spectrum                                   */
00715     /* ------------------------------------------------------------ */
00716 
00717     if (SPLINE) {                            /* Spline */
00718 
00719       /* In the discussion of E02BEF, the threshold is suggested in the range
00720          Variance(w*y)*(npts+-sqrt(npts)), i.e.
00721          SQ(weight[0])*fnoise*(inspec.npts +- sqrt(2*inspec.npts)).  However,
00722          this doesn't work at all!!! (by a factor ~10^9) */
00723 
00724       threshold = smooth_param*inspec.npts;
00725 
00726       coefspline= (double *)malloc((inspec.npts+4)*sizeof(double));
00727       xknots    = (double *)malloc((inspec.npts+4)*sizeof(double));
00728       weight    = (double *)malloc(    inspec.npts*sizeof(double));
00729 
00730       for (i=0; i<inspec.npts; i++) {
00731         if (NOISE) {
00732           fnoise = RD_spec(&noisespec, i);
00733           if (fnoise<=0) {
00734             print_error("ARGH!!! noise = %f <= 0 for lens #%d, pix. %d",
00735                         fnoise, nolens, i);
00736             fnoise = 1;
00737           }          
00738           weight[i] = norm/sqrt(fnoise);     /* weight = 1/normed_sigma */
00739         
00740           if (singlelens) {
00741             WR_spec(&dbg_noiseflat, i, fnoise);
00742             WR_spec(&dbg_wflat, i, weight[i]);
00743           }
00744         } 
00745         else weight[i] = norm;               /* Previously to sqrt(norm) (why???) */
00746       }
00747 
00748       status = fit_spline_nag(&outspec, threshold, weight, coefspline, 
00749                               &nknots, xknots, &khi2);
00750 
00751     }
00752     else {                                   /* Polynomial fit */
00753 
00754       start_cont = outspec.start * 0.95;
00755       end_cont   = outspec.end   * 1.05;
00756     
00757       fit_poly_rej_nag(&outspec, ncoef_cont, siginf_cont, sigsup_cont, 
00758                        itermax_cont, coef_cont, &niter, &error, &nrej, 
00759                        start_cont, end_cont, &status);
00760 
00761     }
00762     
00763     if (!status) {
00764       /* ------------------------------------------------------------ */
00765       /* Write the output spectrum (fitted continuum spectrum)        */
00766       /* ------------------------------------------------------------ */
00767 
00768       for (i = 0; i < inspec.npts; i++) {
00769         lbda = coord_spec(&outspec, i);
00770         if (SPLINE) fval = val_spline_nag(lbda, coefspline, nknots, xknots);
00771         else        fval = val_poly_nag(lbda, coef_cont, ncoef_cont, 
00772                                         start_cont, end_cont, &status);
00773         WR_spec(&outspec, i, fval);
00774       }
00775 
00776       /* ------------------------------------------------------------ */
00777       /* Close all the stuff before looping                           */
00778       /* ------------------------------------------------------------ */
00779 
00780       if (!singlelens) put_tiger_spec(&flatcube, &outspec, NULL, nolens); 
00781     } 
00782     else
00783       print_warning("Error %d in fit procedure, spectrum %d discarded",
00784                     status, nolens);
00785 
00786     free_spec_mem(&inspec);
00787     if (NOISE) free_spec_mem(&noisespec);
00788     if (SPLINE) { free(coefspline); free(xknots); free(weight); }
00789   }
00790 
00791   /* =============================================================== */
00792   /* This is the end my friend the end                               */
00793   /* =============================================================== */
00794 
00795   free(mean_cont);
00796   if (SKY) free(mean_sky);
00797   
00798   if (singlelens) {
00799     write_file_class(&outspec, DEBUG_FILE); 
00800     close_spec(&outspec);
00801     write_file_class(&dbg_flat, DEBUG_FILE); 
00802     close_spec(&dbg_flat);
00803     write_file_class(&dbg_nflat, DEBUG_FILE); 
00804     close_spec(&dbg_nflat);
00805     if (NOISE) {
00806       write_file_class(&dbg_noiseflat, DEBUG_FILE); 
00807       close_spec(&dbg_noiseflat);
00808       write_file_class(&dbg_wflat, DEBUG_FILE); 
00809       close_spec(&dbg_wflat);
00810     } 
00811   } 
00812   else {
00813     CP_non_std_desc(&contcube, &flatcube);
00814 
00815     /* write type of cube (Low freq. flat field)) */
00816     write_file_class(&flatcube, FLAT_CUBE); 
00817     close_tiger_frame(&flatcube);
00818     close_tiger_frame(&contcube);
00819   }
00820 
00821   exit_session(OK);
00822   return(OK);
00823 }
00824 

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