/*                                                                            *
 *   This file is part of the ESO UVES Pipeline                               *
 *   Copyright (C) 2004,2005 European Southern Observatory                    *
 *                                                                            *
 *   This library 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     *
 *                                                                           */

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif

/*---------------------------------------------------------------------------*/
/**
 * @defgroup flames_utils  Utility functions
 */
/*---------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------
                            Includes
 ----------------------------------------------------------------------------*/
#include <flames_utils.h>
#include <flames_def_drs_par.h>
#include <flames_midas_def.h>

#include <uves_dfs.h>
#include <uves_utils_wrappers.h>
#include <uves_pfits.h>
#include <uves_error.h>
#include <uves_dump.h>

#include <fitsio.h>

#include <ctype.h>
#include <string.h>
/*-----------------------------------------------------------------------------
                            Defines
 ----------------------------------------------------------------------------*/


/**@{*/
/*-----------------------------------------------------------------------------
                            Functions prototypes
 ----------------------------------------------------------------------------*/

/*-----------------------------------------------------------------------------
                            Implementation
 ----------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
/**
  @brief get file base name
  @param filename 
  @return base filename
 **/
/*---------------------------------------------------------------------------*/
char * 
flames_get_basename(const char *filename)
{
    char *p ;
    p = strrchr (filename, '/');
    return p ? p + 1 : (char *) filename;
}


/*---------------------------------------------------------------------------*/
/**
  @brief get 'flat size' half width and its shift to and order table
  @param frm input master flat frame
  @param frm input table name

  This test exists just to verify that function works
 **/
/*---------------------------------------------------------------------------*/
cpl_error_code
flames_get_mff_hw_and_yshift(cpl_frame* frm, 
                             const int save_flat_size,
                             const char* name_tbl)
{

    uves_propertylist* h=NULL;
    cpl_type type;
    cpl_image* ima=NULL;
    cpl_image* sim=NULL;
    cpl_image* ima_fit=NULL;
    cpl_image* ima_tmp=NULL;

    cpl_table* tbl=NULL;
    cpl_table* cor=NULL;
    cpl_table* fit=NULL;
    cpl_table* hwt=NULL;
    cpl_table* flat_size=NULL;

    float* psim=NULL;

    float* pint=NULL;
    float* psmo=NULL;
    float* prat=NULL;
    float* py=NULL;

    float* pima=NULL;
    float* pfit=NULL;
    float* pcor=NULL;
    float* ps=NULL;
    int* psize=NULL;

    int unit=0;
    float rms=0;
    const char* name_ima=NULL;
    int *coeffi = NULL;
    double *coeffd = NULL;
    float avg=0;
    float max=0;

    int ima_id=0;

    int ord_min=0;
    int ord_max=0;
    int ord_cen=0;
    int ord=0;
    int index=0;

    int xc_pix=0;
    int yc_pix=0;
    int yc=0;

    int defpol[2];

    int length=0;
    int i=0;
    int j=0;

    int sx=0;
    int sy=0;

    double fct=0;
    double yc_trace=0;
    double flux=0;
    int y_box_size=30; //updated on the fly to 1.5*halfwidth
    int lly=0;
    int ury=0;
    int s=5;
    int size=0;
    int sum=0;
    int count=0;
    int d=0;
    int xmin=0;
    int xmax=0;
    int found=0;
    int k=0;
    cpl_size row=0;
    int status=0;
    float yshift=0;
    float halfwidth=0;
    int ymin=0;
    int ymax=0;
    int rmin=0;


    check_nomsg(name_ima=cpl_frame_get_filename(frm));
    check_nomsg(ima=cpl_image_load(name_ima,CPL_TYPE_FLOAT,0,0));
    check_nomsg(tbl=cpl_table_load(name_tbl,1,0));
    check_nomsg(ord_min=cpl_table_get_column_min(tbl,"ORDER"));
    check_nomsg(ord_max=cpl_table_get_column_max(tbl,"ORDER"));
    //uves_msg_warning("norders=%d",ord_max-ord_min+1);

    check(h=uves_propertylist_load(name_tbl,0),"Loading %s header",name_tbl);

    check(coeffi=uves_read_midas_array(h,"COEFFI",&length,&type,NULL),
          "Error reading COEFFI from %s", name_tbl);

    assure( type == CPL_TYPE_INT, CPL_ERROR_TYPE_MISMATCH,
            "Type of COEFFI is %s, int expected",
            uves_tostring_cpl_type(type) );


    check(coeffd=uves_read_midas_array(h,"COEFFD",&length,&type,NULL),
          "Error reading COEFFD from %s", name_tbl);

    assure( type == CPL_TYPE_DOUBLE, CPL_ERROR_TYPE_MISMATCH,
            "Type of COEFFD is %s, int expected",
            uves_tostring_cpl_type(type) );

    defpol[0] = coeffi[5];
    defpol[1] = coeffi[6];
    //uves_msg_warning("defpol1 = %d, defpol2 = %d", defpol[0], defpol[1]);

    check_nomsg(sx=cpl_image_get_size_x(ima));
    check_nomsg(sy=cpl_image_get_size_y(ima));

    xc_pix=sx/2;
    ord_cen=(ord_min+ord_max)/2;

    check_nomsg(hwt=cpl_table_new(sy));
    check_nomsg(cpl_table_new_column(hwt,"INT",CPL_TYPE_FLOAT));
    check_nomsg(cpl_table_fill_column_window_float(hwt,"INT",0,sy,0));
    check_nomsg(pint=cpl_table_get_data_float(hwt,"INT"));
    check_nomsg(pima=cpl_image_get_data_float(ima));

    //Smooth a bit (radii=5 pix) the central image cross order trace
    for(j=s;j<sy-s-1;j++){
        avg=0;
        for(i=-s;i<=s;i++) {
            avg +=pima[(j+i)*sx+xc_pix];
        }
        pint[j]=avg/(2*s+1);
        //uves_msg_warning("int=%f",pint[j]);
    }
    size=(int)(sy/(ord_max-ord_min+1)/2);
    uves_msg_warning("size=%d",size);

    //Get the approximate master flat intensty profile
    check_nomsg(cpl_table_duplicate_column(hwt,"SMO",hwt,"INT"));
    check_nomsg(psmo=cpl_table_get_data_float(hwt,"SMO"));
    for(j=size;j<sy-size-1;j++){
        max=-10;
        for(i=-size/2;i<size/2;i++) {
            if(pint[j+i]>max) max=pint[j+i];
        }
        psmo[j]=max;
    }

    //Normalize the master flat cross order profile
    check_nomsg(cpl_table_duplicate_column(hwt,"RAT",hwt,"INT"));
    check_nomsg(cpl_table_divide_columns(hwt,"RAT","SMO"));
    check_nomsg(prat=cpl_table_get_data_float(hwt,"RAT"));
    /*
  check_nomsg(cpl_table_save(hwt, NULL, NULL,
                             uves_sprintf("%s%s",name_ima,"_hwt.fits"), CPL_IO_DEFAULT));
     */
    //Remove coss profile intensity oscillations of the top 20%
    for(j=0;j<sy;j++){
        prat[j]=(prat[j]<0.6) ? 0:1;
    }

    //Get the average master flat 'flat size'
    sum=0;
    count=0;
    found=0;
    check_nomsg(flat_size=cpl_table_new(0));
    check_nomsg(cpl_table_new_column(flat_size,"size",CPL_TYPE_INT));
    check_nomsg(psize=cpl_table_get_data_int(flat_size,"size"));

    //uves_msg_warning("size=%d-%d",size,sy-size);
    for(j=size;j<sy-size;j++){
        if((found==0) && (prat[j]<prat[j+1])) {
            xmin=j;
            found=1;
        }
        if((found==1) && (prat[j]>prat[j+1])) {
            xmax=j;
            found=2;
        }
        if((found==2) && (xmax > xmin)) {
            d=xmax-xmin;

            //uves_msg_warning("d=%d",d);
            count++;
            cpl_table_set_size(flat_size,count);
            //psize[count-1];
            check_nomsg(cpl_table_set_int(flat_size,"size",count-1,d));
            sum+=d;
            found=0;
        }
    }
    /*
  check_nomsg(cpl_table_save(flat_size,NULL,NULL,
  uves_sprintf("%s%s",name_ima,"size.fits"),CPL_IO_DEFAULT));
     */
    check_nomsg(d=(int)cpl_table_get_column_median(flat_size,"size"));
    uves_free_table(&flat_size);

    //d=sum/count;

    //Subtracts 1 to be on safe side
    d-=save_flat_size;
    halfwidth=d/2.;
    y_box_size=(int)(halfwidth+5);
    uves_msg_warning("halfwidth=%f",halfwidth);
    uves_msg_warning("y_box_size=%d",y_box_size);
    /*
  check_nomsg(cpl_table_save(hwt, NULL, NULL,
  uves_sprintf("%s%s",name_ima,"_hwt.fits"), CPL_IO_DEFAULT));
     */

    //uves_msg_warning("defpol[0]=%d,defpol[1]=%d",defpol[0],defpol[1]);


    //build a synthetic flat of 'flat size' equal to the input one,
    //centred on the order trace
    check_nomsg(sim=cpl_image_new(sx,sy,CPL_TYPE_FLOAT));
    check_nomsg(psim=cpl_image_get_data_float(sim));
    for (ord=ord_min;ord<ord_max; ord++){

        for(xc_pix=1;xc_pix<=sx;xc_pix++) {
            yc_trace=0;
            index=0;
            // get Y position of trace at central X pix
            for(j=0;j<=defpol[1];j++){

                fct = pow(ord,j);
                for(i=0;i<=defpol[0];i++){

                    yc_trace += coeffd[index]*pow(xc_pix,i)*fct;
                    ++index;
                    //uves_msg_warning("index=%d yc_trace[%d][%d]=%g",index,j,i,yc_trace);
                }

            }
            yc_pix=(int)(yc_trace+0.5);
            lly=(yc_pix-y_box_size>1)  ? yc_pix-y_box_size: 1;
            lly=(lly<sy)  ? lly: sy;
            ury=(yc_pix+y_box_size<sy) ? yc_pix+y_box_size: sy;
            ury=(ury>1) ? ury: 1;
            //uves_msg_warning("xc_pix=%d lly=%d",xc_pix,lly);
            //uves_msg_warning("yc_pix=%d ury=%d",yc_pix,ury);
            check_nomsg(flux=cpl_image_get_max_window(ima,xc_pix,lly,xc_pix,ury));
            //uves_msg_warning("ok min=%d max=%d",yc_pix-d/2,yc_pix+d/2);
            ymin=((yc_pix-d/2)>0)  ? -d/2 : -yc_pix;
            ymax=((yc_pix+d/2)<sy) ? d/2 : sy-yc_pix;
            //uves_msg("min=%d max=%d",ymin,ymax);
            for(yc=ymin;yc<ymax;yc++) {
                psim[xc_pix+(yc_pix+yc)*sx]=flux;
            }
            //uves_msg_warning("ord=%d xc_pix=%d yc_trace=%f f=%f",
            //	     ord,xc_pix,yc_trace,flux);

        }
    }

    uves_msg_warning("xc_pix=%dm yc_pix=%d",xc_pix,yc_pix);
    check_nomsg(ima_tmp=uves_image_smooth_x(sim,50));
    uves_free_image(&sim);
    check_nomsg(sim=cpl_image_duplicate(ima_tmp));
    uves_free_image(&ima_tmp);
    /*
  cpl_image_save(sim,uves_sprintf("%s%s",name_ima,"_sim.fits"), 
  CPL_BPP_IEEE_FLOAT,NULL,CPL_IO_DEFAULT);
     */
    //cross correlate the synthetic flat with the input flat
    check_nomsg(cor=cpl_table_new(4*size));
    cpl_table_new_column(cor,"Y",CPL_TYPE_FLOAT);
    cpl_table_new_column(cor,"COR",CPL_TYPE_FLOAT);

    check_nomsg(cpl_table_fill_column_window_float(cor,"COR",0,4*size,0));
    check_nomsg(cpl_table_fill_column_window_float(cor,"Y",0,4*size,0));
    check_nomsg(pcor=cpl_table_get_data_float(cor,"COR"));
    check_nomsg(ps=cpl_table_get_data_float(cor,"Y"));
    check_nomsg(psim=cpl_image_get_data_float(sim));
    check_nomsg(pima=cpl_image_get_data_float(ima));


    k=0;
    for(s=-2*size;s<2*size;s++){
        k++;
        flux=0;
        for(j=0;j<sy;j++) {
            for(i=0;i<sx;i++) {
                if((j-s>0) && (j-s)<sy) {
                    flux+=pima[i+j*sx]*psim[i+(j-s)*sx];
                }
            }
        }
        ps[k]=(float)s;
        pcor[k]=flux;
        //uves_msg_warning("s=%d,flux=%g",s,flux);

    }
    /*
  check_nomsg(cpl_table_save(cor, NULL, NULL,
  uves_sprintf("%s%s",name_ima,"_cor.fits"), CPL_IO_DEFAULT));
     */

    //Find the maximum of the cross correlation function
    check_nomsg(cpl_table_get_column_maxpos(cor,"COR",&row));
    check_nomsg(yshift=(float)cpl_table_get(cor,"Y",row,&status));
    rmin=( (row-d/2) >0)? (row-d/2):0;

    //uves_msg("row=%d row-d/4=%d d/2=%d",row,row-d/4,d/2);
    //check_nomsg(fit=cpl_table_extract(cor,row-d/4,d/2));
    check_nomsg(fit=cpl_table_extract(cor,rmin,d/2));
    check_nomsg(ima_fit=cpl_image_new(1,d/2,CPL_TYPE_FLOAT));
    check_nomsg(pcor=cpl_table_get_data_float(fit,"COR"));
    check_nomsg(py=cpl_table_get_data_float(fit,"Y"));
    check_nomsg(pfit=cpl_image_get_data_float(ima_fit));

    for(i=0;i<d/2;i++) {
        pfit[i]=pcor[i];
    }
    check_nomsg(yshift=(float)cpl_image_get_centroid_y_window(ima_fit,1,1,1,d/2));
    yshift+=py[0];
    //cpl_image_save(ima_fit,uves_sprintf("%s%s",name_ima,"_ima_fit.fits"),
    //CPL_BPP_IEEE_FLOAT,NULL,CPL_IO_DEFAULT);
    //uves_msg_warning("name_ima=%s",name_ima);
    SCFOPN(name_ima,D_R4_FORMAT,0,F_IMA_TYPE,&ima_id); // open input frame
    rms=0.5;
    SCDWRR(ima_id,"YSHIFT", &yshift,1,1,&unit);
    SCDWRR(ima_id,"ESO QC YSHIFT RMS", &rms,1,1,&unit);
    SCDWRR(ima_id,"HALFWIDTH", &halfwidth,1,1,&unit);
    SCDWRR(ima_id,"ESO QC HALFWIDTH RMS", &rms,1,1,&unit);
    SCFCLO(ima_id);


    cleanup:


    //uves_free_table(&tbl);
    uves_free_image(&ima);
    uves_free_image(&sim);
    uves_free_image(&ima_fit);
    uves_free_table(&flat_size);
    uves_free_table(&cor);
    uves_free_table(&fit);

    return cpl_error_get_code();


}


/**
@brief    
@param             input master flat frame
@param             input order table file name
@param             UVES chip id of frame
 */

cpl_error_code
msffsz_flames(cpl_frame *in_ima,
              const int save_flat_size,
              const char *in_tab,
              enum uves_chip chip)

{

    const char* name=NULL;
    uves_propertylist* header=NULL;
    cpl_image* img=NULL;

    check_nomsg(flames_get_mff_hw_and_yshift(in_ima, save_flat_size, in_tab) );
    check_nomsg(name=cpl_frame_get_filename(in_ima));

    check( header = uves_propertylist_load(name,0),
           "Could not load %s header", name);

    check_nomsg(flames_set_header_char_1(header,"CHIPCHOICE",
                    uves_chip_tochar(chip)));
    check_nomsg(img=cpl_image_load(name,CPL_TYPE_FLOAT,0,0));
    check_nomsg(uves_save_image(img,name,header,true,true));
    uves_free_image(&img);


    cleanup:
    uves_free_propertylist(&header);
    uves_free_image(&img);
    return cpl_error_get_code();


}


void
msffsz_flames2(cpl_frame *in_ima,
               const int save_flat_size,
               const char *in_tab,
               enum uves_chip chip)
{

    double slit_pix;
    int *coeffi = NULL;
    cpl_table *order_table = NULL;

    uves_propertylist *header = NULL;
    uves_propertylist *header_tab = NULL;
    int ord_min, ord_max;
    cpl_image* img=NULL;


    check( header = uves_propertylist_load(cpl_frame_get_filename(in_ima),
                    0),
           "Could not load %s header", cpl_frame_get_filename(in_ima));

    check_nomsg( slit_pix = 
                    uves_pfits_get_slitlength_pixels(header,
                                    chip) - 2.0);

    check( order_table = cpl_table_load(in_tab,
                    1, 1),
           "Failed to load order table %s", in_tab);

    ord_min = cpl_table_get_column_min(order_table, "ORDER");
    ord_max = cpl_table_get_column_max(order_table, "ORDER");


    {
        const char *IN_A = cpl_frame_get_filename(in_ima);
        const char *IN_B = in_tab;

        float slit = slit_pix;


        int defpol[2];
        cpl_type coeffi_type;
        int coeffi_length;

        check( header_tab = uves_propertylist_load(in_tab,
                        0),
               "Could not load %s header", in_tab);

        check( coeffi = uves_read_midas_array(
                        header_tab, "COEFFI", &coeffi_length,
                        &coeffi_type, NULL),
               "Error reading COEFFI from %s", in_tab);

        assure( coeffi_type == CPL_TYPE_INT, CPL_ERROR_TYPE_MISMATCH,
                "Type of COEFFI is %s, int expected",
                uves_tostring_cpl_type(coeffi_type) );

        defpol[0] = coeffi[5];
        defpol[1] = coeffi[6];

        uves_msg_debug("defpol1 = %d, defpol2 = %d", defpol[0], defpol[1]);

        uves_msg_debug("image = %s; table = %s", IN_A, IN_B);

        check_nomsg(flames_set_header_char_1(header,"CHIPCHOICE",
                        uves_chip_tochar(chip)));
        check_nomsg(img=cpl_image_load(IN_A,CPL_TYPE_FLOAT,0,0));
        check_nomsg(uves_save_image(img,IN_A,header,true,true));
        uves_free_image(&img);

        assure( 0 == flames_get_flat_size(IN_A,
                        IN_B,
                        &slit,
                        &ord_min,
                        &ord_max,
                        defpol,
                        save_flat_size,
                        FLAMES_DRS_SFF_HW_MIN,
                        FLAMES_X_WIND_SIZE,
                        FLAMES_Y_WIND_SIZE,
                        FLAMES_Y_SEARCH_WIND,
                        FLAMES_ORD_TRESH,
                        FLAMES_N_CLIP_MED,
                        FLAMES_N_CLIP_AVG,
                        FLAMES_INT_TRESH),
                CPL_ERROR_ILLEGAL_OUTPUT,
                "flames_get_flat_size failed");



    }

    cleanup:
    uves_free_table(&order_table);
    uves_free_int(&coeffi);
    uves_free_propertylist(&header);
    uves_free_propertylist(&header_tab);
    uves_free_image(&img);
    return;
}




/*---------------------------------------------------------------------------*/
/**
   @brief    Get FLAMES library binary version number
   @return   Binary version number

   This functions gives knowledge of the compile time version number of
   the FLAMES library at runtime.
   UVES_BINARY_VERSION is the same for the entire package.
 */
/*---------------------------------------------------------------------------*/
int
flames_get_version_binary(void)
{
    float version;
    float version_returnvalue = fits_get_version(&version);

    uves_msg_debug("Linking against CFITSIO version %g (%g)", version, 
                   version_returnvalue);
    return UVES_BINARY_VERSION;
}

/*---------------------------------------------------------------------------*/
/**
   @brief    Create frame from image buffer and header
   @param    filename          to save to
   @param    image             data to save
   @param    header            FITS header
   @return   newly allocated frame
 */
/*---------------------------------------------------------------------------*/
cpl_frame *
flames_new_table(const char *filename, 
                 const cpl_table *table,
                 const uves_propertylist *header)
{
    cpl_frame *f = NULL;

    f = cpl_frame_new();
    cpl_frame_set_filename(f, filename);
    cpl_frame_set_type(f, CPL_FRAME_TYPE_TABLE);

    check( uves_table_save(table, header, NULL, filename, CPL_IO_DEFAULT),
           "Error creating file %s from image", filename);

    cleanup:
    return f;
}

/*---------------------------------------------------------------------------*/
/**
   @brief    Create frame from image buffer and header
   @param    filename          to save to
   @param    image             data to save
   @param    header            FITS header
   @return   newly allocated frame
 */
/*---------------------------------------------------------------------------*/
cpl_frame *
flames_new_frame(const char *filename, 
                 const cpl_image *image,
                 const uves_propertylist *header)
{
    cpl_frame *f = NULL;

    f = cpl_frame_new();
    cpl_frame_set_filename(f, filename);
    cpl_frame_set_type(f, CPL_FRAME_TYPE_IMAGE);

    check( uves_save_image(image, filename, header,true, true),
           "Error creating file %s from image", filename);

    cleanup:
    return f;
}


/*---------------------------------------------------------------------------*/
/**
   @brief    Create frame from image buffer and header
   @param    filename          to save to
   @param    imagelist             data to save
   @param    header            FITS header
   @return   newly allocated frame
 */
/*---------------------------------------------------------------------------*/
cpl_frame *
flames_new_framelist(const char *filename, 
                     const cpl_imagelist *imagelist,
                     const uves_propertylist *header)
{
    cpl_frame *f = NULL;

    f = cpl_frame_new();
    cpl_frame_set_filename(f, filename);
    cpl_frame_set_type(f, CPL_FRAME_TYPE_IMAGE);

    check( uves_save_imagelist(imagelist, filename, header),
           "Error creating file %s from imagelist", filename);

    cleanup:
    return f;
}


/*---------------------------------------------------------------------------*/
/**
   @brief    Create frame from image and headers
   @param    filename          new filename or prefix
   @param    f                 frame to duplicate
   @param    prefix            whether the first parameter is the new filename
                               (false) or a prefix to attach (true)
   @param    reset_crval       if true, the frame is not just duplicated, but
                               CRVAL and CRPIX are set to one as in SPLIT/UVES

   @return   newly allocated frame
 */
/*---------------------------------------------------------------------------*/
cpl_frame *
flames_image_duplicate(const char *filename_prefix,
                       const cpl_frame *frame,
                       bool prefix,
                       bool reset_crval)
{
    cpl_frame *f = NULL;
    cpl_image *image = NULL;
    uves_propertylist *header = NULL;

    const char *filename = cpl_frame_get_filename(frame);
    const char *new_filename = prefix ?
                    uves_sprintf("%s%s", filename_prefix, filename) :
                    uves_sprintf("%s", filename_prefix);

    int plane = 0;
    int extension = 0;

    check( image = uves_load_image(frame, plane, extension, &header),
           "Could not load image");

    if (reset_crval) {
        check_nomsg( flames_reset_crval_to_one(&header));
    }

    f = cpl_frame_new();
    cpl_frame_set_filename(f, new_filename);
    cpl_frame_set_type(f, CPL_FRAME_TYPE_IMAGE);
    cpl_frame_set_tag(f, cpl_frame_get_tag(frame));

    check( uves_save_image(image, new_filename, header,true, true),
           "Error creating file %s from image", new_filename);

    cleanup:
    uves_free_image(&image);
    uves_free_propertylist(&header);
    uves_free_string_const(&new_filename);
    return f;
}

/*---------------------------------------------------------------------------*/
/**
   @brief    Get image median
   @param    f                 frame pointing to an image
   @return   CPL median
 */
/*---------------------------------------------------------------------------*/
double
flames_image_get_median(const cpl_frame *f)
{
    cpl_image *image = NULL;
    int plane = 0;
    int extension = 0;
    double result = 0;

    check( image = uves_load_image(f, plane, extension, NULL),
           "Could not load image");

    check( result = cpl_image_get_median(image),
           "Could not get median");

    cleanup:
    uves_free_image(&image);
    return result;
}

/*---------------------------------------------------------------------------*/
/**
   @brief    Create frame from table and headers
   @param    filename_prefix   string to prefix to first filename
   @param    f1                frame to be subtracted
   @param    f2                frame to subtract
   @return   newly allocated frame
 */
/*---------------------------------------------------------------------------*/
cpl_frame *
flames_image_subtract_create(const char *filename_prefix,
                             const cpl_frame *f1,
                             const cpl_frame *f2)
{
    cpl_frame *f = NULL;
    cpl_image *image1 = NULL;
    cpl_image *image2 = NULL;
    uves_propertylist *header1 = NULL;
    const char *filename = cpl_frame_get_filename(f1);
    const char *new_filename = uves_sprintf("%s%s", filename_prefix,
                    filename);

    int plane = 0;
    int extension = 0;

    check( image1 = uves_load_image(f1, plane, extension, &header1),
           "Could not load image");

    check( image2 = uves_load_image(f2, plane, extension, NULL),
           "Could not load image");

    check( cpl_image_subtract(image1, image2),
           "Error subtracting images");

    f = cpl_frame_new();
    cpl_frame_set_filename(f, new_filename);
    cpl_frame_set_type(f, CPL_FRAME_TYPE_IMAGE);

    check( uves_save_image(image1, new_filename, header1,true, true),
           "Error creating file %s from image", new_filename);

    cleanup:
    uves_free_image(&image1);
    uves_free_image(&image2);
    uves_free_propertylist(&header1);
    uves_free_string_const(&new_filename);
    return f;
}

/*---------------------------------------------------------------------------*/
/**
   @brief    Create frame from table and headers
   @param    filename_prefix   string to prefix to first filename
   @param    f1                frame to be subtracted
   @param    f2                frame to subtract
   @return   newly allocated frame
 */
/*---------------------------------------------------------------------------*/
cpl_frame *
flames_image_subtract_scalar_create(const char *filename_prefix,
                                    const cpl_frame *f,
                                    double value)
{
    cpl_frame *frame = NULL;
    cpl_image *image = NULL;
    uves_propertylist *header = NULL;
    const char *filename = cpl_frame_get_filename(f);
    const char *new_filename = uves_sprintf("%s%s", filename_prefix,
                    filename);

    int plane = 0;
    int extension = 0;

    check( image = uves_load_image(f, plane, extension, &header),
           "Could not load image");

    check( cpl_image_subtract_scalar(image, value),
           "Error subtracting images");

    frame = cpl_frame_new();
    cpl_frame_set_filename(frame, new_filename);
    cpl_frame_set_type(frame, CPL_FRAME_TYPE_IMAGE);

    check( uves_save_image(image, new_filename, header,true, true),
           "Error creating file %s from image", new_filename);

    cleanup:
    uves_free_image(&image);
    uves_free_propertylist(&header);
    uves_free_string_const(&new_filename);
    return frame;
}


/*---------------------------------------------------------------------------*/
/**
   @brief    Create frame from table and headers
   @param    filename          to save to
   @param    image             data to save
   @param    header            primary FITS header
   @param    theader            next FITS header
   @return   newly allocated frame
 */
/*---------------------------------------------------------------------------*/
cpl_frame *
flames_new_frame_table(const char *filename, 
                       const cpl_table *table,
                       const uves_propertylist *header,
                       const uves_propertylist *theader)
{
    cpl_frame *f = NULL;
    uves_propertylist *full_header;
    uves_propertylist *history;

    f = cpl_frame_new();
    cpl_frame_set_filename(f, filename);
    cpl_frame_set_type(f, CPL_FRAME_TYPE_TABLE);

    /* Copy all HISTORY keywords from extension header
       to primary header
     */
    check_nomsg( full_header = uves_propertylist_duplicate(header) );
    check_nomsg( history     = uves_propertylist_duplicate(theader) );
    check_nomsg( uves_propertylist_erase_regexp(history, "^HISTORY$", 1));
    check_nomsg( uves_propertylist_append(full_header, history));

    check( uves_table_save(table, full_header, full_header, filename, CPL_IO_DEFAULT),
           "Error creating file %s from table", filename);

    cleanup:
    uves_free_propertylist(&full_header);
    uves_free_propertylist(&history);
    return f;
}

/*---------------------------------------------------------------------------*/
/**
   @brief    Concatenate tables
   @param    n1       first table filename
   @param    n2       second table filename

   The combined table is stored in the first filename using the header from
   the first table
 */
/*---------------------------------------------------------------------------*/
void
flames_merge_table(const char *n1, const char *n2)
{
    cpl_table *t1 = NULL;
    cpl_table *t2 = NULL;
    uves_propertylist *t1_header = NULL;

    check( t1 = cpl_table_load(n1, 1, 1),
           "Could not load table %s", n1);

    check( t1_header = uves_propertylist_load(n1, 0),
           "Could not load table %s header", n1);

    check( t2 = cpl_table_load(n2, 1, 1),
           "Could not load table %s", n2);

    check_nomsg(cpl_table_insert(t1, t2, 0));

    check( uves_table_save(t1, t1_header, NULL, n1, CPL_IO_DEFAULT),
           "Could not save table to %s", n1);

    cleanup:
    uves_free_table(&t1);
    uves_free_table(&t2);
    uves_free_propertylist(&t1_header);
    return;
}

/*---------------------------------------------------------------------------*/
/**
   @brief    Rename table
   @param    n1       current name
   @param    n2       new name
 */
/*---------------------------------------------------------------------------*/
void
flames_rename_table(const char *n1, const char *n2)
{
    cpl_table *t = NULL;
    uves_propertylist *header = NULL;

    check( t = cpl_table_load(n1, 1, 1),
           "Could not load table %s", n1);
    check( header = uves_propertylist_load(n1, 0),
           "Could not load table %s header", n1);

    check( uves_table_save(t, header, NULL, n2, CPL_IO_DEFAULT),
           "Could not save table to %s", n2);
    cleanup:
    uves_free_table(&t);
    uves_free_propertylist(&header);
    return;           
}

/*---------------------------------------------------------------------------*/
/**
   @brief    Sort table after 2 columns
   @param    filename table filename
   @param    c1       first column
   @param    c2       second column
 */
/*---------------------------------------------------------------------------*/
void
flames_sort_table_2(const char *filename, const char *c1, const char *c2)
{
    cpl_table *t = NULL;
    uves_propertylist *header = NULL;

    check( t = cpl_table_load(filename, 1, 1),
           "Could not load table %s", filename);
    check( header = uves_propertylist_load(filename, 0),
           "Could not load table %s header", filename);

    /* Sort ascending */
    check( uves_sort_table_2(t, c1, c2,
                    false, false),
           "Sorting table %s failed", filename);

    check( uves_table_save(t, header, NULL, filename, CPL_IO_DEFAULT),
           "Could not save table to %s", filename);
    cleanup:
    uves_free_table(&t);
    uves_free_propertylist(&header);
    return;           
}

/*---------------------------------------------------------------------------*/
/**
   @brief    Select table rows
   @param    filename    table filename
   @param    column      a row is selected iff this column is non-NULL
 */
/*---------------------------------------------------------------------------*/
void
flames_select_non_null(const char *filename, const char *column)
{
    cpl_table *t = NULL;
    uves_propertylist *header = NULL;
    int i;

    check( t = cpl_table_load(filename, 1, 1),
           "Could not load table %s", filename);
    check( header = uves_propertylist_load(filename, 0),
           "Could not load table %s header", filename);

    assure_nomsg( cpl_table_has_column(t, "Select"), CPL_ERROR_DATA_NOT_FOUND );
    assure_nomsg( cpl_table_has_column(t, column), CPL_ERROR_DATA_NOT_FOUND);

    for (i = 0; i < cpl_table_get_nrow(t); i++)
    {
        if (cpl_table_is_valid(t, column, i)) {
            cpl_table_set_int(t, "Select", i, 1);
        }
        else {
            cpl_table_set_int(t, "Select", i, 0);
        }
    }

    check( uves_table_save(t, header, NULL, filename, CPL_IO_DEFAULT),
           "Could not save table to %s", filename);

    cleanup:
    uves_free_table(&t);
    uves_free_propertylist(&header);
    return;           
}

/*---------------------------------------------------------------------------*/
/**
   @brief    Select all table rows
   @param    filename    table filename
 */
/*---------------------------------------------------------------------------*/
void
flames_select_all(const char *filename)
{
    cpl_table *t = NULL;
    uves_propertylist *header = NULL;
    int i;

    check( t = cpl_table_load(filename, 1, 1),
           "Could not load table %s", filename);
    check( header = uves_propertylist_load(filename, 0),
           "Could not load table %s header", filename);

    assure_nomsg( cpl_table_has_column(t, "Select"), CPL_ERROR_DATA_NOT_FOUND );

    for (i = 0; i < cpl_table_get_nrow(t); i++)
    {
        cpl_table_set_int(t, "Select", i, 1);
    }

    check( uves_table_save(t, header, NULL, filename, CPL_IO_DEFAULT),
           "Could not save table to %s", filename);

    cleanup:
    uves_free_table(&t);
    uves_free_propertylist(&header);
    return;           
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Set crval keyword to 1
  @param    h    header to write to
 **/
/*---------------------------------------------------------------------------*/
void
flames_reset_crval_to_one(uves_propertylist** h)
{
    /* FLAMES MIDAS C-code wants the START (CRVAL) to be 1 */
    {
        double crval_old, shift;
        check_nomsg( crval_old = uves_pfits_get_crval1(*h) );
        uves_msg_debug("Old crval1=%f",crval_old);

        shift = 1 - crval_old;

        check_nomsg( uves_pfits_set_crval1(*h,uves_pfits_get_crval1(*h) + shift) );
        check_nomsg( uves_pfits_set_crpix1(*h,uves_pfits_get_crpix1(*h) + shift) );

        check_nomsg( crval_old = uves_pfits_get_crval2(*h) );
        uves_msg_debug("Old crval2=%f",crval_old);
        shift = 1 - crval_old;
        check_nomsg( uves_pfits_set_crval2(*h,uves_pfits_get_crval2(*h) + shift) );
        check_nomsg( uves_pfits_set_crpix2(*h,uves_pfits_get_crpix2(*h) + shift) );
    }

    /* Also force REFNPIX (CRPIX) to one */
    check_nomsg( uves_pfits_set_crpix1(*h, 1) );
    check_nomsg( uves_pfits_set_crpix2(*h, 1) );

    cleanup:
    return;
}

/*---------------------------------------------------------------------------*/
/**
   @brief    Sort table after 3 columns
   @param    filename table filename
   @param    c1       first column
   @param    c2       second column
   @param    c3       third column
 */
/*---------------------------------------------------------------------------*/
void
flames_sort_table(const char *filename, const char *c1, const char *c2, const char *c3)
{
    cpl_table *t = NULL;
    uves_propertylist *header = NULL;

    check( t = cpl_table_load(filename, 1, 1),
           "Could not load table %s", filename);
    check( header = uves_propertylist_load(filename, 0),
           "Could not load table %s header", filename);

    /* Sort ascending */
    check( uves_sort_table_3(t, c1, c2, c3,
                    false, false, false),
           "Sorting table %s failed", filename);

    check( uves_table_save(t, header, NULL, filename, CPL_IO_DEFAULT),
           "Could not save table to %s", filename);
    cleanup:
    uves_free_table(&t);
    uves_free_propertylist(&header);
    return;           
}


/*---------------------------------------------------------------------------*/
/**
   @brief    Identify the instrument setting
   @param    frame             reads from primary FITS header
   @param    chip              CCD chip
   @param    wc                (output) central wavelength (nm)
   @return   instrument conf identifier as a string, dynamically allocated

!.OUTPUT        string identifying: wcen, SimCal/OzPoz setting, 
!               Plate No., ccd id, binning

 */
/*---------------------------------------------------------------------------*/
const char *
flames_get_frmid(const cpl_frame *frame, enum uves_chip chip,
                 int *wc)
{

    const char *fid = NULL;        /* Result */

    const char *det;
    int bin1, bin2;
    const char *mod;
    const char *mod_dynamic = NULL;
    int plt;
    //    const char *modext;

    int wcenter;
    int idx;
    const char *dprtype;
    uves_propertylist *header = NULL;

    check( header = uves_propertylist_load(cpl_frame_get_filename(frame), 0),
           "Could not load header from %s", cpl_frame_get_filename(frame));

    assure_nomsg( wc != NULL, CPL_ERROR_NULL_INPUT );

    /* Not sure what q2 is, but 'det' ends being re or rm */
    det = uves_chip_get_det(chip);


    check( wcenter = uves_pfits_get_gratwlen(header, chip),
           "Error reading grating central wavelength");

    *wc = uves_round_double(wcenter);

    check( mod = uves_pfits_get_insmode(header),
           "Error reading instrument mode");

    check( plt = uves_flames_pfits_get_plateid(header),
           "Error reading plate ID");

    idx = strlen(mod) - 1;
    while(idx >= 0 && mod[idx] == ' ') idx--;   /* finds last non-blank character */

    assure( idx >= 0, CPL_ERROR_ILLEGAL_INPUT,
            "Illegal instrument mode string: '%s'", mod);

    if (isdigit(mod[idx]))
    {

        mod_dynamic = uves_sprintf("d%c", mod[idx]);
    }
    else
    {

        if (strstr(mod, "FIBRE") != NULL)
        {

            check( dprtype = uves_pfits_get_dpr_type(header),
                            "Error reading DPR TYPE");
            if (strstr(dprtype, "OzPoz") != NULL)
            {
                mod_dynamic = uves_sprintf("o%d", plt);                
            }
            else
            {
                mod_dynamic = uves_sprintf("s%d", plt);
            }
        }
        else
        {
            mod_dynamic = uves_sprintf("d0");
        }
    }

    //! The reversal of bin(2)xbin(1) serves for the
    //! naming convention before data rotation

    check( bin2 = uves_pfits_get_biny(header),
           "Error reading binning");
    check( bin1 = uves_pfits_get_binx(header),
           "Error reading binning");

    fid = uves_sprintf("%d%s%s%dx%d", *wc, mod_dynamic, det, bin2, bin1);

    cleanup:
    uves_free_propertylist(&header);
    uves_free_string_const(&mod_dynamic);
    return fid;
}
/**@}*/