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

calib/source/apply_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 
00023 /* === Doxygen Comment ======================================= */
00033 /* =========================================================== */
00034 
00035 void Apply_corrflux(TIGERfile *incube, SPECTRUM *corrflux, TIGERfile *outcube, 
00036                     float exptime, float airmass) {
00037   SPECTRUM inspec,outspec;
00038   SPECTRUM innoise,outnoise;
00039 
00040   int *nolens;
00041   int i,j,npix,k,l,nrej,ntrunc, NOISE;
00042 
00043   float Fabs, frej, ftrunc;
00044   double lbda, val, var, kext, coef;
00045 
00046   nolens = (int *)malloc(incube->nbspec*sizeof(int));
00047   get_lenses_no(incube,nolens);
00048   nrej = ntrunc = 0;
00049 
00050   NOISE = exist_noise(incube);
00051 
00052   reset_print_progress();
00053   for (i=0; i<incube->nbspec; i++) {         /* Loop on input spectra */
00054 
00055     print_progress("Flux calibration ", (i+1.)/incube->nbspec*100.,5);
00056     
00057     get_tiger_spec(incube,&inspec,(NOISE? &innoise:NULL),nolens[i]);
00058 
00059     /* Find intersection between spectrum and flux-correction */
00060     inter_spec(&inspec, corrflux);
00061     npix = inspec.iwend - inspec.iwstart + 1;
00062     if (npix <= 0) {
00063       print_warning("Spectrum %d rejected (no common bounds)", nolens[i]);
00064       free_spec_mem(&inspec);
00065       if (NOISE) free_spec_mem(&innoise);
00066       nrej++;
00067       continue;
00068     }
00069     if (npix < inspec.npts) ntrunc++;
00070 
00071     init_new_tiger_spec(outcube,&outspec,npix,inspec.wstart);
00072     if (NOISE) 
00073       init_new_tiger_spec(outcube,&outnoise,npix,inspec.wstart);
00074 
00075     /* Flux correction */
00076     for (k=inspec.iwstart, l=corrflux->iwstart, j=0; j<npix; j++, k++, l++) {
00077       lbda = coord_spec(&outspec,j);
00078       kext = compute_default_extinction(lbda);
00079       Fabs = pow(10.,0.4*kext*airmass);
00080       if (RD_spec(corrflux,l) == 0) {
00081         val = 0.;
00082         if (NOISE) var = -1.;
00083       } 
00084       else {
00085         coef = Fabs/(RD_spec(corrflux,l) * exptime);
00086         val  = (double)RD_spec(&inspec,k) * coef;
00087         if (NOISE) {
00088           if (RD_spec(&innoise,k) < 0) var = -1;
00089           else var = (double)RD_spec(&innoise,k) * SQ(coef);
00090         }
00091       }
00092       WR_spec(&outspec,j,val);
00093       if (NOISE) WR_spec(&outnoise,j,var);
00094     }
00095     free_spec_mem(&inspec);
00096     if (NOISE) free_spec_mem(&innoise);
00097     put_tiger_spec(outcube,&outspec,(NOISE? &outnoise:NULL),nolens[i]);
00098 
00099   } /* Loop on input spectra */
00100 
00101   free(nolens);
00102 
00103   /* Some statistics */
00104   if ((ntrunc > 0) || (nrej > 0)) {
00105     frej   = 100.*nrej/incube->nbspec;     /* Fraction of rejected spectra */
00106     ftrunc = 100.*ntrunc/incube->nbspec;   /* Fraction of truncated spectra */
00107     print_warning("%d (%5.2f%%) spectra rejected, "
00108                   "%d (%5.2f%%) spectra truncated", 
00109                   nrej, frej, ntrunc, ftrunc);
00110   }
00111 }
00112 
00113 /* === Doxygen Comment ======================================= */
00129 /* =========================================================== */
00130 
00131 int main(int argc,char **argv)
00132 {
00133   TIGERfile incube,outcube;
00134   SPECTRUM corrflux;
00135 
00136   char **label,**param;
00137   char *inname,*outname,*fluxname;
00138   char flxunits[lg_unit+1];
00139 
00140   int npix = -1;
00141   float exptime, airmass;
00142   double start=0,end=0;
00143   
00144   set_purpose("Flux calibration (application)");
00145   set_arglist("-in none -cflux none -out none");
00146 
00147   init_snifs("$Name:  $","$Revision: 1.12 $");
00148   init_session(argv, argc, &label, &param);
00149 
00150   if (DEBUG) {
00151     print_msg("$Id: apply_flux.c,v 1.12 2005/10/21 12:17:45 ycopin Exp $");
00152     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00153   }
00154 
00155   inname   = param[0];
00156   fluxname = param[1];
00157   outname  = param[2];
00158   
00159   /* Open input datacube */
00160   print_msg("Opening input datacube %s", inname);
00161   if (open_tiger_frame(&incube,inname,"I")<0) {
00162     print_error("Cannot open datacube %s",inname);
00163     exit_session(ERR_OPEN);
00164   }
00165 
00166   /* Test input fclass */
00167   switch (read_file_class(&incube)) {
00168   case SUPER_CLASS:
00169   case FLA_OBJ_CUBE:
00170   case COS_OBJ_CUBE:
00171   case SKY_OBJ_CUBE: break;
00172   default:
00173     print_error("Wrong fclass for %s (should be object datacube)",inname);
00174     exit_session(ERR_BAD_TYPE);
00175   }
00176 
00177   /* Read airmass descriptor */
00178   if ((airmass = read_airmass(&incube)) < 1) {
00179     print_error("Cannot read descriptor %s in datacube %s", AIRMASS, inname);
00180     exit_session(ERR_BAD_DESC);
00181   }
00182   /* Read exposure time descriptor */
00183   if ((exptime = read_exptime(&incube)) < 0) {
00184     print_error("Cannot read descriptor %s in datacube %s", EFFTIME, inname);
00185     exit_session(ERR_BAD_DESC);
00186   }
00187   print_msg("   Integration time: %.1fs, Airmass: %.2f", exptime, airmass);
00188 
00189   /* Open flux-correction spectrum */
00190   print_msg("Opening flux correction spectrum %s", fluxname);
00191   if (open_spec(&corrflux, fluxname,"I") < 0) {
00192     print_error("Cannot open spectrum %s",fluxname);
00193     exit_session(ERR_OPEN);
00194   }
00195 
00196   disable_user_warnings();
00197   if (RD_desc(&corrflux, FLXUNITS, CHAR, lg_unit+1, flxunits) < 0) {
00198     print_warning("Cannot read flux units (from %s) in %s", FLXUNITS, fluxname);
00199     strcpy(flxunits,"");
00200   } 
00201   else {
00202     print_msg("Flux units: %s",flxunits);
00203   }
00204   restore_user_warnings();
00205 
00206   /* Test lambda-step coherence between input datacube and flux correction */
00207   if (ABS(incube.step - corrflux.step) > 1e-7) {
00208     print_error("Input DataCube and Input Spectrum do not have the same step");
00209     close_tiger_frame(&incube);
00210     close_spec(&corrflux);
00211     exit_session(ERR_BAD_PARAM);
00212   }
00213 
00214   /* Create output datacube */
00215   print_msg("Creating output datacube %s", outname);
00216   if (has_common_bounds(&incube)) {
00217     get_common_param(&incube,&npix,&start,&end);
00218     if (create_tiger_frame(&outcube,outname,npix,start,incube.step,incube.data_type,
00219                            incube.table_name,incube.ident,incube.cunit)<0) {
00220       print_error("Cannot create datacube %s",outname);
00221       exit_session(ERR_CREAT);
00222     }
00223   } 
00224   else {
00225     if (create_tiger_frame(&outcube,outname,-1,-1,incube.step,incube.data_type,
00226                           incube.table_name,incube.ident,incube.cunit)<0) {
00227       print_error("Cannot create datacube %s",outname);
00228       exit_session(ERR_CREAT);
00229     }
00230   }
00231 
00232   /* Apply flux correction */
00233   Apply_corrflux(&incube, &corrflux, &outcube, exptime, airmass);
00234 
00235   /* Output descriptors */
00236   CP_non_std_desc(&incube, &outcube);
00237   WR_desc(&outcube, FLXUNITS, CHAR, lg_unit+1, flxunits);
00238   WR_desc(&outcube, FLXSPEC,  CHAR, lg_name+1, corrflux.name);
00239   switch (read_file_class(&incube)) {
00240   case SUPER_CLASS:
00241   case FLA_OBJ_CUBE:
00242   case COS_OBJ_CUBE:
00243   case SKY_OBJ_CUBE: write_file_class(&outcube, FLX_OBJ_CUBE); break;
00244   }
00245 
00246   /* Closing */
00247   close_spec(&corrflux);
00248   close_tiger_frame(&incube);
00249   close_tiger_frame(&outcube);
00250 
00251   exit_session(OK);
00252   return(OK);
00253 }

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