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

calib/source/apply_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 
00022 /* === Doxygen Comment ======================================= */
00028 /* =========================================================== */
00029 
00030 int inter_wspec(SPECTRUM *spec1,SPECTRUM *spec2)                        
00031 {
00032   double l1, l2;
00033 
00034   l1 = MAX(spec1->wstart,spec2->wstart);
00035   l2 = MIN(spec1->wend,spec2->wend);
00036   if (l1 >= l2) return(ERR_BAD_PARAM);
00037   spec1->iwstart = ((l1 - spec1->start) / spec1->step) + 0.5;
00038   spec1->iwend   = ((l2 - spec1->start) / spec1->step) + 0.5;
00039   spec1->wstart  = spec1->start + spec1->iwstart*spec1->step;
00040   spec1->wend    = spec1->start + spec1->iwend*spec1->step;
00041   spec2->iwstart = ((l1 - spec2->start) / spec2->step) + 0.5;
00042   spec2->iwend   = ((l2 - spec2->start) / spec2->step) + 0.5;
00043   spec2->wstart  = spec2->start + spec2->iwstart*spec2->step;
00044   spec2->wend    = spec2->start + spec2->iwend*spec2->step;
00045   return(OK);
00046 }
00047 
00048 /* === Doxygen Comment ======================================= */
00062 /* =========================================================== */
00063 
00064 void Apply_lfff(TIGERfile *in, TIGERfile *lfff, TIGERfile *out, int COL, 
00065                 int TSKY, char *colname1, char *colname2)
00066 {
00067   SPECTRUM sin, slfff, sout;
00068   SPECTRUM nin,nout;
00069   TABLE tab;
00070   int *nolens,  NOISE, k1, k2, flag, *row;
00071   int nbsel;
00072   int i, j, k, l, npts, nrej, ntrunc, nzero, col_no, col_flag,col_lb1, col_lb2;
00073   float var, res, valff, lb1, lb2, *A, *D;
00074   char table_name[lg_name+1];
00075 
00076   nrej = 0;
00077   ntrunc = 0;
00078 
00079   nolens = (int *)malloc(in->nbspec*sizeof(int));
00080   A    = (float *)malloc(in->nbspec*sizeof(float));
00081   D    = (float *)malloc(in->nbspec*sizeof(float));
00082   row  =   (int *)malloc(in->nbspec*sizeof(int));
00083   get_lenses_coord_select(in,"A","D",nolens,A,D,row,&nbsel);
00084 
00085   if (COL) {
00086     strcpy(table_name, lfff->table_name);
00087     print_msg("Opening Table %s Columns %s, %s", table_name,colname1, colname2);
00088     open_table(&tab, table_name, "I");
00089     col_no   = get_col_ref(&tab, "No");
00090     col_lb1  = get_col_ref(&tab, colname1);
00091     col_lb2  = get_col_ref(&tab, colname2);
00092     col_flag = get_col_ref(&tab, "FLAG");
00093   }
00094 
00095   set_control_level(NONE);
00096   handle_select_flag(&tab, 'W', NULL);
00097   NOISE = exist_noise(in);
00098 
00099   /* Loop over input lenses */
00100   for (i=0; i<nbsel; i++) {
00101     print_progress("Flat-Fielding",100.*(i+1)/in->nbspec,1);
00102 
00103     if (NOISE) get_tiger_spec(in, &sin, &nin, nolens[i]);
00104     else       get_tiger_spec(in, &sin, NULL, nolens[i]);
00105 
00106     if (get_tiger_spec(lfff, &slfff, NULL, nolens[i]) != 0) {
00107       print_warning("Spectra %d discarded (no flat field spectrum)",nolens[i]);
00108       free_spec_mem(&sin);
00109       if (NOISE) free_spec_mem(&nin);
00110       nrej++;
00111       continue;
00112     }
00113 
00114     if (COL) {
00115       l = search_in_col(&tab, col_no, &(nolens[i])); 
00116       if (RD_tbl(&tab, l, col_flag, &flag) == ERR_NODATA) {
00117         print_warning("Spectra %d discarded (no common pixel)",nolens[i]);
00118         free_spec_mem(&sin);
00119         if (NOISE) free_spec_mem(&nin);
00120         free_spec_mem(&slfff);
00121         nrej++;
00122         continue;
00123       }
00124       /* special treatment if object or (TSKY and sky) */
00125       if ((flag == 1) || (TSKY && flag == 2)) {
00126         if ((RD_tbl(&tab, l, col_lb1, &lb1) == ERR_NODATA) ||
00127             (RD_tbl(&tab, l, col_lb2, &lb2) == ERR_NODATA)) 
00128           npts = 0;
00129         else {
00130           k1  = pixel_spec(&slfff, lb1); 
00131           lb1 = coord_spec(&slfff, k1);
00132           slfff.iwstart = k1; 
00133           slfff.wstart  = lb1;
00134           k2  = pixel_spec(&slfff, lb2); 
00135           lb2 = coord_spec(&slfff, k2);
00136           slfff.iwend = k2; 
00137           slfff.wend  = lb2;
00138           
00139           inter_wspec(&sin, &slfff);
00140           npts = sin.iwend - sin.iwstart + 1; 
00141         } 
00142       } else { /* usual treatment otherwise */
00143         inter_spec(&sin, &slfff);
00144         npts = sin.iwend - sin.iwstart + 1; 
00145       }
00146     } else { /* usual treatment */
00147       inter_spec(&sin, &slfff);
00148       npts = sin.iwend - sin.iwstart + 1; 
00149     }
00150     
00151     if (npts <= 0) {
00152       print_warning("Spectra %d discarded (no common pixel)", nolens[i]);
00153       free_spec_mem(&sin);
00154       if (NOISE) free_spec_mem(&nin);
00155       free_spec_mem(&slfff);
00156       nrej++;
00157       continue;
00158     }
00159     if (npts < sin.npts) ntrunc++; /* truncated spectra */
00160   
00161     init_new_tiger_spec(out,&sout,npts,sin.wstart);
00162     if (NOISE) init_new_tiger_spec(out,&nout,npts,sin.wstart);
00163 
00164     /* Flat-fielding */
00165     nzero = 0;
00166     for(j=0, k=sin.iwstart, l=slfff.iwstart; j<npts; j++, k++, l++) {
00167       valff = RD_spec(&slfff,l);
00168       if (valff == 0) { /* cannot flat-field */
00169         nzero++;
00170         WR_spec(&sout,j,0);
00171         if (NOISE) WR_spec(&nout,j,-1);
00172       } else {
00173         res = RD_spec(&sin,k)/valff;
00174         WR_spec(&sout,j,res);
00175         if (NOISE) {
00176           var = RD_spec(&nin,k)/(valff*valff);
00177           WR_spec(&nout,j,var);
00178         }
00179       }
00180     }
00181     if (nzero > 0) 
00182       print_warning("%d flagged bad points in spectrum %d",nzero, nolens[i]);
00183 
00184     /* Store the result */
00185     if (NOISE) put_tiger_spec(out, &sout, &nout, nolens[i]);
00186     else       put_tiger_spec(out, &sout, NULL, nolens[i]);
00187 
00188     free_spec_mem(&sin);
00189     if (NOISE) free_spec_mem(&nin);
00190     free_spec_mem(&slfff);        
00191   }
00192   print_msg("%d spectra truncated and %d spectra rejected", ntrunc, nrej);
00193 
00194   free(nolens);
00195 }
00196 
00197 /* === Doxygen Comment ======================================= */
00210 /* =========================================================== */
00211 
00212 int main(int argc,char **argv)
00213 {
00214   TIGERfile c_in, c_lfff, c_out;
00215   char **label, **param, **colname, *in, *out, *lfff;
00216   int npix = -1, class, COL, TSKY;
00217   double start=0,end=0;
00218 
00219   set_arglist("-in none -lfff none -out none -col null -tsky");
00220   set_purpose("Low-frequency spectro-spatial flat-field (application)");
00221   
00222   init_snifs("$Name:  $");
00223   init_session(argv, argc, &label, &param);
00224 
00225   if (DEBUG) {
00226     print_msg("$Id: apply_lfff.c,v 1.5 2004/10/06 16:56:35 ycopin Exp $");
00227     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00228   }
00229 
00230   /* Inputs */
00231 
00232   in   = param[0];
00233   lfff = param[1];
00234   out  = param[2];
00235 
00236   print_msg("Datacubes Input: %s Output: %s", in, out);
00237   print_msg("References Flatfield Datacube: %s", lfff);
00238 
00239   colname    = (char **)malloc(2*sizeof(char *));
00240   colname[0] = (char *)malloc(16*sizeof(char));
00241   colname[1] = (char *)malloc(16*sizeof(char));     
00242 
00243   if ((COL = is_set(param[3]))) {
00244     if (decode_argval_char(param[3], colname) != 2) {
00245       print_error("Error in input parameter -col %s", param[3]);
00246       exit_session(ERR_BAD_PARAM);
00247     }
00248     if ((TSKY = is_true(param[4]))) 
00249       print_msg("Object AND Sky spectra truncated according to %s,%s",
00250                 colname[0], colname[1]);
00251     else
00252       print_msg("Only Object spectra truncated according to %s,%s",
00253                 colname[0], colname[1]);
00254   }
00255   
00256   if (open_tiger_frame(&c_in,in,"I") < 0){
00257     print_error("Unable to access the input datacube");
00258     exit_session(ERR_OPEN);
00259   }
00260 
00261   /* check input file class */
00262   class = read_file_class(&c_in);
00263   switch (class) {
00264   case SUPER_CLASS:
00265   case WAV_OBJ_CUBE: 
00266   case WAV_SKY_CUBE:
00267   case WAV_DOM_CUBE:
00268   case WAV_CON_CUBE: break;
00269   default:
00270     print_error("Wrong file class (should be lambda-calibrated datacube)");
00271     exit_session(ERR_BAD_PARAM);
00272   }
00273 
00274   if (open_tiger_frame(&c_lfff, lfff, "I") < 0){
00275     close_tiger_frame(&c_in);
00276     print_error("Unable to access the Flat datacube");
00277     exit_session(ERR_OPEN);
00278   }
00279 
00280   /* check FF file class */
00281   switch (read_file_class(&c_lfff)) {
00282   case SUPER_CLASS:
00283   case FLAT_CUBE: break;
00284   default:
00285     print_error("Wrong file class (should be flat-field datacube)");
00286     exit_session(ERR_BAD_PARAM);
00287   }
00288 
00289   if (ABS(c_in.step - c_lfff.step) > 0.0000001){ /* strange test... */
00290     close_tiger_frame(&c_in);
00291     close_tiger_frame(&c_lfff);
00292     print_error("Input datacube and Flat datacube have different steps");
00293     exit_session(ERR_BAD_PARAM);
00294   }
00295 
00296   if (has_common_bounds(&c_in)) get_common_param(&c_in,&npix,&start,&end);
00297   else {
00298     npix = -1;
00299     start = -1;
00300   }
00301 
00302   if (create_tiger_frame(&c_out,out,npix,start,c_in.step,c_in.data_type,
00303                          c_in.table_name,c_in.ident,c_in.cunit) < 0){
00304     close_tiger_frame(&c_in);
00305     close_tiger_frame(&c_lfff);
00306     print_error("Cannot create output datacube");
00307     exit_session(ERR_OPEN);
00308   }
00309 
00310   Apply_lfff(&c_in, &c_lfff, &c_out, COL, TSKY, colname[0], colname[1]);
00311 
00312   /* write output datacube descriptors */
00313   CP_non_std_desc(&c_in, &c_out);  
00314   switch (class) {
00315   case WAV_OBJ_CUBE: write_file_class(&c_out, FLA_OBJ_CUBE); break;
00316   case WAV_SKY_CUBE: write_file_class(&c_out, FLA_SKY_CUBE); break;
00317   case WAV_DOM_CUBE: write_file_class(&c_out, FLA_DOM_CUBE); break;
00318   case WAV_CON_CUBE: write_file_class(&c_out, FLA_CON_CUBE); break;
00319   case SUPER_CLASS:  write_file_class(&c_out, SUPER_CLASS);  break;
00320   }
00321 
00322   close_tiger_frame(&c_in);
00323   close_tiger_frame(&c_lfff);
00324   close_tiger_frame(&c_out);
00325 
00326   exit_session(OK);
00327   return(OK);
00328 }

Generated on Tue Nov 23 18:04:18 2004 for Snifs by doxygen 1.3.3