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

calib/source/extract_flux.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00014 /* =========================================================== */
00015 
00016 #include <IFU_io.h>
00017 #include <IFU_math.h>
00018 #include <snifs.h>
00019 #include <calib.h>
00020 #include <gsl/gsl_sf_erf.h>
00021 
00022 /* === Doxygen Comment ======================================= */
00049 /* =========================================================== */
00050 
00051 int main(int argc,char **argv) 
00052 {
00053   TIGERfile incube;
00054   SPECTRUM inspec, outspec;
00055 
00056   char **label, **param, **tabcar;
00057   char *inname, *outname, *colXname, *colYname;
00058 
00059   int *no, *il, i, l, k, nspec, nsum, nlbda, PSF;
00060 
00061   float *X, *Y;
00062   float xc, yc, dx, dy, radius, r2, r2max, sig1, sig2, ratio;
00063 
00064   double *sum, bounds[4];
00065   double LbdaStart, LbdaEnd, corr;
00066 
00067   /* ====================================================================== */
00068   /* Starting the session                                                   */
00069   /* ====================================================================== */
00070 
00071   set_purpose("Flux calibration (crude standard spectrum extraction)");
00072   set_arglist("-in none -out none -col A,D -center none -radius none -psf null");
00073   
00074   init_snifs("$Name:  $","$Revision: 1.3 $");
00075   init_session(argv, argc, &label, &param);
00076 
00077   if (DEBUG) {
00078     print_msg("$Id: extract_flux.c,v 1.3 2005/10/05 22:44:30 ycopin Exp $");
00079     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00080   }
00081 
00082   /* ====================================================================== */
00083   /* Reading the input parameters                                           */
00084   /* ====================================================================== */
00085 
00086   inname  = param[0];                        /* Input datacube */
00087   outname = param[1];                        /* Output flux correction */
00088 
00089   /* Input position columns */
00090   alloc2d(&tabcar,2,lg_label+1,CHAR);
00091   if (decode_argval_char(param[2], tabcar) != 2) {
00092     print_error("Wrong argument -col col(X),col(Y) (got %s)",param[2]);
00093     exit_session(ERR_BAD_PARAM);
00094   }
00095   colXname = tabcar[0];
00096   colYname = tabcar[1];
00097 
00098   /* Summation region (center and radius) */
00099   sscanf(param[3], "%f,%f", &xc, &yc);       /* -center x,y */
00100   get_argval(4, "%f", &radius);              /* -radius r */
00101 
00102   /* PSF estimates */
00103   if ((PSF = is_set(param[5]))) {
00104     if (sscanf(param[5],"%f,%f,%f", &sig1, &sig2, &ratio) != 3) {
00105       print_error("Wrong argument -psf S1,S2,I2/I1 (got %s)",param[5]);
00106       exit_session(ERR_BAD_PARAM);
00107     }
00108     if (sig1 <= sig2 || sig2 < 0 || ratio < 0) {
00109       print_error("Wrong PSF parameters (S1 > S2 >= 0, I2/I1 >= 0)");
00110       exit_session(ERR_BAD_PARAM);
00111     }
00112   }
00113 
00114   print_msg("o Input flux standard datacube: %s", inname);
00115   print_msg("o Output flux standard spectrum: %s", outname);
00116   print_msg("o Flux standard position: %s,%s = %.2f,%.2f",
00117             colXname, colYname, xc, yc);
00118   if (PSF)
00119     print_msg("o Input PSF estimate: S1=%.2f, S2=%.2f, I2/I1=%.2f",
00120               sig1, sig2, ratio);
00121 
00122   /* ======================================================================= */
00123   /* Open the input datacube and check its class                             */
00124   /* ======================================================================= */
00125 
00126   print_msg("Opening input datacube %s", inname);
00127   if (open_tiger_frame(&incube, inname, "I") < 0) {
00128     print_error("Cannot open datacube %s", inname);
00129     exit_session(ERR_OPEN);
00130   }
00131 
00132   switch (read_file_class(&incube)) {
00133   case SUPER_CLASS:
00134   case FLA_OBJ_CUBE:
00135   case COS_OBJ_CUBE:
00136   case SKY_OBJ_CUBE: break;
00137   default:
00138     print_error("Wrong fclass for %s (should be object datacube)",inname);
00139     exit_session(ERR_BAD_TYPE);
00140   }
00141 
00142   nspec = incube.nbspec;
00143 
00144   /* ======================================================================= */
00145   /* Compute the wavelength range                                            */
00146   /* 1st case  : if the wavelength range is the same for all spectra, uses   */
00147   /*             directly the output of the get_common_param() routine       */
00148   /* 2nd case  : uses the intersection of the all the wavelength ranges of   */
00149   /*             spectra within the summation region                         */
00150   /* ======================================================================= */
00151 
00152   no =   (int *)malloc(nspec*sizeof(int));
00153   il =   (int *)malloc(nspec*sizeof(int));
00154   X  = (float *)malloc(nspec*sizeof(float));
00155   Y  = (float *)malloc(nspec*sizeof(float));
00156 
00157   get_lenses_coord(&incube, colXname, colYname, no, X, Y, il);
00158   r2max = radius*radius;
00159 
00160   if (has_common_bounds(&incube)) {
00161     /* --------------------------------------------------------------------- */
00162     /* Case 1 : the spectra have common boundaries                           */
00163     /* --------------------------------------------------------------------- */
00164     get_common_param(&incube, &nlbda, &LbdaStart, &LbdaEnd);
00165   }
00166   else {
00167     /* --------------------------------------------------------------------- */
00168     /* Case 2 : the spectra cab have different lengths, takes the            */
00169     /*          intersection of all the wavelength ranges within aperture    */
00170     /* --------------------------------------------------------------------- */
00171     if (RD_desc(&incube,LRANGE4,DOUBLE,4,bounds) < 0) {
00172       print_error("Cannot read descriptor %s in datacube %s",LRANGE4,inname);
00173       exit_session(ERR_NODESC);
00174     }
00175 
00176     LbdaStart = bounds[0];
00177     LbdaEnd = bounds[1];
00178     for (l=0; l<nspec; l++) {
00179       dx = X[l] - xc;
00180       dy = Y[l] - yc;
00181       if (dx*dx + dy*dy > r2max) continue;
00182       get_tiger_spec(&incube, &inspec, NULL, no[l]);
00183       LbdaStart = MAX(LbdaStart,inspec.start);
00184       LbdaEnd   = MIN(LbdaEnd,  inspec.end);
00185       free_spec_mem(&inspec);
00186     }
00187     nlbda = (int)((LbdaEnd - LbdaStart)/incube.step + 1.5);
00188   }
00189 
00190   print_msg("Spectral range: %.2f-%.2f AA", LbdaStart, LbdaEnd);
00191 
00192   /* ======================================================================= */
00193   /* Summing the spectra within a given radius, the summed spectrum is       */
00194   /* first stored in the sum[] double array.                              */
00195   /* ======================================================================= */
00196 
00197   sum = (double *)calloc(nlbda, sizeof(double));
00198 
00199   reset_print_progress();
00200   for (nsum=0, l=0; l<nspec; l++) {
00201     print_progress("Summing spectra within given radius ", 
00202                    (l+1.)/nspec*100,5.);
00203     dx = X[l] - xc;
00204     dy = Y[l] - yc;
00205     r2 = dx*dx + dy*dy;
00206     if (r2 > r2max) continue;
00207     else if (DEBUG) {
00208       print_msg("%3d: Lens #%3d (%.2f x %.2f) at %.2f",l+1,no[l],X[l],Y[l],r2);
00209     }
00210 
00211     get_tiger_spec(&incube, &inspec, NULL, no[l]);
00212     if (set_subspec(&inspec, LbdaStart, LbdaEnd) != 0) {
00213       print_error("Lens #%d, set_subspec() != 0",no[l]);
00214       exit_session(UNKNOWN);
00215     }
00216 
00217     for (i=0, k=inspec.iwstart; k<=inspec.iwend; i++, k++) 
00218       sum[i] += RD_spec(&inspec, k);
00219     free_spec_mem(&inspec);
00220     nsum++;
00221   }
00222 
00223   print_msg("-> %d spectra summed within radius %.2f", nsum, radius);
00224 
00225   if (!nsum) {
00226     print_error("No valid spectrum inside aperture");
00227     exit_session(ERR_BAD_PARAM);
00228   }
00229 
00230   /* ======================================================================= */
00231   /* Using the PSF estimate to compute a correction factor on the summed     */
00232   /* flux spectrum (how much flux did we miss by stopping the summation at   */
00233   /* a given radius)                                                         */
00234   /* ======================================================================= */
00235 
00236   if (PSF) {
00237     /* --------------------------------------------------------------------- */
00238     /* Computing the correction                                              */
00239     /* --------------------------------------------------------------------- */
00240     if (ratio != 0) {
00241       corr = (sig2 + ratio*sig1) /
00242         (       sig2*gsl_sf_erf( radius/M_SQRT2/sig1 ) + 
00243           ratio*sig1*gsl_sf_erf( radius/M_SQRT2/sig2) );
00244     }
00245     else {
00246       corr = 1/gsl_sf_erf( radius/M_SQRT2/sig1 );
00247     }
00248     print_msg("Computed PSF correction factor: 1 + %g", corr-1);
00249 
00250     /* --------------------------------------------------------------------- */
00251     /* Applying the correction                                               */
00252     /* --------------------------------------------------------------------- */
00253     for (i=0; i<nlbda; i++) sum[i] *= corr;
00254   } 
00255 
00256   /* ======================================================================= */
00257   /* Create the output spectrum                                              */
00258   /* ======================================================================= */
00259 
00260   create_spec(&outspec, outname, nlbda, LbdaStart, incube.step, FLOAT, 
00261               incube.ident, incube.cunit);
00262   for (i=0; i<nlbda; i++) WR_spec(&outspec, i, sum[i]);
00263   CP_non_std_desc(&incube,&outspec);
00264   /* Should we define a special FCLASS? */
00265   /* Should we add keywords about center, coord., PSF? */
00266   close_spec(&outspec);
00267 
00268   /* ======================================================================= */
00269   /* The end my friend, the end...                                           */
00270   /* ======================================================================= */
00271 
00272   close_tiger_frame(&incube);
00273 
00274   free(X); free(Y); free(no); free(il); free(sum);
00275   free2d(&tabcar,CHAR);
00276 
00277   exit_session(OK);
00278   return(OK);
00279 }

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