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

calib/source/calib_flux.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 #include "../../extract/incl/extract.h" /* Probably not the right way... */
00022 #include <gsl/gsl_statistics.h>
00023 #include <gsl/gsl_spline.h>
00024 
00025 #define DBG_PREFIX "dbg_flux"                
00026 #define GAUSSNSIG  3
00027 
00028 int Interpolate(int norig, double xorig[], double yorig[], 
00029                 int ntarg, double xtarg[], double ytarg[],
00030                 const gsl_interp_type *interp_type)
00031 {
00032   bool ordered = 1;
00033   int i;
00034   gsl_interp_accel *acc = gsl_interp_accel_alloc();
00035   gsl_spline *spline = gsl_spline_alloc(interp_type, norig);
00036 
00037   if (DEBUG) {                  /* Test ordering of input array */
00038     for (i=1; i<norig; i++) {
00039       if (xorig[i] <= xorig[i-1]) {
00040         ordered = 0;
00041         break;
00042       }
00043     }
00044     if (!ordered) {
00045       print_error("Input array should be strictly increasing in Interpolate "
00046                   "(x[%d]=%f < x[%d]=%f)",i,xorig[i],i-1,xorig[i-1]);
00047       return(UNKNOWN);
00048     }
00049   }
00050   
00051   gsl_spline_init(spline, xorig, yorig, norig);
00052   for (i=0; i<ntarg; i++) 
00053     ytarg[i] = gsl_spline_eval(spline, xtarg[i], acc);
00054 
00055   gsl_spline_free(spline);
00056   gsl_interp_accel_free(acc);
00057   
00058   return(OK);
00059 }
00060 
00061 int fillHalfGauss(double halfGauss[], int npts, double sigma) 
00062 {
00063   int i, nmax;
00064   double fact, norm;
00065 
00066   /* Max extension of the gaussian */
00067   nmax = ceil(GAUSSNSIG*sigma);
00068 
00069   /* Compute unnormalized gaussian on 1st nmax pixels */
00070   fact = 1./(2*SQ(sigma));
00071   halfGauss[0] = norm = 1;
00072   for (i=1; i<nmax; i++) {
00073     halfGauss[i] = exp(-SQ(i)*fact);
00074     norm += 2*halfGauss[i];
00075   }
00076   /* Clean the rest of the gauss array */
00077   for (i=nmax; i<npts; i++) halfGauss[i] = 0;
00078   /* Numerically renormalize gaussian (sum[-nmax+1..nmax-1] = 1) */
00079   for (i=0; i<nmax; i++) halfGauss[i] /= norm;
00080 
00081   return(nmax);
00082 }
00083 
00084 void degradeSpecRes(double array[], int npts, double sigma[])
00085 {
00086   int i,j, hgNpts, nmax;
00087   double *copy, *halfGauss;
00088   double sigMax;
00089 
00090   copy = (double *)calloc(npts, sizeof(double));
00091   
00092   sigMax = gsl_stats_max(sigma, 1, npts);
00093   hgNpts = ceil(GAUSSNSIG*sigMax);
00094   print_msg("Max. sigma = %.2f px -> allocating %d px for halfGauss",
00095             sigMax, hgNpts);
00096   halfGauss = (double *)calloc(hgNpts, sizeof(double));
00097   
00098   for (i=0; i<npts; i++) {
00099     /* Fill in the Gaussian profile. 
00100     *  TODO: skip this step if the sigma is the same as before. */
00101     nmax = fillHalfGauss(halfGauss, hgNpts, sigma[i]);
00102     /* Bestial convolution */
00103     copy[i] += array[i]*halfGauss[0];
00104     for (j=1; j<nmax; j++) {
00105       /* For convolution, extend array on each side by the
00106          first/last-pixel value. */
00107       copy[i] += halfGauss[j]*array[MAX(i-j,0)];
00108       copy[i] += halfGauss[j]*array[MIN(i+j,npts-1)];
00109     }
00110   }
00111 
00112   /* Store result in input array */
00113   memcpy(array, copy, npts*sizeof(double));
00114 
00115   free(copy); free(halfGauss);
00116 }
00117 
00118 
00119 /* === Doxygen Comment ======================================= */
00145 /* =========================================================== */
00146 
00147 int main(int argc,char **argv) 
00148 {
00149   SPECTRUM inspec, outspec;
00150 
00151   bool EXCLUDE, EXTINCT, SMOOTH, MATCHRES, MATCHBIN;
00152 
00153   char **label, **param, **tabcar1, **tabcar3, **tabcar4;
00154   char *inname, *outname;
00155   char *RefTabName, *ColLbdaName, *ColFluxName;
00156   char *ExcTabName, *ColWavsName, *ColWaveName;
00157   char *ExtTabName, *ColLextName, *ColKextName;
00158   char dbgname[lg_name+1], refUnits[lg_unit +1];
00159 
00160   int i, j;
00161   int obsNpts, ncoef, niter, nrej, obsNExcl, refNpts, extNpts, exclNpts;
00162 
00163   float exptime, airmass, specres;
00164 
00165   double *refLbda, *refFlux, *refRes, *obsRef;
00166   double *obsLbda, *obsVal, *obsK, *obsDeg, *wLbda, *weight, *coef;
00167   double *extLbda, *extK, *exclStart, *exclEnd;
00168   double refLMax, refLMin, extLMin, extLMax, obsLMin, obsLMax, sigerr;
00169 
00170   /* ====================================================================== */
00171   /* Starting the session                                                   */
00172   /* ====================================================================== */
00173 
00174   set_purpose("Flux calibration (computation)");
00175   set_arglist("-in none -out none -ref none -smooth null "
00176               "-exclude null -extinct null -matchres -matchbin");
00177   
00178   init_snifs("$Name:  $","$Revision: 1.4 $");
00179   init_session(argv, argc, &label, &param);
00180 
00181   if (DEBUG) {
00182     print_msg("$Id: calib_flux.c,v 1.4 2005/10/05 22:44:30 ycopin Exp $");
00183     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00184   }
00185 
00186   /* ====================================================================== */
00187   /* Reading the input parameters                                           */
00188   /* ====================================================================== */
00189 
00190   inname  = param[0];                        /* Input standard spectrum */
00191   outname = param[1];                        /* Output flux correction */
00192 
00193   /* Reference flux table */
00194   alloc2d(&tabcar1,3,lg_name+1,CHAR);
00195   if (decode_argval_char(param[2], tabcar1) != 3) {
00196     print_error("Wrong argument -ref table,col(lbda),col(flux) (got %s)",
00197                 param[2]);
00198     exit_session(ERR_BAD_PARAM);
00199   }
00200   RefTabName  = tabcar1[0];
00201   ColLbdaName = tabcar1[1];
00202   ColFluxName = tabcar1[2];
00203 
00204   /* Smoothing polynom */
00205   if ((SMOOTH = is_set(param[3]))) {
00206     get_argval(3, "%d", &ncoef);
00207     if (ncoef <= 0) {
00208       print_error("Number of coeff. of smoothing polynomial must be > 0");
00209       exit_session(ERR_BAD_PARAM);
00210     }
00211   }
00212 
00213   /* Exclude wavelength regions from smoothing */
00214   if ((EXCLUDE = is_set(param[4]))) {
00215     alloc2d(&tabcar3,3,lg_name+1,CHAR);
00216     if (decode_argval_char(param[4], tabcar3) != 3) {
00217       print_error("Wrong argument -exclude table,col(start),col(end) (got %s)",
00218                   param[4]);
00219       exit_session(ERR_BAD_PARAM);
00220     }
00221     ExcTabName  = tabcar3[0];
00222     ColWavsName = tabcar3[1];
00223     ColWaveName = tabcar3[2];
00224   }
00225 
00226   if (EXCLUDE && !SMOOTH) 
00227     print_warning("Excluded spectral domains are used during smoothing only");
00228 
00229   /* Extinction table */
00230   if ((EXTINCT = is_set(param[5]))) {
00231     alloc2d(&tabcar4,3,lg_name+1,CHAR);
00232     if (decode_argval_char(param[5], tabcar4) != 3) {
00233       print_error("Wrong argument -extinct table,col(lbda),col(ext) (got %s)",
00234                   param[5]);
00235       exit_session(ERR_BAD_PARAM);
00236     }
00237     ExtTabName  = tabcar4[0];
00238     ColLextName = tabcar4[1];
00239     ColKextName = tabcar4[2];
00240   }
00241 
00242   MATCHRES = is_true(param[6]);  /* Match spectral resolution */
00243   MATCHBIN = is_true(param[7]);  /* Match spectral bins */
00244 
00245   if (MATCHBIN) {
00246     print_error("Option -matchbin not yet implemented");
00247     exit_session(ERR_NOIMPL);
00248   }
00249 
00250   /* ======================================================================= */
00251   /* Open the input standard spectrum and read few descriptors               */
00252   /* ======================================================================= */
00253 
00254   print_msg("Opening flux standard spectrum %s", inname);
00255   if (open_spec(&inspec, inname,"I") < 0) {
00256     print_error("Cannot open spectrum %s",inname);
00257     exit_session(ERR_OPEN);
00258   }
00259 
00260   /* Read airmass descriptor */
00261   if ((airmass = read_airmass(&inspec)) < 1) {
00262     print_error("Cannot read descriptor %s in %s", AIRMASS, inname);
00263     exit_session(ERR_BAD_DESC);
00264   }
00265   /* Read exposure time descriptor */
00266   if ((exptime = read_exptime(&inspec)) < 0) {
00267     print_error("Cannot read descriptor %s in %s", EFFTIME, inname);
00268     exit_session(ERR_BAD_DESC);
00269   }
00270   print_msg("   Integration time: %.1fs, Airmass: %.2f", exptime, airmass);
00271 
00272   obsLMin = inspec.start;
00273   obsLMax = inspec.end;
00274   obsNpts = inspec.npts;
00275   print_msg("   Spectral range: %.2f-%.2f AA", obsLMin, obsLMax);
00276 
00277   obsLbda = (double *)malloc(obsNpts*sizeof(double));
00278   obsVal  = (double *)malloc(obsNpts*sizeof(double));
00279 
00280   for (i=0; i<obsNpts; i++) {
00281     obsLbda[i]= obsLMin + i*inspec.step;
00282     obsVal[i] = RD_spec(&inspec, i);
00283   }
00284 
00285   if (MATCHRES) {
00286     if (RD_desc(&inspec, SPECRES, FLOAT, 1, &specres) < 0) {
00287       print_error("Cannot read descriptor %s in %s", SPECRES, inname);
00288       exit_session(ERR_BAD_DESC);
00289     }
00290     print_msg("   Spectral resolution: %.3f px (%.3f AA)", 
00291               specres, specres*inspec.step);
00292   }
00293 
00294   close_spec(&inspec);
00295 
00296   /* ======================================================================= */
00297   /* Opening and reading the absolute flux table                             */
00298   /* ======================================================================= */
00299 
00300   print_msg("Reading reference table %s", RefTabName);
00301   if ((refNpts = read_fluxref_table(RefTabName, ColLbdaName, ColFluxName,
00302                                     &refLbda, &refFlux, refUnits,
00303                                     MATCHRES?"RESOLUTION":NULL,&refRes)) < 0) {
00304     exit_session(ERR_BAD_TBL);
00305   }
00306   print_msg("   %d lines read from absolute flux table", refNpts);
00307   print_msg("   Flux units: '%s'",(strlen(refUnits))? refUnits:"None given");
00308 
00309   /* ----------------------------------------------------------------------- */
00310   /* Check the consistency between the wavelength range of the absolute      */
00311   /* flux table and that of our spectra                                      */
00312   /* ----------------------------------------------------------------------- */
00313 
00314   gsl_stats_minmax(&refLMin,&refLMax,refLbda,1,refNpts);
00315 
00316   /* check that we are within the wavelength range                           */
00317   if ((refLMin > obsLMin) || (refLMax < obsLMax)) {
00318     print_error("Flux std spectrum %s wavelength range [%.2f-%.2f] outside "
00319                 "reference table range [%.2f-%.2f]", 
00320                 inname, obsLMin, obsLMax, refLMin, refLMax);
00321     exit_session(ERR_BAD_PARAM);
00322   }
00323 
00324   /* ======================================================================= */
00325   /* Correcting the observed summed spectrum from the integration time       */
00326   /* (to get a spectrum in ADCU per s) and the airmass (extinction)          */
00327   /* ======================================================================= */
00328 
00329   obsK = (double *)malloc(obsNpts*sizeof(double));
00330 
00331   /* ----------------------------------------------------------------------- */
00332   /* Preparing the extinction spectrum                                       */
00333   /* ----------------------------------------------------------------------- */
00334 
00335   if (EXTINCT) { 
00336 
00337     /* ===================================================================== */
00338     /* Read Table with extinction values                                     */
00339     /* ===================================================================== */
00340 
00341     print_msg("Reading extinction table %s (columns %s,%s)", 
00342               ExtTabName, ColLextName, ColKextName);
00343     if ((extNpts = read_2columns(ExtTabName, ColLextName, ColKextName, 
00344                                   TBL_FILE, &extLbda, &extK)) < 0) {
00345       exit_session(ERR_BAD_TBL);
00346     }
00347     print_msg("   %d points in extinction table", extNpts);
00348 
00349     /* --------------------------------------------------------------------- */
00350     /* Check the consistency between the wavelength range of the extinction  */
00351     /* table and that of our spectra                                         */
00352     /* --------------------------------------------------------------------- */
00353     gsl_stats_minmax(&extLMin,&extLMax,extLbda,1,extNpts);
00354     if ((extLMin > obsLMin) || (extLMax < obsLMax)) {
00355       print_error("Flux std spectrum %s wavelength range [%f-%f] outside "
00356                   "extinction table limits [%f-%f]", 
00357                   inname, obsLMin, obsLMax, extLMin, extLMax);
00358       exit_session(ERR_BAD_PARAM);
00359     }
00360 
00361     /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00362     /* Interpolation                                                      */
00363     /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00364     Interpolate(extNpts, extLbda, extK, 
00365                 obsNpts, obsLbda, obsK, gsl_interp_cspline);
00366 
00367     free(extLbda); free(extK);
00368 
00369   } /* EXTINCT */
00370   else {
00371     print_msg("Using default extinction for Mauna-Kea "
00372               "(CFHT Bulletin, number 19, p. 16, 1988)");
00373     for (i=0; i<obsNpts; i++) 
00374       obsK[i] = compute_default_extinction(obsLbda[i]);
00375   }
00376 
00377   /* ----------------------------------------------------------------------- */
00378   /* Computing the corrected spectrum                                        */
00379   /* ----------------------------------------------------------------------- */
00380 
00381   for (i=0; i<obsNpts; i++) 
00382     obsVal[i] /= ( pow(10.,(-0.4*obsK[i]*airmass)) * exptime );
00383 
00384   free(obsK);
00385 
00386   if (DEBUG) {
00387     sprintf(dbgname,"%s_norm",DBG_PREFIX);
00388     print_msg("DEBUG: Saving normalized flux std spectrum as %s",dbgname);
00389     dump_dbgArray((void *)obsVal, DOUBLE, obsNpts, dbgname, &inspec);
00390   }
00391 
00392   if (MATCHRES) {
00393     print_warning("DEBUG DEBUG DEBUG: assuming refRes is a sigma!!!");
00394     
00395     /* Compute spectral resolution degradation factor [px] */
00396     for (i=0; i<refNpts; i++) {
00397       refRes[i] /= inspec.step; /* AA -> px */
00398       if (DEBUG && refRes[i]<specres)
00399         print_warning("Ref. specres at %.2f AA (%.3f AA=%.3f px) smaller than "
00400                       "observed one", refLbda[i], 
00401                       refRes[i]*inspec.step, refRes[i]);
00402       refRes[i] = sqrt(MAX(SQ(refRes[i]) - SQ(specres),0));
00403     }
00404 
00405     /* Interpolate degradation factor [px] */
00406     obsDeg = (double *)malloc(obsNpts*sizeof(double));
00407     Interpolate(refNpts, refLbda, refRes, 
00408                 obsNpts, obsLbda, obsDeg, gsl_interp_linear);
00409 
00410     free(refRes);
00411 
00412     /* Degrade resolution to reference resolution */
00413     /* TODO: skip telluric features */
00414     degradeSpecRes(obsVal, obsNpts, obsDeg);
00415 
00416     if (DEBUG) {
00417       sprintf(dbgname,"%s_degr",DBG_PREFIX);
00418       print_msg("DEBUG: Saving degraded flux std spectrum as %s",dbgname);
00419       dump_dbgArray((void *)obsVal, DOUBLE, obsNpts, dbgname, &inspec);
00420     }
00421   }
00422   
00423 
00424   /* TODO:
00425    * - !matchbin:
00426    *   - interpolate reflFlux on obsLbda -> obsRef
00427    *   - compare obsRef and obsVal -> obsCalib
00428    * - matchbin:
00429    *   - rebin obsVal on refLbda -> refObs 
00430    *   - compare refObs and refFlux -> refCalib
00431    *   - interpolate on obsLbda -> obsCalib 
00432    * - smooth (define and exclude telluric features) 
00433    */
00434 
00435   /* ======================================================================= */
00436   /* Compute the absolute spectrum from interpolation of the absolute flux   */
00437   /* values stored in the refFlux[] array                                    */
00438   /* ======================================================================= */
00439 
00440   obsRef = (double *)malloc(obsNpts*sizeof(double));
00441 
00442   /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00443   /* Interpolating                                                           */
00444   /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
00445   Interpolate(refNpts, refLbda, refFlux, 
00446               obsNpts, obsLbda, obsRef, gsl_interp_cspline);
00447   for (i=0; i<obsNpts; i++) {
00448     if (obsRef[i] <= 0) { /* Paranoid check... */
00449       print_error("ARGH! Reference spectrum <0: Flux(%.2f) = %.2f",
00450                   obsLbda[i],obsRef[i]);
00451       exit_session(ERR_BAD_PARAM);
00452     }
00453   }
00454   free(refLbda); free(refFlux);
00455 
00456   if (DEBUG) {
00457     sprintf(dbgname,"%s_ref",DBG_PREFIX);
00458     print_msg("DEBUG: Saving flux reference spectrum as %s",dbgname);
00459     dump_dbgArray((void *)obsRef, DOUBLE, obsNpts, dbgname, &inspec);
00460   }
00461 
00462   /* ======================================================================= */
00463   /* Compute and save the correction spectrum in ADCU s-1 per flux unit      */
00464   /* ======================================================================= */
00465 
00466   for (i=0; i<obsNpts; i++) 
00467     obsVal[i] /= obsRef[i];
00468 
00469   free(obsRef);
00470 
00471   /* ======================================================================= */
00472   /* Filtering the result. Pixels in discarded intervals are assigned        */
00473   /* negligible weights                                                      */
00474   /* ======================================================================= */
00475 
00476   if (SMOOTH) {
00477 
00478     /* ===================================================================== */
00479     /* Read Table with excluded domains                                      */
00480     /* ===================================================================== */
00481 
00482     if (EXCLUDE) { 
00483       print_msg("Reading exclusion table %s", ExcTabName);
00484       if ((exclNpts = read_2columns(ExcTabName, ColWavsName, ColWaveName, 
00485                                     TBL_FILE, &exclStart, &exclEnd)) < 0) {
00486         exit_session(ERR_BAD_TBL);
00487       }
00488       print_msg("   %d wavelength intervals in exclusion table", exclNpts);
00489 
00490       /* One should test that exclStart[i] < exclEnd[i] */
00491 
00492       if (DEBUG) {
00493         for (i=0; i<exclNpts; i++) {
00494           if ((exclEnd[i] >= obsLMin) && (exclStart[i] <= obsLMax)) { 
00495             print_msg("   %d: %.1f-%.1f AA", i+1, exclStart[i], exclEnd[i]); 
00496           }
00497         }
00498       }
00499       
00500     } /* EXCLUDE */
00501 
00502     print_msg("Polynomial fitting of result spectrum (ncoef=%d)", ncoef);
00503 
00504     /* --------------------------------------------------------------------- */
00505     /* Allocating memory                                                     */
00506     /* --------------------------------------------------------------------- */
00507 
00508     wLbda  = (double *)malloc(obsNpts*sizeof(double));
00509     weight = (double *)malloc(obsNpts*sizeof(double));
00510 
00511     /* --------------------------------------------------------------------- */
00512     /* Filling the wLbda[] array and the weight array                        */
00513     /* --------------------------------------------------------------------- */
00514 
00515     for (obsNExcl=0, i=0; i<obsNpts; i++) {
00516       wLbda[i]  = obsLbda[i];
00517       weight[i] = 1.;
00518       if (EXCLUDE) {
00519         for (j = 0; j<exclNpts; j++) {
00520           if ((wLbda[i] >= exclStart[j]) && (wLbda[i] <= exclEnd[j])) {
00521             weight[i] = 1.e-99; 
00522             obsNExcl++;
00523             break; /* Discarded areas cannot overlap */
00524           }
00525         }
00526       }
00527     }
00528 
00529     /* --------------------------------------------------------------------- */
00530     /* DEBUG mode stuff                                                      */
00531     /* --------------------------------------------------------------------- */
00532 
00533     if (DEBUG) {
00534       sprintf(dbgname,"%s_prefit",DBG_PREFIX);
00535       print_msg("DEBUG: Saving pre-fit spectrum as %s",dbgname);
00536       dump_dbgArray((void *)obsVal, DOUBLE, obsNpts, dbgname, &inspec);
00537 
00538       if (obsNExcl > 0) {
00539         print_msg("%d points (over %d) discarded from polynomial fit",
00540                   obsNExcl, obsNpts);
00541       }
00542       else 
00543         print_msg("No points discarded from polynomial fit");
00544 
00545       sprintf(dbgname,"%s_weight",DBG_PREFIX);
00546       print_msg("DEBUG: Saving weight spectrum as %s",dbgname);
00547       dump_dbgArray((void *)weight, DOUBLE, obsNpts, dbgname, &inspec);
00548     }
00549 
00550     /* --------------------------------------------------------------------- */
00551     /* Let's smooth...                                                       */
00552     /* --------------------------------------------------------------------- */
00553     coef = (double *)malloc(ncoef*sizeof(double));
00554 
00555     fit_xpoly_rej_nag_tab(wLbda, obsVal, weight, 
00556                           obsNpts, ncoef,-3., 3., 20, 
00557                           coef, &niter, &sigerr, &nrej);
00558 
00559     for (i=0; i<obsNpts; i++) 
00560       obsVal[i] = val_minipoly(wLbda[i], coef, ncoef, 
00561                                wLbda[0], wLbda[obsNpts-1]);
00562 
00563     free(wLbda); free(weight); free(coef);
00564 
00565   } /* SMOOTH */
00566 
00567   /* ======================================================================= */
00568   /* Saving the stuff...                                                     */
00569   /* ======================================================================= */
00570 
00571   print_msg("Saving flux correction spectrum in %s", outname);
00572 
00573   create_spec(&outspec, outname, obsNpts, obsLMin, inspec.step, FLOAT, 
00574               inspec.ident, inspec.cunit);
00575   for (i=0; i<obsNpts; i++) 
00576     WR_spec(&outspec, i, obsVal[i]);
00577 
00578   if (strlen(refUnits)) {
00579     print_msg("Adding flux units keyword %s = '%s' to %s", 
00580               FLXUNITS, refUnits, outname);
00581     WR_desc(&outspec, FLXUNITS, CHAR, lg_unit+1, refUnits);
00582   }
00583   else {
00584     print_warning("Please add a %s string-keyword to spectrum %s "
00585                   "to precise flux units.",FLXUNITS,outname);
00586   }
00587 
00588   /* Write type of the throughput spectrum */
00589   write_file_class(&outspec, THR_SPEC); 
00590   close_spec(&outspec);
00591 
00592   /* ======================================================================= */
00593   /* The end my friend, the end...                                           */
00594   /* ======================================================================= */
00595 
00596   free2d(&tabcar1,CHAR);
00597   if (EXCLUDE) free2d(&tabcar3,CHAR); 
00598   if (EXTINCT) free2d(&tabcar4,CHAR); 
00599 
00600   free(obsLbda); free(obsVal);
00601   if (EXCLUDE) {
00602     free(exclStart);
00603     free(exclEnd);
00604   }
00605   exit_session(OK);
00606   return(OK);
00607 }

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