~ubuntu-branches/debian/sid/cpl-plugin-sinfo/sid

/*
 * This file is part of the ESO SINFONI Pipeline
 * Copyright (C) 2004,2005 European Southern Observatory
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
 */
/*----------------------------------------------------------------------------

   File name    :       sinfo_find_distortions.c
   Author       :    A. Modigliani
   Created on   :    Aug 12, 2004
   Description  :

 ---------------------------------------------------------------------------*/
#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif
#include <sinfo_cpl_size.h>

/*----------------------------------------------------------------------------
                                Includes
 ---------------------------------------------------------------------------*/

#include "sinfo_new_find_distortions.h"
#include "sinfo_pro_save.h"
#include "sinfo_pro_types.h"
#include "sinfo_finddist_ini_by_cpl.h"
#include "sinfo_wave_calibration.h"
#include "sinfo_cube_construct.h"
#include "sinfo_absolute.h"
#include "sinfo_distortion.h"
#include "sinfo_utilities.h"
#include "sinfo_recipes.h"
#include "sinfo_utils_wrappers.h"
#include "sinfo_error.h"
#include "sinfo_globals.h"

/*----------------------------------------------------------------------------
                                Defines
 ---------------------------------------------------------------------------*/
#define SINFO_ARC_SATURATION                 100000
#define SINFO_ARC_MAX_WIDTH                  64


/**@{*/
/**
 * @addtogroup sinfo_rec_mflat Master flat generation
 *
 * TBD
 */


/*----------------------------------------------------------------------------
                             Function Definitions
 ---------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------
   Function     :       sinfo_find_distortions()
   In           :
   Out          :
   Job          :

 ---------------------------------------------------------------------------*/
static double
new_compute_shift(double x,
                  double y,
                  double c00,
                  double c10,
                  double c01,
                  double c11,
                  double c20,
                  double c02,
                  double c21,
                  double c12,
                  double c30,
                  double c03);

int
sinfo_new_find_distortions(const char* plugin_id,
                           cpl_parameterlist* config,
                           cpl_frameset* sof,
                           cpl_frameset* set_fibre_ns)
{

    finddist_config * cfg=NULL ;
    cpl_image * imonind=NULL ;
    cpl_image * impoly=NULL ;
    cpl_image * im=NULL ;
    cpl_image *  map=NULL ;
    cpl_image *  mask=NULL ;

    int*            degx=NULL;
    int*            degy=NULL;
    int* n_found_lines=NULL;
    int* sum_pointer=NULL;
    int** row_clean=NULL;
    int pdensity=0;


    int x_l=SIZEX/4*1;
    int x_u=SIZEX/4*3;
    int y_l=SIZEY/4*1;
    int y_u=SIZEY/4*3;
    int x_c=SIZEX/4*2;
    int y_c=SIZEY/4*2;
    const int pdx=4;
    const int pdy=4;
    int check = 0;
    int i = 0;
    int lx=0;
    //int ly=0;


    int sum=0;
    int n_lines=0;
    int fit=0;
    int n=0;
    int j=0;
    cpl_size coef_pow[2];
    int arcs_window_size=0;

    float value=0;
    float newval=0;
    float total_dist=0;
    float shift=0;

    float* distances=NULL;
    float** wavelength_clean=NULL;
    float** slit_pos=NULL;
    float** acoefs=NULL;
    float* wave=NULL;
    float* intens=NULL;
    float* first =NULL;

    double*         coef=NULL;

    double xshift=0.;
    double arcs_thres_factor=0;

    double arcs_min_arclen_factor=0;
    double pcf[pdx][pdy];


    char tbl_name[FILE_NAME_SZ];
    char key_name[FILE_NAME_SZ];

    cpl_table* poly_tbl=NULL;
    cpl_table* tbl_spos=NULL;
    cpl_table* tbl_line_list=NULL;

    FitParams** par=NULL;

    cpl_frameset* stk=NULL;
    cpl_table* qclog_tbl=NULL;

    cpl_polynomial* distor_poly=NULL;
    cpl_apertures       *   arcs=NULL ;
    cpl_parameter* p=NULL;
    int smooth_rad=0;


    check_nomsg(p=cpl_parameterlist_find(config,"sinfoni.product.density"));
    check_nomsg(pdensity=cpl_parameter_get_int(p));

    /*
    parse the file names and parameters to the finddist_config
    data structure cfg
     */
    check_nomsg(stk=cpl_frameset_new());
    cknull(cfg = sinfo_parse_cpl_input_finddist(config,sof,&stk),
           "could not parse CPL input!");

    if(sinfo_is_fits_file (cfg->inFrame) != 1) {
        sinfo_msg_error("Input file %s is not FITS",cfg->inFrame);
        goto cleanup;
    }

    /* load the emission line frame--- */
    check(im = cpl_image_load(cfg->inFrame,CPL_TYPE_FLOAT,0,0),
          "could not load arc image");

    /* load the fake on - fake off frame--- */
    check(imonind = cpl_image_load(cfg->nsFrame,CPL_TYPE_FLOAT,0,0),
          "could not load on-off fake image");

    check(impoly  = cpl_image_load(cfg->nsFrame,CPL_TYPE_FLOAT,0,0),
          "could not load on-off fake image");

    /* load the fake on - fake off frame--- */
    check(mask = cpl_image_load(cfg->mask,CPL_TYPE_FLOAT,0,0),
          "could not load mask image");

    check(cpl_image_multiply (im,mask),
          "Failing to correct arc ima by bp mask");

    check(cpl_image_multiply (imonind,mask),
          "Failing to correct on-off fake ima by bp mask");

    check(cpl_image_multiply(impoly,mask),
          "Failing to correct on-off fake ima by bp mask");

    sinfo_free_image(&mask);

    check_nomsg(lx = cpl_image_get_size_x(im));
    //check_nomsg(ly = cpl_image_get_size_y(im));
    /* TO BE CHENGED THE FOLLOWING INPUT */
    /* open the line list and read the number of lines */
    check(tbl_line_list = cpl_table_load(cfg->lineList,1,0),
          "problems loading table %s",cfg->lineList);

    check_nomsg(n_lines = cpl_table_get_nrow(tbl_line_list));
    check_nomsg(wave=cpl_table_get_data_float(tbl_line_list,"wave"));
    cpl_type type=cpl_table_get_column_type(tbl_line_list,"int");
    if(type==CPL_TYPE_INT) {
        check_nomsg(cpl_table_cast_column(tbl_line_list,"int", "fint",CPL_TYPE_FLOAT));
        check_nomsg(intens = cpl_table_get_data_float(tbl_line_list,"fint"));
    } else {
        check_nomsg(intens = cpl_table_get_data_float(tbl_line_list,"int"));
    }
    /*
       #---------------------------------------------------------------------
       #---------------------------FINDLINES---------------------------------
       #---------------------------------------------------------------------
     */

    sinfo_msg("Find Lines");
    /*
       do the wavelength calibration, that means:
       find the dispersion relation and parameterize its coefficients
     */

    acoefs  = sinfo_new_2Dfloatarray(cfg->nrDispCoefficients, lx);
    /* allocate memory */
    n_found_lines    = sinfo_new_intarray(lx);
    row_clean        = sinfo_new_2Dintarray(lx, n_lines);
    wavelength_clean = sinfo_new_2Dfloatarray(lx, n_lines);
    sum_pointer      = sinfo_new_intarray(1);

    /* find the emission lines in each image column */
    sinfo_new_intarray_set_value(sum_pointer, 0, 0);

    ck0((check = sinfo_new_find_lines(im,
                    wave,
                    intens,
                    n_lines,
                    row_clean,
                    wavelength_clean,
                    cfg->guessBeginWavelength,
                    cfg->guessDispersion1,
                    cfg->guessDispersion2,
                    cfg->mindiff,
                    cfg->halfWidth,
                    n_found_lines,
                    cfg->sigma,
                    sum_pointer )),
        "FindLines failed!");

    sinfo_free_table(&tbl_line_list);

    /*
   #-------------------------------------------------------------------------
   #---------------------------WAVECALIB-------------------------------------
   #-------------------------------------------------------------------------
     */
    sinfo_msg("Do wave calib");
    sum = sinfo_new_intarray_get_value(sum_pointer,0);

    /* #allocate memory for the fit parameters */
    cknull(par = sinfo_new_fit_params( sum ),
           "sinfo_new_fit_params failed!");

    /*
      fit each line, make a polynomial fit and fit the resulting fit
       coefficients across the columns of the slitlet
     */
    cknull(map = sinfo_new_wave_cal(im,
                    par,
                    acoefs,
                    cfg->nslitlets,
                    row_clean,
                    wavelength_clean,
                    n_found_lines,
                    cfg->guessDispersion1,
                    cfg->halfWidth,
                    cfg->minAmplitude,
                    cfg->maxResidual,
                    cfg->fwhm,
                    cfg->nrDispCoefficients,
                    cfg->nrCoefCoefficients,
                    cfg->sigmaFactor,
                    cfg->pixeldist,
                    cfg->pixel_tolerance),
           "sinfo_waveCal failed!");
    /* here mem leak */

    sinfo_msg("Check line positions");
    shift = sinfo_new_check_line_positions (im, acoefs,
                    cfg->nrDispCoefficients,
                    cfg->guessDispersion1,par);

    if (FLAG == shift) {
        sinfo_msg_error ( "checkForLinePositions failed!");
        goto cleanup;
    }
    sinfo_free_image(&map);


    /* free memory */
    sinfo_new_destroy_2Dfloatarray ( &wavelength_clean, lx );
    sinfo_new_destroy_2Dintarray (&row_clean, lx);
    sinfo_new_destroy_intarray(&n_found_lines );
    sinfo_new_destroy_intarray(&sum_pointer );
    sinfo_new_destroy_2Dfloatarray ( &acoefs, cfg->nrDispCoefficients );
    /*
  sinfo_new_destroy_array (& wave );
  sinfo_new_destroy_array (& intens );
     */
    /*
  #----------------------------------------------------------------------------
  #-------------------------SLITFITS-------------------------------------------
  #----------------------------------------------------------------------------
  #--fit the slitlet sinfo_edge positions if desired--
     */

    sinfo_msg("Find slit pos");
    /*allocate memory for the slitlet position array */
    slit_pos = sinfo_new_2Dfloatarray(32,2);

    fit = sinfo_new_fit_slits_boltz( im, par, slit_pos, cfg->boxLength,
                    cfg->yBox, cfg->diffTol );


    if (fit < 0) {
        sinfo_msg_error( "sinfo_fitSlitsBoltz failed" );
        goto cleanup;
    }
    sinfo_new_destroy_fit_params(&par);

    check_nomsg(tbl_spos=cpl_table_new(32));
    check_nomsg(cpl_table_new_column(tbl_spos,"pos1", CPL_TYPE_DOUBLE));
    check_nomsg(cpl_table_new_column(tbl_spos,"pos2", CPL_TYPE_DOUBLE));
    check_nomsg(cpl_table_set_column_format(tbl_spos,"pos1", "15.9f"));
    check_nomsg(cpl_table_set_column_format(tbl_spos,"pos2", "15.9f"));

    for (i =0; i< 32; i++) {
        check_nomsg(cpl_table_set_double(tbl_spos,"pos1",i,
                        sinfo_new_array2D_get_value(slit_pos,i,0)));
        check_nomsg(cpl_table_set_double(tbl_spos,"pos2",i,
                        sinfo_new_array2D_get_value(slit_pos,i,1)));

    }

    if(pdensity > 2) {
        strcpy(tbl_name,"out_slitlets_pos_predist.fits");
        ck0(sinfo_pro_save_tbl(tbl_spos,set_fibre_ns,sof,tbl_name,
                        PRO_SLITLETS_POS_PREDIST,NULL,plugin_id,config),
            "cannot save tbl %s", tbl_name);
    }
    sinfo_free_table(&tbl_spos);

    /*
  #---------------------------------------------------------
  # do the north - south - test
  #---------------------------------------------------------
     */
    sinfo_msg("Do north south test");
    /*
  sinfo_msg("cfg->nslits =%d\n", cfg->nslits);
  sinfo_msg("cfg->nshalfWidth =%d\n", cfg->nshalfWidth);
  sinfo_msg("cfg->nsfwhm=%f\n",cfg->nsfwhm );
  sinfo_msg("cfg->minDiff=%f\n", cfg->minDiff);
  sinfo_msg("cfg->estimated_dist=%f\n", cfg->estimated_dist);
  sinfo_msg("cfg->devtol=%f\n", cfg->devtol);
  sinfo_msg("cfg->loPos=%d\n", cfg->loPos);
  sinfo_msg("cfg->hiPos=%d\n",cfg->hiPos);
     */
    cknull(distances = sinfo_north_south_test(imonind,
                    cfg->nslits,
                    cfg->nshalfWidth,
                    cfg->nsfwhm ,
                    cfg->minDiff,
                    cfg->estimated_dist,
                    cfg->devtol,
                    cfg->loPos,
                    cfg->hiPos),
           "north_south_test failed");

    sinfo_free_image(&imonind);

    /*
  #---------------------------------------------------------
  # get firstcol
  #---------------------------------------------------------
     */
    sinfo_msg("get first col");
    first = sinfo_new_floatarray(61);
    total_dist=0;
    n=0;

    for (i=0; i< 31; i++) {
        total_dist=total_dist+sinfo_new_array_get_value(distances,i);

        for (j=0; j< 2; j++) {
            if (j == 0) {
                if (i != 30) {
                    newval = sinfo_new_array2D_get_value(slit_pos,i+1,j) -
                                    total_dist;
                    sinfo_new_array_set_value(first, newval, n);
                    n = n+1;
                }
            }
            if (j == 1) {
                newval = sinfo_new_array2D_get_value(slit_pos,i,j) - total_dist;
                sinfo_new_array_set_value(first, newval, n);
                n = n+1;
            }
        }
    }
    value = sinfo_new_f_median(first,61);
    sinfo_msg("Firstcol =%f", value);
    sinfo_new_destroy_array(&first);
    sinfo_new_destroy_array(&distances);
    sinfo_new_destroy_2Dfloatarray ( &slit_pos, 32 );

    /*
  #---------------------------------------------------------
  # Determine distortions
  #---------------------------------------------------------
     */
    sinfo_msg("Determine distortions");
    degx=cpl_calloc(8,sizeof(int));
    degy=cpl_calloc(8,sizeof(int));
    coef=cpl_calloc(8,sizeof(double));

    check_nomsg(p=cpl_parameterlist_find(config,
                    "sinfoni.distortion.arcs_thresh_factor"));
    check_nomsg(arcs_thres_factor=cpl_parameter_get_double(p));

    check_nomsg(p=cpl_parameterlist_find(config,
                    "sinfoni.distortion.arcs_min_arclen_factor"));
    check_nomsg(arcs_min_arclen_factor=cpl_parameter_get_double(p));

    check_nomsg(p=cpl_parameterlist_find(config,
                    "sinfoni.distortion.arcs_window_size"));
    check_nomsg(arcs_window_size=cpl_parameter_get_int(p));

    check_nomsg(p=cpl_parameterlist_find(config,
                    "sinfoni.distortion.smooth_rad"));
    check_nomsg(smooth_rad=cpl_parameter_get_int(p));


    cknull(distor_poly=sinfo_distortion_estimate_new(impoly,
                    0,
                    0,
                    cpl_image_get_size_x(impoly),
                    cpl_image_get_size_y(impoly),
                    FALSE,
                    SINFO_ARC_SATURATION,
                    SINFO_ARC_MAX_WIDTH,
                    arcs_thres_factor,
                    arcs_min_arclen_factor,
                    arcs_window_size,
                    smooth_rad,
                    3,
                    (double)value,
                    &arcs),
           "cannot estimage distortion") ;
    /*AMo: additional mem leaks */

    sinfo_free_apertures(&arcs);
    coef_pow[0]=0;
    coef_pow[1]=0;
    check_nomsg(pcf[0][0]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
    sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",0,0,pcf[0][0]);
    /*
  pcf[0][0]+=value;
     */
     sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",0,0,pcf[0][0]);
     check_nomsg(cpl_polynomial_set_coeff(distor_poly,coef_pow,pcf[0][0]));

     /* up to here ok */
     /* To check mem leaks */
     sinfo_free_image(&impoly);

     //sinfo_msg("Polynomial fit results: deg=%d",
     //      cpl_polynomial_get_degree(distor_poly));


     coef_pow[0]=0;
     coef_pow[1]=0;
     check_nomsg(pcf[0][0]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
     //sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",0,0,pcf[0][0]);

     coef_pow[0]=1;
     coef_pow[1]=0;
     check_nomsg(pcf[1][0]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
     //sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",1,0,pcf[1][0]);

     coef_pow[0]=0;
     coef_pow[1]=1;
     check_nomsg(pcf[0][1]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
     //sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",0,1,pcf[0][1]);

     coef_pow[0]=1;
     coef_pow[1]=1;
     check_nomsg(pcf[1][1]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
     //sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",1,1,pcf[1][1]);

     coef_pow[0]=2;
     coef_pow[1]=0;
     check_nomsg(pcf[2][0]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
     //sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",2,0,pcf[2][0]);

     coef_pow[0]=0;
     coef_pow[1]=2;
     check_nomsg(pcf[0][2]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
     //sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",0,2,pcf[0][2]);

     coef_pow[0]=2;
     coef_pow[1]=1;
     check_nomsg(pcf[2][1]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
     //sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",2,1,pcf[2][1]);

     coef_pow[0]=1;
     coef_pow[1]=2;
     check_nomsg(pcf[1][2]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
     //sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",1,2,pcf[1][2]);

     coef_pow[0]=3;
     coef_pow[1]=0;
     check_nomsg(pcf[3][0]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
     //sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",3,0,pcf[3][0]);

     coef_pow[0]=0;
     coef_pow[1]=3;
     check_nomsg(pcf[0][3]=cpl_polynomial_get_coeff(distor_poly,coef_pow));
     //sinfo_msg("Polynomial fit results: coeff[%d][%d]=%g",0,3,pcf[0][3]);

     sinfo_msg("Save results");

     check_nomsg(poly_tbl=cpl_table_new(10));
     check_nomsg(cpl_table_new_column(poly_tbl,"degx", CPL_TYPE_INT));
     check_nomsg(cpl_table_new_column(poly_tbl,"degy", CPL_TYPE_INT));
     check_nomsg(cpl_table_new_column(poly_tbl,"coeff", CPL_TYPE_DOUBLE));
     check_nomsg(cpl_table_set_int(poly_tbl,"degx",0,0));
     check_nomsg(cpl_table_set_int(poly_tbl,"degx",1,1));
     check_nomsg(cpl_table_set_int(poly_tbl,"degx",2,0));
     check_nomsg(cpl_table_set_int(poly_tbl,"degx",3,1));
     check_nomsg(cpl_table_set_int(poly_tbl,"degx",4,2));
     check_nomsg(cpl_table_set_int(poly_tbl,"degx",5,0));
     check_nomsg(cpl_table_set_int(poly_tbl,"degx",6,2));
     check_nomsg(cpl_table_set_int(poly_tbl,"degx",7,1));
     check_nomsg(cpl_table_set_int(poly_tbl,"degx",8,3));
     check_nomsg(cpl_table_set_int(poly_tbl,"degx",9,0));

     check_nomsg(cpl_table_set_int(poly_tbl,"degy",0,0));
     check_nomsg(cpl_table_set_int(poly_tbl,"degy",1,0));
     check_nomsg(cpl_table_set_int(poly_tbl,"degy",2,1));
     check_nomsg(cpl_table_set_int(poly_tbl,"degy",3,1));
     check_nomsg(cpl_table_set_int(poly_tbl,"degy",4,0));
     check_nomsg(cpl_table_set_int(poly_tbl,"degy",5,2));
     check_nomsg(cpl_table_set_int(poly_tbl,"degy",6,1));
     check_nomsg(cpl_table_set_int(poly_tbl,"degy",7,2));
     check_nomsg(cpl_table_set_int(poly_tbl,"degy",8,0));
     check_nomsg(cpl_table_set_int(poly_tbl,"degy",9,3));

     check_nomsg(cpl_table_set_double(poly_tbl,"coeff",0,pcf[0][0]));
     check_nomsg(cpl_table_set_double(poly_tbl,"coeff",1,pcf[1][0]));
     check_nomsg(cpl_table_set_double(poly_tbl,"coeff",2,pcf[0][1]));
     check_nomsg(cpl_table_set_double(poly_tbl,"coeff",3,pcf[1][1]));
     check_nomsg(cpl_table_set_double(poly_tbl,"coeff",4,pcf[2][0]));
     check_nomsg(cpl_table_set_double(poly_tbl,"coeff",5,pcf[0][2]));
     check_nomsg(cpl_table_set_double(poly_tbl,"coeff",6,pcf[2][1]));
     check_nomsg(cpl_table_set_double(poly_tbl,"coeff",7,pcf[1][2]));
     check_nomsg(cpl_table_set_double(poly_tbl,"coeff",8,pcf[3][0]));
     check_nomsg(cpl_table_set_double(poly_tbl,"coeff",9,pcf[0][3]));

     //sinfo_msg("Polynomial distortion coefficients \n");
     //sinfo_msg("c= %g   %g  %g  %g  %g  %g  %g  %g %g %g\n",
     //       pcf[0][0],pcf[1][0],pcf[0][1],pcf[1][1],
     //           pcf[2][0],pcf[0][2],pcf[2][1],pcf[1][2],pcf[3][0],pcf[0][3]);

     /* QC LOG */
     check_nomsg(qclog_tbl = sinfo_qclog_init());
     snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",0,0);
     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[0][0],
                     "Polynomial distortion coefficient","%g"));

     snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",1,0);
     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[1][0],
                     "Polynomial distortion coefficient","%g"));

     snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",0,1);
     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[0][1],
                     "Polynomial distortion coefficient","%g"));

     snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",1,1);
     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[1][1],
                     "Polynomial distortion coefficient","%g"));

     snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",2,0);
     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[2][0],
                     "Polynomial distortion coefficient","%g"));

     snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",0,2);
     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[0][2],
                     "Polynomial distortion coefficient","%g"));

     snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",2,1);
     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[2][1],
                     "Polynomial distortion coefficient","%g"));

     snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",1,2);
     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[1][2],
                     "Polynomial distortion coefficient","%g"));

     snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",3,0);
     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[3][0],
                     "Polynomial distortion coefficient","%g"));

     snprintf(key_name,MAX_NAME_SIZE-1,"%s%d%d","QC COEFF",0,3);
     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,pcf[0][3],
                     "Polynomial distortion coefficient","%g"));

     snprintf(key_name,MAX_NAME_SIZE-1,"%s","QC OFFSET");
     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,value,
                     "Polynomial distortion coefficient","%g"));

     xshift=new_compute_shift(x_c,y_c,pcf[0][0],pcf[1][0],pcf[0][1],
                     pcf[1][1],pcf[2][0],pcf[0][2],
                     pcf[2][1],pcf[1][2],pcf[3][0],pcf[0][3]);

     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC XSHIFT CC",xshift,
                     "X shift in x_c,y_c","%g"));


     xshift=new_compute_shift(x_l,y_l,pcf[0][0],pcf[1][0],pcf[0][1],
                     pcf[1][1],pcf[2][0],pcf[0][2],
                     pcf[2][1],pcf[1][2],pcf[3][0],pcf[0][3]);

     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC XSHIFT LL",xshift,
                     "X shift in x_l,y_l","%g"));

     xshift=new_compute_shift(x_l,y_u,pcf[0][0],pcf[1][0],pcf[0][1],
                     pcf[1][1],pcf[2][0],pcf[0][2],
                     pcf[2][1],pcf[1][2],pcf[3][0],pcf[0][3]);

     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC XSHIFT UL",xshift,
                     "X shift in x_l,y_u","%g"));

     xshift=new_compute_shift(x_u,y_u,pcf[0][0],pcf[1][0],pcf[0][1],
                     pcf[1][1],pcf[2][0],pcf[0][2],
                     pcf[2][1],pcf[1][2],pcf[3][0],pcf[0][3]);

     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC XSHIFT UR",xshift,
                     "X shift in x_u,y_u","%g"));


     xshift=new_compute_shift(x_u,y_l,pcf[0][0],pcf[1][0],pcf[0][1],
                     pcf[1][1],pcf[2][0],pcf[0][2],
                     pcf[2][1],pcf[1][2],pcf[3][0],pcf[0][3]);

     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC XSHIFT LR",xshift,
                     "X shift in x_u,y_l","%g"));


     ck0(sinfo_pro_save_tbl(poly_tbl,set_fibre_ns,sof,cfg->outName,
                     PRO_DISTORTION,qclog_tbl,plugin_id,config),
         "cannot dump tbl %s", cfg->outName);

     sinfo_free_table(&poly_tbl);
     sinfo_free_table(&qclog_tbl);
     sinfo_free_polynomial(&distor_poly);
     sinfo_free_int(&degx);
     sinfo_free_int(&degy);
     sinfo_free_double(&coef);
     sinfo_free_image(&im);
     sinfo_free_frameset(&stk);
     sinfo_finddist_free (&cfg);

     return 0;

     cleanup:
     sinfo_free_table(&poly_tbl);
     sinfo_free_table(&qclog_tbl);
     sinfo_free_polynomial(&distor_poly);
     sinfo_free_int(&degx);
     sinfo_free_int(&degy);
     sinfo_free_double(&coef);
     sinfo_free_apertures(&arcs);



     /*if(wave != NULL) sinfo_new_destroy_array (& wave );*/
     /*if(intens != NULL) sinfo_new_destroy_array (& intens );*/
     if(first != NULL) sinfo_new_destroy_array(&first);
     sinfo_free_table(&tbl_spos);
     if(slit_pos != NULL) sinfo_new_destroy_2Dfloatarray (&slit_pos,32);
     if(distances != NULL) sinfo_new_destroy_array(&distances);
     if(par!=NULL) sinfo_new_destroy_fit_params(&par);
     if(n_found_lines != NULL) sinfo_new_destroy_intarray(&n_found_lines );
     if(row_clean != NULL) sinfo_new_destroy_2Dintarray(&row_clean, lx);
     if(wavelength_clean != NULL) {
         sinfo_new_destroy_2Dfloatarray(&wavelength_clean,lx);
     }
     if(sum_pointer != NULL) sinfo_new_destroy_intarray(&sum_pointer);
     if(acoefs != NULL) {
         sinfo_new_destroy_2Dfloatarray(&acoefs, cfg->nrDispCoefficients );
     }
     sinfo_free_table(&tbl_line_list);
     sinfo_free_image(&mask);
     sinfo_free_image(&impoly);
     sinfo_free_image(&imonind);
     sinfo_free_image(&map);
     sinfo_free_image(&im);
     sinfo_finddist_free (&cfg);
     sinfo_free_frameset(&stk);
     return -1;

}


static double
new_compute_shift(double x,
                  double y,
                  double c00,
                  double c10,
                  double c01,
                  double c11,
                  double c20,
                  double c02,
                  double c21,
                  double c12,
                  double c30,
                  double c03)
{

    double x_shift=0;
    double shift=0;
    //sinfo_msg("c= %g   %g  %g  %g  %g  %g  %g  %g %g %g\n",
    //       c00,c10,c01,c11,
    //       c20,c02,c21,c12,c30,c03);

    shift=c00+c10*x+c01*y+
                    c11*x*y+c20*x*x+c02*y*y+
                    c21*x*x*y+c12*x*y*y+c30*x*x*x+c03*y*y*y;
    x_shift=x-shift;
    return x_shift;



}
/**@}*/