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

calib/source/comp_flux.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 #include <gsl/gsl_statistics.h>
00023 #include <gsl/gsl_sf_erf.h>
00024 
00025 #ifdef UNDER
00026 #define e01bef e01bef_
00027 #define e01bff e01bff_
00028 #endif
00029 
00030 #define DBG_PREFIX "dbg_flux"                
00031 
00032 /* ==================================================================== */
00033 /* Miscellaneous routines                                               */
00034 /* ==================================================================== */
00035 
00036 void Interpol_init(int n, double X[], double Y[], double **work)
00037 {
00038   long ifail=1, npts=n;
00039   void e01bef();
00040   
00041   *work = (double *)malloc(n*sizeof(double));
00042 
00043   e01bef(&npts, X, Y, *work, &ifail);
00044 }
00045 
00046 double Interpol(int n, double X[], double Y[], double work[], double lambda)
00047 {
00048   long npts=n, m=1, ifail=0;
00049   double dlambda=lambda, extinction;
00050   void e01bff();
00051 
00052   e01bff(&npts, X, Y, work, &m, &dlambda, &extinction, &ifail);
00053 
00054   return(extinction);
00055 }
00056 
00057 /* === Doxygen Comment ======================================= */
00088 /* =========================================================== */
00089 
00090 int main(int argc,char **argv) 
00091 {
00092   TIGERfile cube;
00093   SPECTRUM Spec, SpecOut;
00094 
00095   char **param, **label, **tabcar1, **tabcar2, **tabcar3, **tabcar4;
00096   char *inname, *outname;
00097   char *ColAName, *ColDName;
00098   char *RefTabName, *ColLbdaName, *ColFluxName;
00099   char *ExcTabName, *ColWavsName, *ColWaveName;
00100   char *ExtTabName, *ColLextName, *ColKextName;
00101   char tmpname[lg_name+1], flxunits[lg_unit +1];
00102 
00103   int *Nolens, *il, i, j, l, k, nlens;
00104   int Nlbda, ncoef, niter, nrej, nexcl, nrefline, nextline, nexcline;
00105   int EXCLUDE, PSF, EXTINCT, SMOOTH;
00106 
00107   float *RealVal, *A, *D, *lbda;
00108   float xc, yc, dx, dy, radius, r2, r2max, exptime, airmass, sig1, sig2, ratio;
00109 
00110   double *Lbda, *Flux, *Data, *ObsVal, *ObsLbda, *Weight, *coef, bounds[4];
00111   double *ExtL, *ExtK, *Kext, *ExclStart, *ExclEnd;
00112   double LbdaMax, LbdaMin, ExtLMin, ExtLMax;
00113   double LbdaStart, LbdaEnd, corr, sigerr;
00114 
00115   /* ====================================================================== */
00116   /* Starting the session                                                   */
00117   /* ====================================================================== */
00118 
00119   set_purpose("Flux calibration (computation)");
00120   set_arglist("-in none -out none -ref none -col A,D -center none -radius none"
00121               " -psf null -smooth null -exclude null -extinct null");
00122   
00123   init_snifs("$Name:  $","$Revision: 1.15 $");
00124   init_session(argv, argc, &label, &param);
00125 
00126   if (DEBUG) {
00127     print_msg("$Id: comp_flux.c,v 1.15 2005/09/15 21:43:14 ycopin Exp $");
00128     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00129   }
00130 
00131   print_warning("+----------------------------------------------------+");
00132   print_warning("| DEPRECATED: use extract_flux (or external program) |");
00133   print_warning("| to extract flux standard spectrum, and calib_flux  |");
00134   print_warning("| to compute flux calibration.                       |");
00135   print_warning("+----------------------------------------------------+");
00136 
00137   set_control_level(NONE);
00138 
00139   /* ====================================================================== */
00140   /* Reading the input parameters                                           */
00141   /* ====================================================================== */
00142 
00143   inname = param[0];                         /* Input datacube */
00144   outname= param[1];                         /* Output flux correction */
00145 
00146   /* Reference flux table */
00147   alloc2d(&tabcar1,3,80,CHAR);
00148   if (decode_argval_char(param[2], tabcar1) != 3) {
00149     print_error("Wrong argument -ref table,col(lbda),col(flux) (got %s)",
00150                 param[2]);
00151     exit_session(ERR_BAD_PARAM);
00152   }
00153   RefTabName  = tabcar1[0];
00154   ColLbdaName = tabcar1[1];
00155   ColFluxName = tabcar1[2];
00156 
00157   /* Input position columns */
00158   alloc2d(&tabcar2,2,80,CHAR);
00159   if (decode_argval_char(param[3], tabcar2) != 2) {
00160     print_error("Wrong argument -col col(X),col(Y) (got %s)",param[3]);
00161     exit_session(ERR_BAD_PARAM);
00162   }
00163   ColAName = tabcar2[0];
00164   ColDName = tabcar2[1];
00165 
00166   /* Summation region (center and radius) */
00167   sscanf(param[4], "%f,%f", &xc, &yc);       /* -center x,y */
00168   get_argval(5, "%f", &radius);              /* -radius r */
00169 
00170   /* PSF estimates */
00171   if ((PSF = is_set(param[6]))) {
00172     if (sscanf(param[6],"%f,%f,%f", &sig1, &sig2, &ratio) != 3) {
00173       print_error("Wrong argument -psf S1,S2,I2/I1 (got %s)",param[6]);
00174       exit_session(ERR_BAD_PARAM);
00175     }
00176     if (sig1 <= 0 || ratio < 0 || sig2 < 0) {
00177       print_error("Wrong PSF parameters (S1>0, S2>=0, I2/I1>=0)");
00178       exit_session(ERR_BAD_PARAM);
00179     }
00180     print_msg("Input PSF estimate: S1=%.2f, S2=%.2f, I2/I1=%.2f",
00181               sig1, sig2, ratio);
00182   }
00183 
00184   /* Smoothing polynom */
00185   if ((SMOOTH = is_set(param[7]))) {
00186     get_argval(7, "%d", &ncoef);
00187     if (ncoef <= 0) {
00188       print_error("Number of smoothing polynomial coefficients must be > 0");
00189       exit_session(ERR_BAD_PARAM);
00190     }
00191   }
00192 
00193   /* Exclude wavelength regions from smoothing */
00194   if ((EXCLUDE = is_set(param[8]))) {
00195     alloc2d(&tabcar3,3,80,CHAR);
00196     if (decode_argval_char(param[8], tabcar3) != 3) {
00197       print_error("Wrong argument -exclude table,col(start),col(end) (got %s)",
00198                   param[8]);
00199       exit_session(ERR_BAD_PARAM);
00200     }
00201     ExcTabName  = tabcar3[0];
00202     ColWavsName = tabcar3[1];
00203     ColWaveName = tabcar3[2];
00204   }
00205 
00206   if (EXCLUDE && !SMOOTH) 
00207     print_warning("Excluded spectral domains are used during smoothing only");
00208 
00209   /* Extinction table */
00210   if ((EXTINCT = is_set(param[9]))) {
00211     alloc2d(&tabcar4,3,80,CHAR);
00212     if (decode_argval_char(param[9], tabcar4) != 3) {
00213       print_error("Wrong argument -extinct table,col(lbda),col(ext) (got %s)",
00214                   param[9]);
00215       exit_session(ERR_BAD_PARAM);
00216     }
00217     ExtTabName  = tabcar4[0];
00218     ColLextName = tabcar4[1];
00219     ColKextName = tabcar4[2];
00220   }
00221 
00222   /* =========================================================================== */
00223   /* Opening the input data cube, checking its class and reading few descriptors */
00224   /* =========================================================================== */
00225 
00226   print_msg("Opening datacube %s", inname);
00227 
00228   if (open_tiger_frame(&cube, inname, "I") < 0) {
00229     print_error("Cannot open datacube %s", inname);
00230     exit_session(ERR_OPEN);
00231   }
00232 
00233   switch (read_file_class(&cube)) {
00234   case SUPER_CLASS:
00235   case FLA_OBJ_CUBE:
00236   case COS_OBJ_CUBE:
00237   case SKY_OBJ_CUBE: break;
00238   default:
00239     print_error("Wrong fclass for %s (should be object datacube)",inname);
00240     exit_session(ERR_BAD_TYPE);
00241   }
00242 
00243   if ((airmass = read_airmass(&cube)) < 1) {
00244     print_error("Cannot read descriptor %s in file %s", AIRMASS, inname);
00245     exit_session(ERR_BAD_DESC);
00246   }
00247 
00248   /* Read exposure time descriptor */
00249   if ((exptime = read_exptime(&cube)) < 0) {
00250     print_error("Cannot read descriptor %s in datacube %s", EFFTIME, inname);
00251     exit_session(ERR_BAD_DESC);
00252   }
00253   print_msg("Integration time: %.1fs",exptime);
00254 
00255   /* ======================================================================= */
00256   /* Compute the wavelength range                                            */
00257   /* 1st case  : if the wavelength range is the same for all spectra, uses   */
00258   /*             directly the output of the get_common_param() routine       */
00259   /* 2nd case  : uses the intersection of the all the wavelength ranges of   */
00260   /*             spectra within the summation region                         */
00261   /* ======================================================================= */
00262 
00263   /* ----------------------------------------------------------------------- */
00264   /* Allocating memory                                                       */
00265   /* ----------------------------------------------------------------------- */
00266 
00267   Nolens = (int *)malloc(cube.nbspec*sizeof(int));
00268   il     = (int *)malloc(cube.nbspec*sizeof(int));
00269   A    = (float *)malloc(cube.nbspec*sizeof(float));
00270   D    = (float *)malloc(cube.nbspec*sizeof(float));
00271 
00272   /* ----------------------------------------------------------------------- */
00273   /* Computing the wavelength domain                                         */
00274   /* ----------------------------------------------------------------------- */
00275 
00276   get_lenses_coord(&cube, ColAName, ColDName, Nolens, A, D, il);
00277 
00278   r2max = radius*radius;
00279 
00280   if (has_common_bounds(&cube)) {
00281     /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00282     /* Case 1 : the spectra have common boundaries                           */
00283     /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00284     get_common_param(&cube, &Nlbda, &LbdaStart, &LbdaEnd);
00285   }
00286   else {                       
00287     /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00288     /* Case 2 : the spectra cab have different lengths, takes the            */
00289     /*          intersection of all the wavelength ranges                    */
00290     /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00291     if (RD_desc(&cube,LRANGE4,DOUBLE,4,bounds) < 0) {
00292       print_error("Cannot read descriptor %s in datacube %s",LRANGE4,cube.name);
00293       exit_session(ERR_NODESC);
00294     }
00295 
00296     LbdaStart = bounds[0];
00297     LbdaEnd = bounds[1];
00298     for (l=0; l<cube.nbspec; l++) {
00299       dx = A[l] - xc;
00300       dy = D[l] - yc;
00301       r2 = dx*dx + dy*dy;
00302       if (r2 > r2max) continue;
00303       get_tiger_spec(&cube, &Spec, NULL, Nolens[l]);
00304       LbdaStart = MAX(LbdaStart,Spec.start);
00305       LbdaEnd   = MIN(LbdaEnd,Spec.end);
00306       free_spec_mem(&Spec);
00307     }
00308     Nlbda = (int)((LbdaEnd - LbdaStart)/cube.step + 1.5);
00309   }
00310 
00311   print_msg("Spectral range: %.1f-%.1f AA", LbdaStart, LbdaEnd);
00312 
00313   /* ----------------------------------------------------------------------- */
00314   /* Allocating and filling the array of the pixel lambda's                  */
00315   /* ----------------------------------------------------------------------- */
00316 
00317   lbda = (float *)malloc(Nlbda*sizeof(float));
00318 
00319   for (i=0; i<Nlbda; i++) lbda[i] = LbdaStart + i*cube.step;
00320 
00321   /* ======================================================================= */
00322   /* Opening and reading the absolute flux table                             */
00323   /* ======================================================================= */
00324 
00325   print_msg("Reading reference table %s", RefTabName);
00326   if ((nrefline = read_fluxref_table(RefTabName, ColLbdaName, ColFluxName,
00327                                      &Lbda, &Flux, flxunits)) < 0) {
00328     exit_session(ERR_BAD_TBL);
00329   }
00330   print_msg("   %d lines read from absolute flux table", nrefline);
00331   print_msg("   Flux units: '%s'",(strlen(flxunits))? flxunits:"None given");
00332 
00333   /* ----------------------------------------------------------------------- */
00334   /* Check the consistency between the wavelength range of the absolute      */
00335   /* flux table and that of our spectra                                      */
00336   /* ----------------------------------------------------------------------- */
00337 
00338   gsl_stats_minmax(&LbdaMin,&LbdaMax,Lbda,1,nrefline);
00339 
00340   /* check that we are within the wavelength range                           */
00341   if ((LbdaStart < LbdaMin) || (LbdaEnd > LbdaMax)) {
00342     print_error("Datacube %s wavelength range [%.2f-%.2f] outside "
00343                 "reference table range [%.2f-%.2f]", 
00344                 inname, LbdaStart, LbdaEnd, LbdaMin, LbdaMax);
00345     exit_session(ERR_BAD_PARAM);
00346   }
00347 
00348   /* ======================================================================= */
00349   /* Summing the spectra within a given radius, the summed spectrum is       */
00350   /* stored in the ObsVal[] array.                                           */
00351   /* ======================================================================= */
00352 
00353   ObsVal = (double *)calloc(Nlbda, sizeof(double));
00354 
00355   /* ----------------------------------------------------------------------- */
00356   /* Summing spectra within the given radius                                 */
00357   /* ----------------------------------------------------------------------- */
00358 
00359   reset_print_progress();
00360   for (nlens=0, l=0; l<cube.nbspec; l++) {
00361     print_progress("Summing spectra within given radius ", (l+1.)/cube.nbspec*100,5.);
00362     dx = A[l] - xc;
00363     dy = D[l] - yc;
00364     r2 = dx*dx + dy*dy;
00365 
00366     if (r2 > r2max) continue;
00367     else if (DEBUG) {
00368       print_msg("%3d: Lens #%3d (%.2f x %.2f) at %.2f",l+1,Nolens[l],A[l],D[l],r2);
00369     }
00370 
00371     get_tiger_spec(&cube, &Spec, NULL, Nolens[l]);
00372 
00373     if (set_subspec(&Spec, LbdaStart, LbdaEnd) != 0) {
00374       print_error("Lens #%d, set_subspec() != 0",Nolens[l]);
00375       exit_session(UNKNOWN);
00376     }
00377 
00378     for (i=0, k=Spec.iwstart; k<=Spec.iwend; i++, k++) 
00379       ObsVal[i] += RD_spec(&Spec, k);
00380 
00381     free_spec_mem(&Spec);
00382     nlens++;
00383   }
00384 
00385   print_msg("-> %d spectra summed within radius %.2f", nlens, radius);
00386 
00387   if (nlens == 0) {
00388     print_error("No valid spectrum was found inside aperture, aborting");
00389     exit_session(ERR_BAD_PARAM);
00390   }
00391 
00392   if (DEBUG) {
00393     sprintf(tmpname,"%s_rawsum",DBG_PREFIX);
00394     print_msg("DEBUG: Saving raw sum spectrum to %s",tmpname);
00395     disable_erase_flag();
00396     create_spec(&SpecOut, tmpname, Nlbda, LbdaStart, cube.step, FLOAT, 
00397                 cube.ident, cube.cunit);
00398     for (i=0; i<Nlbda; i++) WR_spec(&SpecOut, i, ObsVal[i]);
00399     write_file_class(&SpecOut, DEBUG_FILE);
00400     close_spec(&SpecOut);
00401     restore_erase_flag();
00402   }
00403 
00404   /* ======================================================================= */
00405   /* Freeing some arrays which won't be used anymore                         */
00406   /* ======================================================================= */
00407 
00408   free(A); free(D); free(Nolens); free(il);
00409 
00410   /* ======================================================================= */
00411   /* Using the PSF estimate to compute a correction factor on the summed     */
00412   /* flux spectrum (how much flux did we miss by stopping the summation at   */
00413   /* a given radius)                                                         */
00414   /* ======================================================================= */
00415 
00416   if (PSF) {
00417     /* --------------------------------------------------------------------- */
00418     /* Computing the correction                                              */
00419     /* --------------------------------------------------------------------- */
00420     if (ratio != 0) {
00421       corr = (sig2 + ratio*sig1) /
00422         (       sig2*gsl_sf_erf( radius/M_SQRT2/sig1 ) + 
00423           ratio*sig1*gsl_sf_erf( radius/M_SQRT2/sig2) );
00424     }
00425     else {
00426       corr = 1/gsl_sf_erf( radius/M_SQRT2/sig1 );
00427     }
00428 
00429     print_msg("Computed PSF correction factor: 1 + %g", corr-1);
00430 
00431     /* --------------------------------------------------------------------- */
00432     /* Applying the correction                                               */
00433     /* --------------------------------------------------------------------- */
00434     for (i=0; i<Nlbda; i++) ObsVal[i] *= corr;
00435   } 
00436 
00437   /* ======================================================================= */
00438   /* Correcting the observed summed spectrum from the integration time       */
00439   /* (to get a spectrum in ADCU per s) and the airmass (extinction)          */
00440   /* ======================================================================= */
00441 
00442   /* ----------------------------------------------------------------------- */
00443   /* Allocating memory                                                       */
00444   /* ----------------------------------------------------------------------- */
00445 
00446   Kext = (double *)malloc(Nlbda*sizeof(double));
00447 
00448   /* ----------------------------------------------------------------------- */
00449   /* Preparing the extinction spectrum                                       */
00450   /* ----------------------------------------------------------------------- */
00451 
00452   if (EXTINCT) { 
00453     /* ======================================================================= */
00454     /* Read Table with extinction values                                       */
00455     /* ======================================================================= */
00456 
00457     print_msg("Reading extinction table %s (columns %s,%s)", 
00458               ExtTabName, ColLextName, ColKextName);
00459     if ((nextline = read_2columns(ExtTabName, ColLextName, ColKextName, TBL_FILE,
00460                                   &ExtL, &ExtK)) < 0) {
00461       exit_session(ERR_BAD_TBL);
00462     }
00463     print_msg("   %d points in extinction table", nextline);
00464 
00465     /* --------------------------------------------------------------------- */
00466     /* Check the consistency between the wavelength range of the extinction  */
00467     /* table and that of our spectra                                         */
00468     /* --------------------------------------------------------------------- */
00469 
00470     gsl_stats_minmax(&ExtLMin,&ExtLMax,ExtL,1,nextline);
00471 
00472     if ((LbdaStart < ExtLMin) || (LbdaEnd > ExtLMax)) {
00473       print_error("Datacube %s wavelength range [%f-%f] outside "
00474                   "extinction table limits [%f-%f]", 
00475                   inname, LbdaStart, LbdaEnd, ExtLMin, ExtLMax);
00476       exit_session(ERR_BAD_PARAM);
00477     }
00478 
00479     /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00480     /* Interpolating                                                      */
00481     /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00482     Interpol_init(nextline, ExtL, ExtK, &Data);
00483     for (i=0; i<Nlbda; i++) 
00484       Kext[i] = Interpol(nextline, ExtL, ExtK, Data, lbda[i]);
00485     free(ExtL); free(ExtK); free(Data);
00486   }
00487   else {
00488     print_msg("Using default extinction for Mauna-Kea "
00489               "(CFHT Bulletin, number 19, p. 16, 1988)");
00490     for (i=0; i<Nlbda; i++) 
00491       Kext[i] = compute_default_extinction(lbda[i]);
00492   }
00493 
00494   /* ----------------------------------------------------------------------- */
00495   /* Computing the correction spectrum                                       */
00496   /* ----------------------------------------------------------------------- */
00497 
00498   for (i=0; i<Nlbda; i++) 
00499     ObsVal[i] /= ( pow(10.,(-0.4*Kext[i]*airmass)) * exptime );
00500 
00501   free(Kext);
00502 
00503   if (DEBUG) {
00504     sprintf(tmpname,"%s_normsum",DBG_PREFIX);
00505     print_msg("DEBUG: Saving normalized sum spectrum as %s",tmpname);
00506     disable_erase_flag();
00507     create_spec(&SpecOut, tmpname, Nlbda, LbdaStart, cube.step, FLOAT, 
00508                 cube.ident, cube.cunit);
00509     for (i=0; i<Nlbda; i++) WR_spec(&SpecOut, i, ObsVal[i]);
00510     write_file_class(&SpecOut, DEBUG_FILE);
00511     close_spec(&SpecOut);
00512     restore_erase_flag();
00513   }
00514 
00515   /* ======================================================================= */
00516   /* Compute the absolute spectrum from interpolation of the absolute flux   */
00517   /* values stored in the Flux[] array                                       */
00518   /* ======================================================================= */
00519 
00520   RealVal = (float *)malloc(Nlbda*sizeof(float));
00521 
00522   /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00523   /* Interpolating                                                           */
00524   /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00525   Interpol_init(nrefline, Lbda, Flux, &Data);
00526   for (i=0; i<Nlbda; i++) {
00527     RealVal[i] = Interpol(nrefline, Lbda, Flux, Data, lbda[i]);
00528 
00529     if (RealVal[i] <= 0) { /* Paranoid check... */
00530       print_error("ARGH! Reference spectrum <0: Flux(%.2f) = %.2f",
00531                   lbda[i],RealVal[i]);
00532       exit_session(ERR_BAD_PARAM);
00533     }
00534   }
00535   free(Lbda); free(Flux); free(Data);
00536 
00537   if (DEBUG) {
00538     sprintf(tmpname,"%s_ref",DBG_PREFIX);
00539     print_msg("DEBUG: Saving flux reference spectrum as %s",tmpname);
00540     disable_erase_flag();
00541     create_spec(&SpecOut, tmpname, Nlbda, LbdaStart, cube.step, FLOAT, 
00542                 cube.ident, cube.cunit);
00543     for (i=0; i<Nlbda; i++) WR_spec(&SpecOut, i, RealVal[i]);
00544     write_file_class(&SpecOut, DEBUG_FILE);
00545     close_spec(&SpecOut);
00546     restore_erase_flag();
00547   }
00548 
00549   /* ======================================================================= */
00550   /* Compute and save the correction spectrum in ADCU s-1 per flux unit      */
00551   /* ======================================================================= */
00552 
00553   for (i=0; i<Nlbda; i++) ObsVal[i] /= RealVal[i];
00554 
00555   free(RealVal);
00556 
00557   /* ======================================================================= */
00558   /* Filtering the result. Pixels in discarded intervals are assigned        */
00559   /* negligible weights                                                      */
00560   /* ======================================================================= */
00561 
00562   if (SMOOTH) {
00563 
00564     /* ======================================================================= */
00565     /* Read Table with excluded domains                                        */
00566     /* ======================================================================= */
00567 
00568     if (EXCLUDE) { 
00569       print_msg("Reading exclusion table %s", ExcTabName);
00570       if ((nexcline = read_2columns(ExcTabName, ColWavsName, ColWaveName, TBL_FILE,
00571                                     &ExclStart, &ExclEnd)) < 0) {
00572         exit_session(ERR_BAD_TBL);
00573       }
00574       print_msg("   %d wavelength intervals in exclusion table", nexcline);
00575 
00576       /* One should test that ExclStart[i] < ExclEnd[i] */
00577 
00578       if (DEBUG) {
00579         for (i=0; i<nexcline; i++) {
00580           if ((ExclEnd[i] >= LbdaStart) && (ExclStart[i] <= LbdaEnd)) { 
00581             print_msg("   %d: %.1f-%.1f AA", i+1, ExclStart[i], ExclEnd[i]); 
00582           }
00583         }
00584       }
00585       
00586     } /* EXCLUDE */
00587 
00588     print_msg("Polynomial fitting of result spectrum (Ncoef %d)", ncoef);
00589 
00590     /* --------------------------------------------------------------------- */
00591     /* Allocating memory                                                     */
00592     /* --------------------------------------------------------------------- */
00593 
00594     ObsLbda = (double *)malloc(Nlbda*sizeof(double));
00595     Weight  = (double *)malloc(Nlbda*sizeof(double));
00596 
00597     /* --------------------------------------------------------------------- */
00598     /* Filling the ObsLbda[] array and the weight array                      */
00599     /* --------------------------------------------------------------------- */
00600 
00601     for (nexcl=0, i=0; i<Nlbda; i++) {
00602       ObsLbda[i] = lbda[i];
00603       Weight[i]  = 1.;
00604       if (EXCLUDE) {
00605         for (j = 0; j<nexcline; j++) {
00606           if ((ObsLbda[i] >= ExclStart[j]) && (ObsLbda[i] <= ExclEnd[j])) {
00607             Weight[i] = 1.e-99; 
00608             nexcl++;
00609             break; /* if it's in a discarded area, it cannot be in another one */
00610           }
00611         }
00612       }
00613     }
00614 
00615     /* --------------------------------------------------------------------- */
00616     /* DEBUG mode stuff                                                      */
00617     /* --------------------------------------------------------------------- */
00618 
00619     if (DEBUG) {
00620       sprintf(tmpname,"%s_prefit",DBG_PREFIX);
00621       print_msg("DEBUG: Saving pre-fit spectrum as %s",tmpname);
00622       disable_erase_flag();
00623       create_spec(&SpecOut, tmpname, Nlbda, LbdaStart, cube.step, FLOAT, 
00624                   cube.ident, cube.cunit);
00625       for (i=0; i<Nlbda; i++) WR_spec(&SpecOut, i, ObsVal[i]);
00626       write_file_class(&SpecOut, DEBUG_FILE);
00627       close_spec(&SpecOut);
00628       restore_erase_flag();
00629 
00630       if (nexcl > 0) {
00631         print_msg("%d points (over %d) discarded from polynomial fit",
00632                   nexcl, Nlbda);
00633       }
00634       else 
00635         print_msg("No points discarded from polynomial fit");
00636 
00637       sprintf(tmpname,"%s_weight",DBG_PREFIX);
00638       print_msg("DEBUG: Saving weight spectrum as %s",tmpname);
00639       disable_erase_flag();
00640       create_spec(&SpecOut, tmpname, Nlbda, LbdaStart, cube.step, FLOAT, 
00641                   cube.ident, cube.cunit);
00642       for (i=0; i<Nlbda; i++) WR_spec(&SpecOut, i, Weight[i]);
00643       write_file_class(&SpecOut, DEBUG_FILE);
00644       close_spec(&SpecOut);
00645       restore_erase_flag();
00646     }
00647 
00648     /* --------------------------------------------------------------------- */
00649     /* Let's smooth...                                                       */
00650     /* --------------------------------------------------------------------- */
00651     coef = (double *)malloc(ncoef*sizeof(double));
00652 
00653     fit_xpoly_rej_nag_tab(ObsLbda, ObsVal, Weight, 
00654                           Nlbda, ncoef,-3., 3., 20, 
00655                           coef, &niter, &sigerr, &nrej);
00656 
00657     for (i=0; i<Nlbda; i++) {
00658       ObsVal[i] = val_minipoly(ObsLbda[i], coef, ncoef, ObsLbda[0], ObsLbda[Nlbda-1]);
00659     }
00660 
00661     /* --------------------------------------------------------------------- */
00662     /* Free memory                                                           */
00663     /* --------------------------------------------------------------------- */
00664     free(ObsLbda); free(Weight); free(coef);
00665   }
00666 
00667   /* ======================================================================= */
00668   /* Saving the stuff...                                                     */
00669   /* ======================================================================= */
00670 
00671   print_msg("Saving flux correction spectrum in %s", outname);
00672 
00673   create_spec(&SpecOut, outname, Nlbda, LbdaStart, cube.step, FLOAT, 
00674               cube.ident, cube.cunit);
00675   for (i=0; i<Nlbda; i++) WR_spec(&SpecOut, i, ObsVal[i]);
00676 
00677   if (strlen(flxunits)) {
00678     print_msg("Adding flux units keyword %s = '%s' to %s", 
00679               FLXUNITS, flxunits, outname);
00680     WR_desc(&SpecOut, FLXUNITS, CHAR, lg_unit+1, flxunits);
00681   }
00682   else {
00683     print_warning("Please add a %s string-keyword to spectrum %s "
00684                   "to precise flux units.",FLXUNITS,outname);
00685   }
00686 
00687   /* Write type of the throughput spectrum */
00688   write_file_class(&SpecOut, THR_SPEC); 
00689   close_spec(&SpecOut);
00690 
00691   /* ======================================================================= */
00692   /* The end my friend, the end...                                           */
00693   /* ======================================================================= */
00694 
00695   close_tiger_frame(&cube);
00696 
00697   free2d(&tabcar1,CHAR);
00698   free2d(&tabcar2,CHAR); 
00699   if (EXCLUDE) free2d(&tabcar3,CHAR); 
00700   if (EXTINCT) free2d(&tabcar4,CHAR); 
00701 
00702   free(lbda); free(ObsVal);
00703   if (EXCLUDE) {
00704     free(ExclStart);
00705     free(ExclEnd);
00706   }
00707   exit_session(OK);
00708   return(OK);
00709 }

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