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

calib/source/wrebin.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00018 /* =========================================================== */
00019 
00020 #include <IFU_io.h>
00021 #include <IFU_math.h>
00022 #include <snifs.h>
00023 
00024 /* === Doxygen Comment ======================================= */
00036 /* =========================================================== */
00037 
00038 double newLambda(double lambda, 
00039                  double *coef, long ncoef, double xmin, double xmax, 
00040                  int *status)
00041 {
00042   double new;
00043 
00044   new = lambda + val_poly_nag(lambda, coef, ncoef, xmin, xmax, status);
00045 
00046   return(new);
00047 }
00048 
00049 /* === Doxygen Comment ======================================= */
00063 /* =========================================================== */
00064 
00065 int main(int argc,char **argv)
00066 {
00067   TABLE table_cal, table_cal2;
00068   TIGERfile cubein, cubeout;
00069   SPECTRUM specin, specout, noisein, noiseout;
00070 
00071   char **arglabel,**argval;
00072   char *inname, *outname, ident[lg_ident+1];
00073   char *tblname, *tblname2;
00074 
00075   long ncoef, ncoef2;
00076   int *col_c,  col_no,  col_xmin,  col_xmax,  col_deg,  col_flag,  col_nrej;
00077   int *col2_c, col2_no, col2_xmin, col2_xmax, col2_deg, col2_flag, col2_nrej;
00078   int ligne,  maxdegre_cal,  degre,  status;
00079   int ligne2, maxdegre_cal2, degre2, status2;
00080   int *nolens, i, j, k, l;
00081   int LOG, CAL2, NOISE, flag, flagmax, fclass, fclass2;
00082   int status_end, status_start, ifail, nb_notrebin=0;
00083   int nbpix, kmin, kmax, npts_rebin, k_start, k_end;
00084 
00085   float rdesc, rdesc2;
00086   double *value_in, *value_out, *value_noisein, *value_noiseout;
00087   double *coef, *coef2;
00088   double xmin_pol,  xmax_pol,  eval_tmp,  xinf_rebin,  xsup_rebin;
00089   double xmin_pol2, xmax_pol2, eval_tmp2, xinf_rebin2, xsup_rebin2;
00090   double range[4], step_in, step_rebin_log, step_rebin_lin;
00091   double lambdaref_start, lambdaref_end, lambda_inf, lambda_sup, lambda;
00092   double xmin_in, xmax_in, start_rebin, end_rebin, xinf_pix, xsup_pix;
00093   double min_end_rebin, max_start_rebin;
00094   double frac, frac_start, frac_end;
00095 
00096   set_purpose("Wavelength calibration (application)");
00097   set_arglist("-in none -cal|caltable none -out none -log -quality 2 "
00098               "-cal2|caltable2 null");
00099 
00100   init_snifs("$Name:  $","$Revision: 1.14 $");
00101   init_session(argv,argc,&arglabel,&argval);
00102 
00103   if (DEBUG) {
00104     print_msg("$Id: wrebin.c,v 1.14 2005/10/23 23:32:16 ycopin Exp $");
00105     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00106   }
00107 
00108   inname  = argval[0];
00109   tblname = argval[1];
00110   outname = argval[2];
00111   LOG = is_true(argval[3]);
00112   get_argval(4,"%d",&flagmax);
00113   if ((CAL2 = is_set(argval[5]))) tblname2 = argval[5];
00114 
00115   /*=========================================================================*/
00116   /* Ouverture du cube TIGER a rebinner */
00117   /*=========================================================================*/
00118   if (open_tiger_frame(&cubein,inname,"I")<0) {
00119     print_error("Cannot open datacube %s", inname);
00120     exit_session(ERR_OPEN);
00121   }
00122   switch ((fclass = read_file_class(&cubein))) {
00123   case SUPER_CLASS:
00124   case RAW_OBJ_CUBE:
00125   case RAW_SKY_CUBE:
00126   case RAW_DOM_CUBE:
00127   case RAW_CON_CUBE:
00128   case RAW_CAL_CUBE: break;
00129   default:
00130     print_error("Wrong fclass (%d) for %s (should be a `Raw Datacube')",
00131                 fclass,cubein.name);
00132     print_warning("WE SHOULD STOP HERE but we let it go for the moment");
00133     //      exit_session(ERR_BAD_TYPE);
00134   }
00135   NOISE = exist_noise(&cubein); /* check if noise exist */
00136 
00137   /*=========================================================================*/
00138   /* On cherche les numeros des lentilles */
00139   /*=========================================================================*/
00140   nolens = (int *)malloc(cubein.nbspec*sizeof(int));
00141   get_lenses_no(&cubein,nolens);
00142 
00143   /*=========================================================================*/
00144   /* Ouverture de la table associee au cube de calibration */
00145   /*=========================================================================*/
00146   print_msg("Opening Calibration table '%s'",tblname);
00147   if (open_table(&table_cal, tblname, "IO") < 0) {
00148     close_tiger_frame(&cubein);
00149     print_error("Cannot open table %s",tblname);
00150     exit_session(ERR_OPEN);
00151   }
00152   if ((fclass2 = read_file_class(&table_cal)) != TBL_FILE_CAL) {
00153     print_error("Wrong fclass (%d) for %s (should be a "
00154                 "`Wavelength Calibrated Table')",fclass2,table_cal.name);
00155     exit_session(ERR_BAD_TYPE);
00156   }
00157 
00158   /*=========================================================================*/
00159   /* Ouverture des colonnes de l'ajustement du polynome de calibration */
00160   /*=========================================================================*/
00161   col_no   = get_col_ref(&table_cal, LAB_COL_NO);
00162   col_xmin = get_col_ref(&table_cal, "XMIN");
00163   col_xmax = get_col_ref(&table_cal, "XMAX");
00164   col_deg  = get_col_ref(&table_cal, "DEGREE");
00165   col_flag = get_col_ref(&table_cal, "FLAG");
00166   col_nrej = get_col_ref(&table_cal, "NREJ");
00167 
00168   /*=========================================================================*/
00169   /* Lecture du degre maximum et ouverture des colonnes des coefs */
00170   /*=========================================================================*/
00171   if (RD_desc(&table_cal, WLMAXDEGRE, INT, 1, &maxdegre_cal) < 0) {
00172     close_tiger_frame(&cubein);
00173     print_error("Cannot read descriptor %s in table %s",
00174                 WLMAXDEGRE,table_cal.name);
00175     exit_session(ERR_NODESC);
00176   }
00177   col_c    =    (int *)malloc(maxdegre_cal*sizeof(int));
00178   coef     = (double *)malloc(maxdegre_cal*sizeof(double));
00179   for (i=0; i<maxdegre_cal; i++) {
00180     sprintf(ident,"C%d", i); /* colonnes des coefficients du polynome */
00181     col_c[i] = get_col_ref(&table_cal, ident);
00182   }
00183 
00184   /* ====================================================================== */
00185   /* Is there any second calibration table ? */
00186   /* ====================================================================== */
00187 
00188   if (CAL2) {
00189     print_msg("Opening 2nd calibration table: '%s'",tblname2);
00190     /* Ouverture de la table associee au 2nd cube de calibration */
00191     if (open_table(&table_cal2, tblname2, "IO") < 0) {
00192       close_tiger_frame(&cubein);
00193       print_error("Cannot open table %s",tblname2);
00194       exit_session(ERR_OPEN);
00195     }
00196     if ((fclass2 = read_file_class(&table_cal2)) != TBL_FILE_CAL) {
00197       print_error("Wrong fclass (%d) for %s (should be a "
00198                   "`Wavelength Calibrated Table')",fclass2,table_cal2.name);
00199       exit_session(ERR_BAD_TYPE);
00200     }
00201 
00202     /* Ouverture des colonnes de l'ajustement du polynome de calibration */
00203     col2_no   = get_col_ref(&table_cal2, LAB_COL_NO);
00204     col2_xmin = get_col_ref(&table_cal2, "XMIN");
00205     col2_xmax = get_col_ref(&table_cal2, "XMAX");
00206     col2_deg  = get_col_ref(&table_cal2, "DEGREE");
00207     col2_flag = get_col_ref(&table_cal2, "FLAG");
00208     col2_nrej = get_col_ref(&table_cal2, "NREJ");
00209 
00210     /* Lecture du degre maximum  et ouverture des colonnes des coefs */
00211     if (RD_desc(&table_cal2, WLMAXDEGRE, INT, 1, &maxdegre_cal2) < 0) {
00212       close_tiger_frame(&cubein);
00213       print_error("Cannot read descriptor %s in table %s",
00214                   WLMAXDEGRE,table_cal2.name);
00215       exit_session(ERR_NODESC);
00216     }
00217     col2_c    =    (int *)malloc(maxdegre_cal2*sizeof(int));
00218     coef2     = (double *)malloc(maxdegre_cal2*sizeof(double));
00219     for (i=0; i<maxdegre_cal2; i++) {
00220       sprintf(ident,"C%d", i); /* colonnes des coeff. du polynome */
00221       col2_c[i] = get_col_ref(&table_cal2, ident);
00222     }
00223   }
00224 
00225   /*=========================================================================*/
00226   /* Option de rebin en log neperien */
00227   /*=========================================================================*/
00228   if (LOG) print_msg("Logarithmic rebin");
00229   else     print_msg("Linear rebin");
00230 
00231   /*=========================================================================*/
00232   /* On recupere le descripteur LRANGE que l'on va utiliser pour lambdastart */
00233   /*=========================================================================*/
00234   if (RD_desc(&cubein,LRANGE4,DOUBLE,4,range) < 0) {
00235     print_error("Cannot read descriptor %s in datacube %s",
00236                 LRANGE4,cubein.name);
00237 
00238     print_warning("Previous versions of extract_spec (bef. 10/11/04, v.1.31)\n"
00239                   "were writing a too long keyword LBDARANGn (n=1...4).\n"
00240                   "See if you can find it back and rename it to LRANGEn.\n");
00241     exit_session(ERR_NODESC);
00242   }
00243   else {
00244     lambdaref_start = range[0];              /* Useful range */
00245     lambdaref_end   = range[1];
00246     nbpix = ((lambdaref_start - lambdaref_end) / cubein.step + 1.5);
00247     print_msg("Spectral domain: %.2f--%.2f",lambdaref_start,lambdaref_end);
00248   }
00249 
00250   /*=========================================================================*/
00251   /* On recupere le step du cube, arrondi pour ne pas avoir de problemes */
00252   /*=========================================================================*/
00253   step_in = cubein.step;
00254   step_rebin_log = (int)(1.e8 * cubein.step * 2. / 
00255                          (lambdaref_start + lambdaref_end)) / 1.e8;
00256   step_rebin_lin = cubein.step;
00257 
00258   /*=========================================================================*/
00259   /* Creation du cube de sortie rebinne */
00260   /*=========================================================================*/
00261   if (LOG) 
00262     create_tiger_frame(&cubeout, outname, -1, -1, step_rebin_log, FLOAT, 
00263                        cubein.table_name, "Rebinned cube", " ");
00264   else 
00265     create_tiger_frame(&cubeout, outname, -1, -1, step_rebin_lin, FLOAT,
00266                        cubein.table_name, "Rebinned cube", " ");
00267 
00268   min_end_rebin = MAXFLOAT;
00269   max_start_rebin = 0;
00270 
00271   /*=========================================================================*/
00272   /* On commence a iterer sur les spectres du cube d'entree */
00273   /*=========================================================================*/
00274 
00275   reset_print_progress();
00276   for (l=0; l<cubein.nbspec; l++) {
00277 
00278     /*=============================================================*/
00279     /* Initialisation pour le test du statut de la calibration */
00280     /*=============================================================*/
00281     ifail = 0;
00282     /* recherche de la ligne correspondante dans la table */
00283     if ((ligne = search_in_col(&table_cal,col_no,&(nolens[l]))) < 0) {
00284       print_warning("Lens #%d not found in the calibration table", 
00285                     nolens[l]);
00286       nb_notrebin++;
00287       continue;
00288     }
00289     if (CAL2) {
00290       /* recherche de la ligne correspondante dans la 2eme table */
00291       if ((ligne2 = search_in_col(&table_cal2,col2_no,&(nolens[l]))) < 0) {
00292         print_warning("Lens #%d not found in the 2nd calibration table",
00293                       nolens[l]);
00294         nb_notrebin++;
00295         continue;
00296       }
00297     }
00298 
00299     /*=============================================================*/
00300     /* progressssssss */
00301     /*=============================================================*/
00302     print_progress("Rebin in progress",100.*(l+1)/cubein.nbspec,1.);
00303 
00304     /*=============================================================*/
00305     /* On teste la valeur du flag */
00306     /*=============================================================*/
00307     if (RD_tbl(&table_cal, ligne, col_flag, &flag) == ERR_NODATA) {
00308       print_warning("No quality flag for lens #%d", nolens[l]);
00309       nb_notrebin++;
00310       continue;
00311     }
00312     if (flag > flagmax) {
00313       print_warning("Lens #%d not rebinned (quality criterion %d>%d)",
00314                     nolens[l],flag,flagmax);
00315       nb_notrebin++;
00316       continue;
00317     }
00318 
00319     if (CAL2) {
00320       if (RD_tbl(&table_cal2, ligne2, col2_flag, &flag) == ERR_NODATA) {
00321         print_warning("No quality flag for lens #%d (2nd table)", nolens[l]);
00322         nb_notrebin++;
00323         continue;
00324       }
00325 
00326       if (flag > flagmax) {
00327         print_warning("Lens #%d not rebinned (quality criterion %d>%d, 2nd table)",
00328                       nolens[l],flag,flagmax);
00329         nb_notrebin++;
00330         continue;
00331       }
00332     }
00333 
00334     /*=============================================================*/
00335     /* On recupere les coefs du polynome */
00336     /*=============================================================*/
00337     if (RD_tbl(&table_cal, ligne, col_nrej, &status)   == ERR_NODATA || 
00338         status == -1                                                 ||
00339         RD_tbl(&table_cal, ligne, col_xmin, &xmin_pol) == ERR_NODATA ||
00340         RD_tbl(&table_cal, ligne, col_xmax, &xmax_pol) == ERR_NODATA ||
00341         RD_tbl(&table_cal, ligne, col_deg,  &degre)    == ERR_NODATA) {
00342       print_warning("Lens #%d not rebinned (invalid calibration)",nolens[l]);
00343       nb_notrebin++;
00344       continue;
00345     }
00346     ncoef = degre + 1;
00347     for (i=0; i<ncoef; i++) {
00348       if (RD_tbl(&table_cal, ligne, col_c[i], &coef[i]) == ERR_NODATA)
00349         ifail = 1;
00350     }
00351 
00352     if (CAL2) {
00353       if (RD_tbl(&table_cal2, ligne2, col2_nrej, &status)    == ERR_NODATA || 
00354           status == -1                                                     ||
00355           RD_tbl(&table_cal2, ligne2, col2_xmin, &xmin_pol2) == ERR_NODATA ||
00356           RD_tbl(&table_cal2, ligne2, col2_xmax, &xmax_pol2) == ERR_NODATA ||
00357           RD_tbl(&table_cal2, ligne2, col2_deg,  &degre2)    == ERR_NODATA) {
00358         print_warning("Lens #%d not rebinned (invalid calibration, 2nd table)",
00359                       nolens[l]);
00360         nb_notrebin++;
00361         continue;
00362       }
00363       ncoef2 = degre2 + 1;
00364       for (i=0; i<ncoef2; i++) {
00365         if (RD_tbl(&table_cal2, ligne2, col2_c[i], &coef2[i]) == ERR_NODATA)
00366           ifail = 1;
00367       }
00368     }
00369 
00370     /*=============================================================*/
00371     /* On ouvre le spectre correspondant */
00372     /*=============================================================*/
00373     get_tiger_spec(&cubein,&specin, (NOISE? &noisein:NULL),nolens[l]);
00374     
00375 //    /* DEBUG DEBUG DEBUG */
00376 //    if (nolens[l] == 126) 
00377 //      for (i=0; i<specin.npts; i++) 
00378 //        print_msg("DEBUG DEBUG DEBUG Lens #%d: lambda=%f -> %f", 
00379 //                  nolens[l], coord_spec(&specin, i), 
00380 //                  newLambda(coord_spec(&specin, i), 
00381 //                            coef, ncoef, xmin_pol, xmax_pol, &status));
00382 //    /* DEBUG DEBUG DEBUG */
00383 
00384     /*=============================================================*/
00385     /* On teste les bornes du filtre comparees a celles du spectre d'entree */
00386     /*=============================================================*/
00387     if (lambdaref_start > specin.end || lambdaref_end < specin.start) {
00388       print_warning("Lens #%d: Spectral domain (%.2f--%.2f) not within "
00389                     "filter bounds",nolens[l],specin.start,specin.end);
00390       ifail = 1;
00391     }
00392 
00393     if (ifail) {
00394       print_msg("Lens #%d not rebinned", nolens[l]);
00395       nb_notrebin++;
00396       free_spec_mem(&specin);
00397       continue;
00398     }
00399 
00400     /*=============================================================*/
00401     /* DETERMINATION DU START DU SPECTRE REBINNE */
00402     /*=============================================================*/
00403     /* On determine d'abord le debut et fin des donnees d'entree */
00404     /*=============================================================*/
00405     xmin_in = specin.start - step_in / 2.;
00406     xmax_in = specin.end   + step_in / 2.;
00407 
00408     /*=============================================================*/
00409 
00410     /* ********** LOG ********** */
00411 
00412     /* On part de la borne inf du filtre jusqu'a rencontrer le premier pixel
00413        du spectre d'entree Puis on fait la meme chose en partant de la fin et
00414        en descendant jusqu'a trouver le premier pixel Cela nous donne le
00415        premier et dernier pixel du spectre rebinne et les indices kmin et kmax
00416        correspondants */
00417 
00418     if (LOG) {
00419 
00420       kmin = floor(log(lambdaref_start) / step_rebin_log);
00421       kmax = ceil(log(lambdaref_end) / step_rebin_log);
00422       if (kmin > kmax) {
00423         print_msg("Lens #%d not rebinned (problems in bounds)",nolens[l]);
00424         free_spec_mem(&specin);
00425         continue;
00426       }
00427 
00428       for (k=kmin; k<=kmax; k++) {
00429         lambda = exp((k -0.5) * step_rebin_log);
00430         eval_tmp = newLambda(lambda, coef, ncoef, xmin_pol, xmax_pol, 
00431                              &status);
00432 
00433         if (CAL2) {
00434           eval_tmp2 = newLambda(lambda, coef2, ncoef2, xmin_pol2, xmax_pol2, 
00435                                 &status2);
00436           eval_tmp = (eval_tmp + eval_tmp2)/2.;
00437           status = status || status2;
00438         }
00439 
00440         if (eval_tmp >= xmin_in && status == 0) {
00441           kmin = k;
00442           break;
00443         }
00444       }
00445 
00446       for (k=kmax; k>=kmin; k--) {
00447         lambda = exp((k +0.5) * step_rebin_log);
00448         eval_tmp = newLambda(lambda, coef, ncoef, xmin_pol, xmax_pol, 
00449                              &status);
00450 
00451         if (CAL2) {
00452           eval_tmp2 = newLambda(lambda, coef2, ncoef2, xmin_pol2, xmax_pol2, 
00453                                 &status2);
00454           eval_tmp = (eval_tmp + eval_tmp2)/2.;
00455           status = status || status2;
00456         }
00457 
00458         if (eval_tmp <= xmax_in && status == 0) {
00459           kmax = k;
00460           break;
00461         }
00462       }
00463 
00464       /*=============================================================*/
00465       /* on calcule le start et end du spectre rebinne (LOG) 
00466          a partir des indices trouves */
00467       /*=============================================================*/
00468       start_rebin = kmin * step_rebin_log;
00469       end_rebin   = kmax * step_rebin_log;
00470 
00471     } /* LOG */
00472 
00473     /* ********** LIN ********** */
00474 
00475     /* meme chose */
00476 
00477     else {
00478       if ( (kmin = floor(lambdaref_start / step_rebin_lin)) > 
00479            (kmax = ceil(lambdaref_end / step_rebin_lin)   ) ) {
00480         print_msg("Lens #%d not rebinned (problems in bounds)", nolens[l]);
00481         nb_notrebin++;
00482         free_spec_mem(&specin);
00483         continue;
00484       }
00485 
00486       for (k=kmin; k<=kmax; k++) {
00487         lambda = (k - 0.5) * step_rebin_lin;
00488         eval_tmp = newLambda(lambda, coef, ncoef, xmin_pol, xmax_pol, 
00489                              &status);
00490 
00491         if (CAL2) {
00492           eval_tmp2 = newLambda(lambda, coef2, ncoef2, xmin_pol2, xmax_pol2, 
00493                                 &status2);
00494           eval_tmp = (eval_tmp + eval_tmp2)/2.;
00495           status = status || status2;
00496         }
00497 
00498         if (eval_tmp >= xmin_in && status == 0) {
00499           kmin = k;
00500           break;
00501         }
00502       }
00503 
00504       for (k=kmax; k>=kmin; k--) {
00505         lambda = (k +0.5) * step_rebin_lin;
00506         eval_tmp = newLambda(lambda, coef, ncoef, xmin_pol, xmax_pol, 
00507                              &status);
00508 
00509         if (CAL2) {
00510           eval_tmp2 = newLambda(lambda, coef2, ncoef2, xmin_pol2, xmax_pol2, 
00511                                 &status2);
00512           eval_tmp = (eval_tmp + eval_tmp2)/2.;
00513           status = status || status2;
00514         }
00515 
00516         if (eval_tmp <= xmax_in && status == 0) {
00517           kmax = k;
00518           break;
00519         }
00520       }
00521 
00522       /*=============================================================*/
00523       /* on calcule le start et end du spectre rebinne (LIN) a partir des
00524          indices trouves */
00525       /*=============================================================*/
00526       start_rebin = kmin * step_rebin_lin;
00527       end_rebin   = kmax * step_rebin_lin;
00528 
00529     } /* LIN */
00530 
00531     /*=============================================================*/
00532     /* on sauvegarde le plus grand intervalle commun a tout les spectres */
00533     /*=============================================================*/
00534     max_start_rebin = MAX(max_start_rebin, start_rebin);
00535     min_end_rebin   = MIN(min_end_rebin,   end_rebin);
00536 
00537     /*=============================================================*/
00538     /* d'ou le nombre de points du spectre rebinne */
00539     /*=============================================================*/
00540     npts_rebin  = kmax - kmin + 1;
00541     value_in  = (double *)malloc(specin.npts*sizeof(double));
00542     value_out = (double *)malloc(npts_rebin *sizeof(double));
00543     if (NOISE) {
00544       value_noisein  = (double *)malloc(specin.npts*sizeof(double));
00545       value_noiseout = (double *)malloc(npts_rebin *sizeof(double));
00546     }
00547     /*=============================================================*/
00548     /* on transfere les donnees du spectre dans un tableau de doubles */
00549     /*=============================================================*/
00550     for (i=0; i<specin.npts; i++) {
00551       value_in[i]= RD_spec(&specin,i);
00552       if (NOISE) value_noisein[i]= RD_spec(&noisein,i);
00553     }
00554 
00555     /*=============================================================*/
00556     /* REBIN DU SPECTRE PAR INTERPOLATION */
00557     /*=============================================================*/
00558 
00559     for (i=0; i<npts_rebin; i++) {
00560       /*=============================================================*/
00561       /* Pixel par pixel du spectre rebinne */
00562       /*=============================================================*/
00563 
00564       /*=============================================================*/
00565       /* Bornes LAMBDA sup et inf du spectre rebinne */
00566       /*=============================================================*/
00567       if (LOG) {
00568         lambda_inf = exp(start_rebin + (i - 0.5) * step_rebin_log);
00569         lambda_sup = exp(start_rebin + (i + 0.5) * step_rebin_log);
00570       }
00571       else {
00572         lambda_inf = start_rebin + (i - 0.5) * step_rebin_lin;
00573         lambda_sup = start_rebin + (i + 0.5) * step_rebin_lin;
00574       }
00575 
00576       /*=============================================================*/
00577       /* Bornes X sup et inf du spectre rebinne */
00578       /*=============================================================*/
00579       xinf_rebin = newLambda(lambda_inf, coef, ncoef, 
00580                              xmin_pol, xmax_pol, &status);
00581       xsup_rebin = newLambda(lambda_sup, coef, ncoef, 
00582                              xmin_pol, xmax_pol, &status);
00583 
00584       if (CAL2) {
00585         xinf_rebin2 = newLambda(lambda_inf, coef2, ncoef2, 
00586                                 xmin_pol2, xmax_pol2, &status2);
00587         xsup_rebin2 = newLambda(lambda_sup, coef2, ncoef2, 
00588                                 xmin_pol2, xmax_pol2, &status2);
00589 
00590         xinf_rebin = (xinf_rebin + xinf_rebin2)/2.;
00591         xsup_rebin = (xsup_rebin + xsup_rebin2)/2.;
00592       }
00593 
00594       /*=============================================================*/
00595       /* On determine maintenant les pixels du spectre d'entree
00596          qui sont intersectes par le pixel courant du spectre rebinne */
00597       /*=============================================================*/
00598       status_start = 1;
00599       status_end   = 1;
00600       for (j=0; j<=specin.npts; j++) {
00601         xinf_pix = specin.start + (j-0.5) * step_in;
00602         xsup_pix = specin.start + (j+0.5) * step_in;
00603         if (status_start && xsup_pix >= xinf_rebin) {
00604           k_start = j;
00605           status_start = 0;
00606         }
00607         if (xinf_pix >= xsup_rebin) {
00608           k_end = j - 1;
00609           status_end = 0;
00610           break;
00611         }
00612       }
00613 
00614       /*=============================================================*/
00615       /* si on n'a pas trouve de pixel correspondant */
00616       /*=============================================================*/
00617       if (status_start || status_end) ifail = 1;
00618 
00619       /*=============================================================*/
00620       /* si le pixel du spectre rebinne courant est entierement
00621          dans un pixel du spectre d'entree */
00622       /*=============================================================*/
00623       if (k_start == k_end) {
00624         frac = (xsup_rebin - xinf_rebin) / step_in;
00625         value_out[i] = frac * value_in[k_start];
00626         if (NOISE)
00627           /* value_noiseout[i] = frac*frac * value_noisein[k_start]; */
00628           value_noiseout[i] = frac * value_noisein[k_start];
00629         continue;
00630       }
00631 
00632       /*=============================================================*/
00633       /* si il y a plus d'un pixel on calcule les fractions
00634          correspondantes aux pixels extremes sup et inf */
00635       /*=============================================================*/      
00636       frac_end  = (xsup_rebin - specin.start - (k_end - 0.5)*step_in)/step_in;
00637       frac_start= (specin.start - xinf_rebin + (k_start+0.5)*step_in)/step_in;
00638       value_out[i] = frac_start*value_in[k_start] + 
00639         frac_end*value_in[k_end];
00640       if (NOISE)
00641         /* value_noiseout[i] = frac_start*frac_start*value_noisein[k_start] + 
00642            frac_end*frac_end*value_noisein[k_end]; */
00643         value_noiseout[i] = frac_start*value_noisein[k_start] +
00644           frac_end*value_noisein[k_end];
00645 
00646       /*=============================================================*/
00647       /* et on ajoute les pixels restants en entier */
00648       /*=============================================================*/
00649       for (k=k_start+1; k<k_end; k++) {
00650         value_out[i] += value_in[k];
00651         if (NOISE) value_noiseout[i] += value_noisein[k];
00652       }
00653     }
00654 
00655     /*=============================================================*/
00656     /* on n'a pas reussi a trouver des pixels correspondants */
00657     /*=============================================================*/
00658     if (ifail) {
00659       print_msg("Lens #%d not rebinned (error in interpolation)",nolens[l]);
00660       nb_notrebin++;
00661       continue;
00662     }
00663 
00664     /*=============================================================*/
00665     /* Sauvegarde du spectre rebine */
00666     /*=============================================================*/
00667 
00668     init_new_tiger_spec(&cubeout, &specout, npts_rebin, start_rebin);
00669     if (NOISE)
00670       init_new_tiger_spec(&cubeout, &noiseout, npts_rebin, start_rebin);
00671     for (i=0; i<npts_rebin; i++) {
00672       WR_spec(&specout, i, value_out[i]);
00673       if (NOISE) WR_spec(&noiseout, i, value_noiseout[i]);
00674     }
00675     put_tiger_spec(&cubeout, &specout, (NOISE? &noiseout:NULL), nolens[l]);
00676 
00677     /*=============================================================*/
00678     /* Fermeture et liberation de la memoire allouee dans la boucle */
00679     /*=============================================================*/
00680     free(value_out);
00681     free(value_in);
00682     free_spec_mem(&specin);
00683     if (NOISE) {
00684       free(value_noisein);
00685       free(value_noiseout);
00686       free_spec_mem(&noisein);
00687     }
00688 
00689   } /* Loop over spectra */
00690 
00691   if (nb_notrebin > 0) 
00692     print_warning("%d spectra not rebinned", nb_notrebin);
00693 
00694   /*=========================================================================*/
00695   /* Transfert des descripteurs non standards */
00696   /*=========================================================================*/
00697   CP_non_std_desc(&cubein, &cubeout);
00698 
00699   /*=========================================================================*/
00700   /* Ecriture du descripteur de class */
00701   /*=========================================================================*/
00702   switch (fclass) {
00703   case RAW_OBJ_CUBE: write_file_class(&cubeout, WAV_OBJ_CUBE); break;
00704   case RAW_SKY_CUBE: write_file_class(&cubeout, WAV_SKY_CUBE); break;
00705   case RAW_DOM_CUBE: write_file_class(&cubeout, WAV_DOM_CUBE); break;
00706   case RAW_CON_CUBE: write_file_class(&cubeout, WAV_CON_CUBE); break;
00707   case RAW_CAL_CUBE: write_file_class(&cubeout, WAV_CAL_CUBE); break;
00708   case SUPER_CLASS:  write_file_class(&cubeout, SUPER_CLASS);  break;
00709   }
00710 
00711   /*=========================================================================*/
00712   /* On met a jour le descripteur LRANGE avec l'intervalle commun max */
00713   /*=========================================================================*/
00714   range[2] = max_start_rebin; 
00715   range[3] = min_end_rebin;
00716   WR_desc(&cubeout, LRANGE4, DOUBLE, 4, range);
00717 
00718   /* save meancalib error and stddev in output tiger */
00719   RD_desc(&table_cal, WLMEANERR, FLOAT, 1, &rdesc);
00720   if (CAL2) {
00721     RD_desc(&table_cal2, WLMEANERR, FLOAT, 1, &rdesc2);
00722     rdesc = (rdesc + rdesc2)/2.;
00723   }
00724   WR_desc(&cubeout, WLMEANERR, FLOAT, 1, &rdesc);
00725 
00726   RD_desc(&table_cal, WLSTDERR,  FLOAT, 1, &rdesc);
00727   if (CAL2) {
00728     RD_desc(&table_cal2, WLSTDERR, FLOAT, 1, &rdesc2);
00729     rdesc = sqrt((rdesc*rdesc + rdesc2*rdesc2)/2.);
00730   }
00731   WR_desc(&cubeout, WLSTDERR,  FLOAT, 1, &rdesc);
00732 
00733   close_tiger_frame(&cubeout);
00734 
00735   /*=========================================================================*/
00736   /* on ferme tout */
00737   /*=========================================================================*/
00738 
00739   close_tiger_frame(&cubein);
00740   close_table(&table_cal);
00741 
00742   free(nolens); free(col_c); free(coef);
00743 
00744   if (CAL2) {
00745     close_table(&table_cal2);
00746     free(col2_c); free(coef2);
00747   }
00748 
00749   exit_session(OK);
00750   return(OK);
00751 }

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