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

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