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

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