sinfo_new_wave_cal_slit2.c

00001 /*
00002  * This file is part of the ESO SINFONI Pipeline
00003  * Copyright (C) 2004,2005 European Southern Observatory
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or
00008  * (at your option) any later version.
00009  *
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00018  */
00019 /*----------------------------------------------------------------------------
00020  
00021    File name    :       sinfo_new_wave_cal_slit2.c
00022    Author       :   A. Modigliani
00023    Created on   :   Sep 17, 2003
00024    Description  : 
00025 
00026  ---------------------------------------------------------------------------*/
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #  include <config.h>
00030 #endif
00031 /*----------------------------------------------------------------------------
00032                                 Includes
00033  ---------------------------------------------------------------------------*/
00034 #include "sinfo_new_wave_cal_slit2.h"
00035 #include "sinfo_pro_save.h"
00036 #include "sinfo_pro_types.h"
00037 #include "sinfo_wavecal_ini_by_cpl.h"
00038 #include "sinfo_wcal_functions.h"
00039 #include "sinfo_absolute.h"
00040 #include "sinfo_wave_calibration.h"
00041 #include "sinfo_wavecal.h"
00042 #include "sinfo_globals.h"
00043 
00044 #include "sinfo_utilities.h"
00045 #include "sinfo_utils_wrappers.h"
00046 #include "sinfo_error.h"
00047 
00048 
00049 /*----------------------------------------------------------------------------
00050                                 Defines
00051  ---------------------------------------------------------------------------*/
00052 
00053 /*----------------------------------------------------------------------------
00054                              Function Definitions
00055  ---------------------------------------------------------------------------*/
00056 int new_dump_fit_params_to_tbl ( FitParams ** params, char * filename );
00057 
00065 /*----------------------------------------------------------------------------
00066    Function     :       sinfo_wave_cal_slit()
00067    In           :       ini_file: file name of according .ini file
00068    Out          :       integer (0 if it worked, -1 if it doesn't) 
00069    Job          :
00070 
00071 
00072   Normal method:
00073 
00074   does the wavelength calibration and the fitting of the slitlet sinfo_edge 
00075   positions (ASCII file 32 x 2 values) if wished
00076   produces an array of the bcoefs and of the fit parameters if wished and a 
00077   wavelength calibration map input is an emission line frame and a line list
00078 
00079 
00080   o searching for lines by cross sinfo_correlation with a line list
00081   o Gaussian fitting of emission lines in each column->positions of the lines->
00082     resulting fit parameters can be stored in an ASCII file
00083   o Fitting of a polynomial to the line positions for each column
00084   o Smoothing: fitting of each polynomial coefficient by another polynomial
00085     across the whole frame -> resulting polynomial coefficients can be stored 
00086     in an ASCII file.
00087   o Wavelength calibration map (micron value for each frame pixel) can be
00088     produced by using these coefficients and a cross sinfo_correlation to the
00089     original frame
00090 
00091   o The slitlet sinfo_edge positions can be fitted:
00092     1) Automatically (not really stable) or by using guess sinfo_edge positions
00093     2) By using a Boltzmann or a linear slope function
00094 
00095   Slit method:
00096 
00097   does the wavelength calibration and the fitting of the slitlet sinfo_edge 
00098   positions (ASCII file 32 x 2 values) if wished produces a list of the fit 
00099   parameters and of the smoothed coefficients if wished and a wavelength 
00100   calibration map input is an emission line frame and a line list
00101 
00102   o Does the same as other method but smoothes the found polynomial 
00103     coefficients within each slitlet and not over the whole frame.
00104 
00105   o Produces always a wavelength calibration map and does not crosscorrelate.
00106 
00107  ---------------------------------------------------------------------------*/
00108 
00109 
00110 int 
00111 sinfo_new_wave_cal_slit2(const char* plugin_id,
00112                         cpl_parameterlist* config, 
00113                         cpl_frameset* sof)
00114 {
00115   wave_config * cfg =NULL;
00116   char col_name[MAX_NAME_SIZE];
00117   char tbl_name[MAX_NAME_SIZE];
00118   char tbl_fitpar_name[MAX_NAME_SIZE];
00119   char tbl_line_list_name[MAX_NAME_SIZE];
00120   char tbl_slitpos_guess_name[MAX_NAME_SIZE];
00121   char key_name[MAX_NAME_SIZE];
00122   char col[MAX_NAME_SIZE];
00123   int check = 0;
00124   int lx = 0;
00125   int ly = 0;
00126   int n_lines=0;
00127   int i = 0;
00128   int j = 0;
00129   int n = 0;
00130   int readsum=0;
00131   int sum=0;
00132   int fit=0;
00133   int sw=0;
00134  
00135   int* status=NULL;
00136   int* n_found_lines=NULL;
00137   int* sum_pointer=NULL;
00138   int** row_clean=NULL;
00139 
00140   float a=0;
00141   float shift=0;
00142   float val_x=0;
00143   float val_y=0;
00144 
00145   float* wave=NULL;
00146   float* intens=NULL;
00147 
00148   float** acoefs=NULL;
00149   float** wavelength_clean=NULL;
00150   float** sinfo_slit_pos=NULL;
00151 
00152   double fwhm_med=0;
00153   double fwhm_avg=0;
00154   double coef_med=0;
00155   double coef_avg=0;
00156  
00157   cpl_image * im=NULL ;
00158 
00159   FitParams** par=NULL;
00160 
00161   cpl_table* tbl_wcal=NULL;
00162   cpl_table* tbl_spos=NULL;
00163   cpl_table* tbl_fitpar = NULL;
00164   cpl_table* tbl_line_list = NULL;
00165   cpl_table* tbl_slitpos_guess=NULL;
00166   cpl_table * tbl_fp =NULL;
00167   cpl_table* qclog_tbl=NULL;
00168 
00169   cpl_image* map_img=NULL;
00170 
00171   cpl_frameset* raw=NULL;
00172   cpl_parameter* p=NULL;
00173 
00174   qc_wcal* qc=sinfo_qc_wcal_new();
00175   /*   -----------------------------------------------------------------
00176        1) parse the file names and parameters to the ns_config data 
00177           structure cfg
00178        -----------------------------------------------------------------
00179   */
00180   sinfo_msg("Parsing cpl input");
00181   check_nomsg(raw=cpl_frameset_new());
00182   cknull(cfg = sinfo_parse_cpl_input_wave(config,sof,&raw),
00183      "could not parse cpl input!") ;
00184  
00185   check_nomsg(p = cpl_parameterlist_find(config,
00186                                          "sinfoni.wavecal.slitpos_boostrap"));
00187   check_nomsg(sw=cpl_parameter_get_bool(p));
00188 
00189   if(sw==1) {
00190      cfg->nslitlets=32;
00191      cfg->calibIndicator=1;
00192      cfg->wavemapInd=0;
00193      cfg->slitposIndicator=1;
00194      sinfo_msg("***********************************");
00195      sinfo_msg("parameter setting for %s",PRO_WAVE_SLITPOS_STACKED);
00196      sinfo_msg("***********************************");
00197   }
00198 
00199   if(sinfo_is_fits_file(cfg->inFrame) != 1) {
00200     sinfo_msg_error("Input file cfg->inFrame %s is not FITS",cfg->inFrame);
00201     goto cleanup;
00202   }
00203 
00204 
00205   if (cfg->slitposIndicator == 1 && cfg->estimateIndicator == 1) {
00206     if (sinfo_is_fits_file(cfg->slitposGuessName) != 1) {
00207       sinfo_msg_error("slitlet position guess list not given!");
00208       goto cleanup;
00209     }
00210   }
00211 
00212   if (cfg->calibIndicator == 0 && cfg->wavemapInd == 1) {
00213     if (sinfo_is_fits_file(cfg->coeffsName) != 1) {
00214       sinfo_msg_error("coefficients list not given!");
00215       goto cleanup;
00216     }
00217   }
00218 
00219   if (cfg->slitposIndicator == 1) {
00220     if (cfg->calibIndicator != 1 && cfg->estimateIndicator != 1) {
00221       if (sinfo_is_fits_file(cfg->paramsList) != 1) {
00222     sinfo_msg_error("parameter list not given!");
00223     goto cleanup;
00224       }
00225     }
00226   }
00227  
00228   /*---load the emission line frame---*/
00229   check(im = cpl_image_load(cfg->inFrame,CPL_TYPE_FLOAT,0,0)
00230         ,"could not load image");
00231   lx = cpl_image_get_size_x(im);
00232   ly = cpl_image_get_size_y(im);
00233 
00234   if (cfg->calibIndicator == 1 || cfg->wavemapInd == 1) {
00235     /*---open the line list and read the number of lines---*/
00236     strcpy(tbl_line_list_name,cfg->lineList);
00237     check_nomsg(tbl_line_list = cpl_table_load(tbl_line_list_name,1,0));
00238     check_nomsg(n = cpl_table_get_nrow(tbl_line_list));
00239     n_lines = n;
00240 
00241     check_nomsg(wave   = cpl_table_get_data_float(tbl_line_list,"wave"));
00242     check_nomsg(intens = cpl_table_get_data_float(tbl_line_list,"int"));
00243   }
00244 
00245   /*
00246   ----------------------------------------------------------------------
00247   ---------------------------FINDLINES----------------------------------
00248   ----------------------------------------------------------------------
00249   */
00250 
00251   /*if not yet done:
00252     do the wavelength calibration, that means: 
00253     find the dispersion relation and parameterize its coefficients
00254   */
00255   sinfo_msg("guessBeginWave=%g",cfg->guessBeginWavelength);
00256   sinfo_msg("guessDisp1=%g",cfg->guessDispersion1);
00257   sinfo_msg("guessDisp2=%g",cfg->guessDispersion2);
00258 
00259   if (cfg->calibIndicator == 1 && cfg->wavemapInd == 0) {
00260     sinfo_msg("Findlines");
00261     acoefs  = sinfo_new_2Dfloatarray(cfg->nrDispCoefficients, lx);
00262 
00263     /*allocate memory*/
00264     n_found_lines    = sinfo_new_intarray(lx); 
00265     row_clean        = sinfo_new_2Dintarray(lx, n_lines);
00266     wavelength_clean = sinfo_new_2Dfloatarray(lx, n_lines);
00267     sum_pointer      = sinfo_new_intarray(1) ;
00268 
00269     /*find the emission lines in each image column*/
00270     sinfo_new_intarray_set_value(sum_pointer, 0, 0);
00271  
00272     ck0(check = sinfo_new_find_lines(im, 
00273                                 wave, 
00274                                 intens, 
00275                                 n_lines, 
00276                                 row_clean, 
00277                                 wavelength_clean, 
00278                                 cfg->guessBeginWavelength, 
00279                         cfg->guessDispersion1, 
00280                                 cfg->guessDispersion2,
00281                                 cfg->mindiff, 
00282                                 cfg->halfWidth, 
00283                                 n_found_lines, 
00284                                 cfg->sigma, 
00285                                 sum_pointer),
00286     "sinfo_findLines failed!");
00287 
00288     /*---------------------------------------------------------------------
00289      *-----------------------WAVE_CALIB-------------------------------------
00290      *---------------------------------------------------------------------
00291     */
00292     sinfo_msg("Wave Calibration");
00293     sum = sinfo_new_intarray_get_value(sum_pointer,0);
00294     /* allocate memory for the fit parameters */
00295     cknull(par = sinfo_new_fit_params( sum ),
00296        "sinfo_newFitParams failed!");
00297 
00298   /*
00299    fit each line, make a polynomial fit and fit the resulting fit 
00300    coefficients across the columns of the slitlet
00301    */
00302    cknull(map_img = sinfo_new_wave_cal(im, 
00303                   par, 
00304                   acoefs,
00305                   cfg->nslitlets, 
00306                   row_clean,
00307                   wavelength_clean,
00308                   n_found_lines,
00309                   cfg->guessDispersion1,
00310                   cfg->halfWidth,
00311                   cfg->minAmplitude,
00312                   cfg->maxResidual,
00313                   cfg->fwhm,
00314                   cfg->nrDispCoefficients, 
00315                   cfg->nrCoefCoefficients,
00316                   cfg->sigmaFactor,
00317                   cfg->pixeldist,
00318                   cfg->pixel_tolerance),
00319       "sinfo_wave_cal failed!");
00320 
00321    sinfo_msg("Check line positions");
00322 
00323    shift=sinfo_new_check_line_positions(im,acoefs,cfg->nrDispCoefficients,
00324                     cfg->guessDispersion1,par);
00325    if (FLAG == shift){ 
00326       sinfo_msg_error("checkForLinePositions failed!\n");
00327    }
00328 
00329 
00330    sinfo_det_ncounts(raw, cfg->qc_thresh_max, qc);
00331 
00332    cknull_nomsg(qclog_tbl = sinfo_qclog_init());
00333    ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC WAVE ALL",n_lines,
00334                                  "Number of found lines","%d"));
00335    ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC WAVE NPIXSAT",qc->nsat,
00336                                  "Number of saturated pixels","%d"));
00337    ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC WAVE MAXFLUX",qc->max_di,
00338                                  "Max int off-lamp subtracted frm","%g"));
00339    ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC WAVE POSERR",shift,
00340                                  "Overall positioning error in mum","%g"));
00341 
00342    ck0(sinfo_pro_save_ima(map_img,raw,sof,cfg->outName,
00343               PRO_WAVE_MAP,qclog_tbl,plugin_id,config),
00344        "cannot save ima %s", cfg->outName);
00345 
00346    sinfo_free_table(&qclog_tbl);
00347    sinfo_free_image(&map_img);
00348 
00349    /*
00350     #store the resulting polynomial fit coefficients in an 
00351      ASCII file if wished
00352    */
00353 
00354    if (cfg->writeCoeffsInd == 1) {
00355      check_nomsg(tbl_wcal = cpl_table_new(lx));
00356      for (i=0; i< cfg->nrDispCoefficients; i++) {
00357        snprintf(col_name,MAX_NAME_SIZE-1,"%s%d","coeff",i);
00358        check_nomsg(cpl_table_new_column(tbl_wcal,col_name, CPL_TYPE_DOUBLE));
00359      }
00360 
00361      cknull_nomsg(qclog_tbl = sinfo_qclog_init());
00362      ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC WAVE ALL",n_lines,
00363                                    "Number found lines","%d"));
00364 
00365      for (j=0; j< lx; j++) {
00366        for (i=0; i< cfg->nrDispCoefficients; i++) {
00367      snprintf(col_name,MAX_NAME_SIZE-1,"%s%d","coeff",i);
00368      a = sinfo_new_array2D_get_value(acoefs, i, j);
00369      /* fprintf(acoefs_file, "%15.13g ", a) ; */
00370      cpl_table_set_double(tbl_wcal,col_name,j,a);
00371        }
00372      }
00373      for (i=0; i< cfg->nrDispCoefficients; i++) {
00374        snprintf(col_name,MAX_NAME_SIZE-1,"%s%d","coeff",i);
00375        check_nomsg(coef_avg=cpl_table_get_column_mean(tbl_wcal,col_name));
00376        check_nomsg(coef_med=cpl_table_get_column_median(tbl_wcal,col_name));
00377 
00378        snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%s","QC COEF",i," AVG");
00379        ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,coef_avg,
00380                                         "Average wavecal Coef","%g"));
00381 
00382        snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%s","QC COEF",i," MED");
00383        ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,coef_med,
00384                                         "Median wavecal Coef","%g"));
00385 
00386      }
00387 
00388      strcpy(tbl_name,cfg->coeffsName);
00389      ck0(sinfo_pro_save_tbl(tbl_wcal,raw,sof,tbl_name,
00390                 PRO_WAVE_COEF_SLIT,qclog_tbl,plugin_id,config),
00391      "cannot save tbl %s", tbl_name);
00392      sinfo_free_table(&tbl_wcal);
00393      sinfo_free_table(&qclog_tbl);
00394 
00395    }
00396 
00397    /*
00398     #store the resulting Gaussian fit parameters in an ASCII file if wished
00399    */
00400    if (cfg->writeParInd == 1) {
00401       sinfo_new_dump_fit_params_to_ascii(par, cfg->paramsList);
00402  
00403       cknull(par,"no fit parameters available!") ;
00404       cknull(cfg->paramsList,"no filename available!") ;
00405       check_nomsg(tbl_fp = cpl_table_new(par[0] -> n_params));
00406       check_nomsg(cpl_table_new_column(tbl_fp,"n_params", CPL_TYPE_INT));
00407       check_nomsg(cpl_table_new_column(tbl_fp,"column", CPL_TYPE_INT));
00408       check_nomsg(cpl_table_new_column(tbl_fp,"line", CPL_TYPE_INT));
00409  
00410       for(j=0;j<4;j++) {
00411          snprintf(col,MAX_NAME_SIZE-1,"%s%d","fpar",j);
00412          cpl_table_new_column(tbl_fp,col, CPL_TYPE_DOUBLE);
00413          snprintf(col,MAX_NAME_SIZE-1,"%s%d","dpar",j);
00414          cpl_table_new_column(tbl_fp,col, CPL_TYPE_DOUBLE);
00415       }
00416 
00417       cknull_nomsg(qclog_tbl = sinfo_qclog_init());
00418       ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC NLINES",n_lines,
00419                                     "Number of found lines","%d"));
00420 
00421      for ( i = 0 ; i < par[0] -> n_params ; i++ ) {
00422       
00423          check_nomsg(cpl_table_set_int(tbl_fp,"n_params",i,par[i]->n_params));
00424          check_nomsg(cpl_table_set_int(tbl_fp,"column",i,par[i]->column));
00425          check_nomsg(cpl_table_set_int(tbl_fp,"line",i,par[i]->line));
00426 
00427          for(j=0;j<4;j++) {
00428        snprintf(col,MAX_NAME_SIZE-1,"%s%d","fpar",j);
00429        if(isnan(par[i]->fit_par[j])) {
00430          cpl_table_set_invalid(tbl_fp,col,i);
00431        } else {
00432          cpl_table_set_double(tbl_fp,col,i,par[i]->fit_par[j]);
00433        }
00434        snprintf(col,MAX_NAME_SIZE-1,"%s%d","dpar",j);
00435        if(isnan(par[i]->derv_par[j])) {
00436          cpl_table_set_invalid(tbl_fp,col,i);
00437        } else {
00438          cpl_table_set_double(tbl_fp,col,i,par[i]->derv_par[j]);
00439        }
00440      }
00441       }
00442 
00443       check_nomsg(fwhm_avg = cpl_table_get_column_mean(tbl_fp,"fpar1"));
00444       check_nomsg(fwhm_med = cpl_table_get_column_median(tbl_fp,"fpar1"));
00445       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC FWHM MED",fwhm_med,
00446                                        "Median FWHM of found lines","%f"));
00447       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC FWHM AVG",fwhm_avg,
00448                                        "Average FWHM of found lines","%f"));
00449 
00450       ck0(sinfo_pro_save_tbl(tbl_fp,raw,sof,cfg->paramsList,
00451                  PRO_WAVE_PAR_LIST,qclog_tbl,plugin_id,config),
00452       "cannot save tbl %s", cfg->paramsList);
00453 
00454       sinfo_free_table(&qclog_tbl);
00455       sinfo_free_table(&tbl_fp) ;
00456 
00457    }
00458 
00459    /* free memory */
00460    sinfo_new_destroy_2Dfloatarray ( &wavelength_clean, lx );
00461    sinfo_new_destroy_2Dintarray (&row_clean, lx);
00462    sinfo_new_destroy_intarray(&n_found_lines );
00463    sinfo_new_destroy_intarray(&sum_pointer );
00464    sinfo_new_destroy_2Dfloatarray ( &acoefs, cfg->nrDispCoefficients );
00465    sinfo_free_table(&tbl_line_list);
00466 
00467    /*----------------------------------------------------------------------
00468     *-------------------WAVEMAP--------------------------------------------
00469     *----------------------------------------------------------------------
00470     */
00471 
00472    /*
00473     #---now do the cross sinfo_correlation and produce a wavelength map---
00474    */
00475 
00476   } else if (cfg->wavemapInd == 1 && cfg->calibIndicator == 0) { 
00477     sinfo_msg("Wavemap");
00478     acoefs = sinfo_new_2Dfloatarray ( cfg->nrDispCoefficients, lx);
00479     /* #read the parameterized dispersion relation */
00480 
00481     strcpy(tbl_name,cfg->coeffsName);
00482     check_nomsg(tbl_wcal = cpl_table_load(tbl_name,1,0));
00483 
00484     for (i =0; i < lx; i++) {
00485       for (j = 0; j< cfg->nrDispCoefficients; j++) {
00486     snprintf(col_name,MAX_NAME_SIZE-1,"%s%d","coeff",j);
00487     acoefs[j][i]=cpl_table_get_double(tbl_wcal,col_name,i,status);
00488       }
00489     }
00490     sinfo_free_table(&tbl_wcal);
00491 
00492     cknull(map_img = sinfo_new_create_shifted_slit_wavemap2(im,
00493                            acoefs,
00494                            cfg->nrDispCoefficients,
00495                            wave,
00496                            intens,
00497                            n_lines,
00498                            cfg->magFactor,
00499                            cfg->guessDispersion1,
00500                            cfg->pixeldist ),
00501            "sinfo_createShiftedSlitWavemap2 failed!");
00502 
00503     par = sinfo_new_fit_params(15*n_lines);
00504     sinfo_msg("Check shifts");
00505 
00506     shift = sinfo_new_check_correlated_line_positions ( im, acoefs, 
00507                          cfg->nrDispCoefficients, 
00508                          wave,
00509                          intens,
00510                          n_lines,
00511                          cfg->fwhm,
00512                          cfg->halfWidth,
00513                          cfg->minAmplitude,
00514                          cfg->guessDispersion1,
00515                          par );
00516 
00517     if (FLAG == shift){
00518       sinfo_msg_error("sinfo_checkCorrelatedLinePositions failed!\n");
00519     }
00520 
00521     sinfo_det_ncounts(raw, cfg->qc_thresh_max,qc);
00522     cknull_nomsg(qclog_tbl = sinfo_qclog_init());
00523 
00524     ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC NLINES",n_lines,
00525                                   "Number of found lines","%d"));
00526 
00527     check_nomsg(fwhm_avg = cpl_table_get_column_mean(tbl_fp,"fpar1"));
00528     check_nomsg(fwhm_med = cpl_table_get_column_median(tbl_fp,"fpar1"));
00529 
00530     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC FWHM MED",fwhm_med,
00531                                   "Median FWHM of found lines","%f"));
00532     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC FWHM AVG",fwhm_avg,
00533                                   "Average FWHM of found lines","%f"));
00534 
00535     ck0(sinfo_pro_save_ima(map_img,raw,sof,cfg->outName,
00536                PRO_WAVE_MAP,qclog_tbl,plugin_id,config),
00537     "cannot save ima %s", cfg->outName);
00538 
00539     sinfo_free_table(&qclog_tbl);
00540     sinfo_free_image(&map_img);
00541 
00542     /* # ---free memory--- */
00543     sinfo_new_destroy_2Dfloatarray ( &acoefs, cfg->nrDispCoefficients );
00544     /* To fix a memory bug we comment the following. But check! */
00545     /* cpl_free ( wave ); */
00546     /* cpl_free ( intens ); */
00547 
00548   } else if (cfg->wavemapInd == 1 && cfg->calibIndicator == 1) {
00549     sinfo_msg_error("give either wavemapIndicator = yes and calibIndicator");
00550     sinfo_msg_error("= no or wavemapIndicator = no and calibIndicator = yes") ;
00551     goto cleanup;
00552   }
00553 
00554   /*-------------------------------------------------------------------
00555    *-------------------SLITFITS----------------------------------------
00556    *-------------------------------------------------------------------
00557    #--fit the slitlet sinfo_edge positions if desired--
00558   */
00559   if (cfg->slitposIndicator == 1) {
00560   sinfo_msg("fit the slitlet sinfo_edge positions");
00561   if (cfg->calibIndicator != 1 && cfg->estimateIndicator != 1) {
00562 
00563     /*
00564       #read the first integer to determine the number of fit 
00565       parameters to allocate
00566     */
00567 
00568      if(sinfo_is_fits_file(cfg->paramsList) !=1 ) {
00569        sinfo_msg_error("cannot read FITS file %s ",cfg->paramsList);
00570        goto cleanup;
00571      }
00572      strcpy(tbl_fitpar_name,cfg->paramsList);
00573      check_nomsg(tbl_fitpar = cpl_table_load(tbl_fitpar_name,1,0));
00574 
00575      check_nomsg(readsum=cpl_table_get_int(tbl_fitpar,"n_params",1,status));
00576      sinfo_free_table(&tbl_fitpar);
00577  
00578      cknull(sinfo_new_fit_params( readsum ),"sinfo_new_fit_params failed!");
00579 
00580      ck0(sinfo_dumpTblToFitParams(par, cfg->paramsList),
00581      "reading tbl %s ", cfg->paramsList);
00582   }
00583 
00584   /* #allocate memory for the slitlet position array */
00585   sinfo_slit_pos = sinfo_new_2Dfloatarray(32,2);
00586   /* #now decide which fit model function you want to use by 
00587      reading a given character
00588   */
00589 
00590   sinfo_msg("cfg->fitBoltzIndicator=%d\n",cfg->fitBoltzIndicator);
00591   sinfo_msg("cfg->estimateIndicator=%d\n",cfg->estimateIndicator);
00592   sinfo_msg("cfg->fitEdgeIndicator=%d\n",cfg->fitEdgeIndicator);
00593 
00594   if (cfg->fitBoltzIndicator == 1) {
00595     if (cfg->estimateIndicator == 1) {
00596       /* #open the ASCII list of the slitlet positions--- */
00597       /*READ TFITS TABLE*/
00598       strcpy(tbl_slitpos_guess_name,cfg->slitposGuessName);
00599       check_nomsg(tbl_slitpos_guess=cpl_table_load(tbl_slitpos_guess_name,1,0));
00600       check_nomsg(n = cpl_table_get_nrow(tbl_slitpos_guess));
00601 
00602       for (i =0 ; i< 32; i++){
00603         val_x=cpl_table_get_double(tbl_slitpos_guess,"pos1",i,status);
00604         val_y=cpl_table_get_double(tbl_slitpos_guess,"pos2",i,status);
00605     sinfo_new_array2D_set_value(sinfo_slit_pos,val_x,i,0);
00606     sinfo_new_array2D_set_value(sinfo_slit_pos,val_y,i,1);
00607       }
00608       sinfo_free_table(&tbl_slitpos_guess);
00609 
00610       sinfo_msg("sinfo_new_fit_slits_boltz_with_estimate");
00611       fit = sinfo_new_fit_slits_boltz_with_estimate(im, 
00612                                             sinfo_slit_pos, 
00613                                 cfg->boxLength, 
00614                                             cfg->yBox, 
00615                                             cfg->diffTol, 
00616                                             cfg->loPos, 
00617                                             cfg->hiPos );
00618       if (fit < 0){
00619     sinfo_msg_error("sinfo_new_fit_slits_boltz_with_estimate failed" );
00620     goto cleanup;
00621       }
00622     } else {
00623     sinfo_msg("sinfo_new_fit_slits_boltz");
00624     fit = sinfo_new_fit_slits_boltz(im, 
00625                                   par, 
00626                                   sinfo_slit_pos, 
00627                                   cfg->boxLength, 
00628                                   cfg->yBox, 
00629                                   cfg->diffTol );
00630 
00631     if (fit < 0) {
00632       sinfo_msg_error ( "sinfo_new_fit_slits_boltz failed" );
00633       goto cleanup;
00634     }
00635     }
00636   } else if (cfg->fitEdgeIndicator == 1) {
00637 
00638     if (cfg->estimateIndicator == 1){
00639       /*READ TFITS TABLE*/
00640       strcpy(tbl_slitpos_guess_name,cfg->slitposGuessName);
00641       check_nomsg(tbl_slitpos_guess=cpl_table_load(tbl_slitpos_guess_name,1,0));
00642       check_nomsg(n = cpl_table_get_nrow(tbl_slitpos_guess));
00643 
00644       for (i =0 ; i< 32; i++){
00645     val_x=cpl_table_get_double(tbl_slitpos_guess,"pos1",i,status);
00646     val_y=cpl_table_get_double(tbl_slitpos_guess,"pos2",i,status);
00647     sinfo_new_array2D_set_value(sinfo_slit_pos,val_x,i,0);
00648     sinfo_new_array2D_set_value(sinfo_slit_pos,val_y,i,1);
00649       }
00650       cpl_table_delete(tbl_slitpos_guess);
00651 
00652       sinfo_msg("sinfo_new_fit_slits_edge_with_estimate");
00653       fit = sinfo_new_fit_slits_edge_with_estimate(im, 
00654                                            sinfo_slit_pos,
00655                                            cfg->boxLength, 
00656                                            cfg->yBox, 
00657                                            cfg->diffTol, 
00658                                            cfg->loPos, 
00659                                cfg->hiPos );
00660       if (fit < 0) {
00661      sinfo_msg_error( "sinfo_new_fit_slits_boltz failed" );
00662          goto cleanup;
00663       }
00664     } else {
00665       sinfo_msg("sinfo_new_fit_slits_edge");
00666       fit = sinfo_new_fit_slits_edge(im, 
00667                                par, 
00668                                sinfo_slit_pos, 
00669                                cfg->boxLength, 
00670                                cfg->yBox, 
00671                                cfg->diffTol );
00672       if (fit < 0) {
00673     sinfo_msg_error("sinfo_new_fit_slits_edge failed" );
00674         goto cleanup;
00675       }
00676     }
00677   } else {
00678     sinfo_msg_error("no indication of desired fit function given" );
00679     goto cleanup;
00680   }
00681   sinfo_free_image(&im);
00682 
00683   /* #store the resulting sitlet positions in an TFITS table */
00684   check_nomsg(tbl_spos = cpl_table_new(32));
00685   check_nomsg(cpl_table_new_column(tbl_spos,"pos1", CPL_TYPE_DOUBLE));
00686   check_nomsg(cpl_table_new_column(tbl_spos,"pos2", CPL_TYPE_DOUBLE));
00687   check_nomsg(cpl_table_set_column_format(tbl_spos,"pos1", "15.9f"));
00688   check_nomsg(cpl_table_set_column_format(tbl_spos,"pos2", "15.9f"));
00689 
00690   for (i =0; i< 32; i++) {
00691     cpl_table_set_double(tbl_spos,"pos1",i,
00692                          sinfo_new_array2D_get_value(sinfo_slit_pos,i,0));
00693     cpl_table_set_double(tbl_spos,"pos2",i,
00694                          sinfo_new_array2D_get_value(sinfo_slit_pos,i,1));
00695      
00696   }
00697   if(sw == 1) {
00698     strcpy(tbl_name,"out_slitpos_guess.fits");
00699     ck0(sinfo_pro_save_tbl(tbl_spos,raw,sof,tbl_name,
00700                PRO_SLIT_POS_GUESS,NULL,plugin_id,config),
00701     "cannot save tbl %s", tbl_name);
00702   } else {
00703     strcpy(tbl_name,cfg->slitposName);
00704     ck0(sinfo_pro_save_tbl(tbl_spos,raw,sof,tbl_name,
00705                PRO_SLIT_POS,NULL,plugin_id,config),
00706     "cannot save tbl %s", tbl_name);
00707   }
00708   sinfo_free_table(&tbl_spos);
00709   sinfo_new_destroy_2Dfloatarray ( &sinfo_slit_pos, 32 );
00710   }
00711 
00712   if ( (cfg->slitposIndicator == 1 && cfg->estimateIndicator != 1) || 
00713        (cfg->calibIndicator == 1)  || (cfg->wavemapInd == 1) ){
00714     sinfo_new_destroy_fit_params(&par);
00715   }
00716 
00717   sinfo_wavecal_free(&cfg);
00718   sinfo_free_frameset(&raw);
00719   sinfo_qc_wcal_delete(&qc);
00720 
00721   return 0;
00722 
00723   cleanup:
00724   sinfo_free_table(&tbl_spos);
00725   sinfo_free_table(&tbl_fitpar);
00726   sinfo_free_image(&map_img);
00727   sinfo_free_table(&tbl_wcal);
00728   sinfo_free_table(&tbl_fp) ;
00729   sinfo_free_table(&tbl_line_list);
00730   sinfo_free_table(&tbl_wcal);
00731   sinfo_free_table(&qclog_tbl);
00732   sinfo_free_image(&map_img);
00733   sinfo_new_destroy_fit_params(&par);
00734   sinfo_new_destroy_2Dfloatarray ( &sinfo_slit_pos, 32 );
00735   sinfo_new_destroy_2Dfloatarray ( &wavelength_clean, lx );
00736   sinfo_new_destroy_2Dintarray (&row_clean, lx);
00737   sinfo_new_destroy_intarray(&n_found_lines );
00738   sinfo_new_destroy_intarray(&sum_pointer );
00739   if(acoefs!=NULL) {
00740      sinfo_new_destroy_2Dfloatarray(&acoefs,cfg->nrDispCoefficients);
00741   }
00742   sinfo_free_table(&tbl_line_list);
00743   sinfo_free_image(&im);
00744   sinfo_wavecal_free(&cfg);
00745   sinfo_free_frameset(&raw);
00746   sinfo_qc_wcal_delete(&qc);
00747   return -1 ;
00748 
00749 }

Generated on Wed Jan 17 08:33:43 2007 for SINFONI Pipeline Reference Manual by  doxygen 1.4.4