/* @(#)flames_fitting     */
/*===========================================================================
  Copyright (C) 2001 European Southern Observatory (ESO)

  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, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.

  Corresponding concerning ESO-MIDAS should be addressed as follows:
  Internet e-mail: midas@eso.org
  Postal address: European Southern Observatory
  Data Management Division 
  Karl-Schwarzschild-Strasse 2
  D 85748 Garching bei Muenchen 
  GERMANY
  ===========================================================================*/
/*-------------------------------------------------------------------------*/
/**
 * @defgroup flames_fitting  Follow orders. Fit the x,y pairs as obtained from flames_tracing 
 *
 */
/*-------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------
  Includes
  --------------------------------------------------------------------------*/
/**@{*/



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

#include <flames_fitting.h>
#include <flames_striptblext.h>

#include <flames_getordpos.h>
#include <flames_midas_tblsys.h>
#include <flames_midas_tblerr.h>
#include <flames_midas_tbldef.h>
#include <flames_midas_macrogen.h>
#include <flames_midas_atype.h>
#include <flames_midas_def.h>
#include <flames_uves.h>
#include <flames_newmatrix.h>
#include <flames_mvfit.h>
#include <uves_msg.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>


static int 
Derivatives(fit_info *loc_fitpar)
{
    int32_t mvalue, xvalue, fvalue, count_row, count_col;
    double mpower=0;

    double xpower=0;


    /* Loop over m-values of the polynomial */
    /* The "zero" degree parts (interfibre distance) are managed as a
     separate variable */

    for (count_row = 1; count_row <= loc_fitpar->n_xy; count_row++) {
        count_col = 0;

        /* treat separately the case xvalue = 0: here mvalue runs on values >0 */
        double prevmpower = 1;
        for(mvalue = 1; mvalue <= loc_fitpar->mdegree; mvalue++) {
            count_col++;
            (loc_fitpar->deriv)[count_row][count_col] = mpower =
                            prevmpower*(double) loc_fitpar->seq_order[count_row][1];
            prevmpower = mpower;
        }

        /* now run the loop for xvalue>0 */
        double prevxpower = 1;
        for (xvalue=1 ; xvalue <= loc_fitpar->xdegree; xvalue++) {
            /* treat separately mvalue=0 */
            count_col++;
            (loc_fitpar->deriv)[count_row][count_col] = xpower =
                            prevxpower*(loc_fitpar->f_xvalue)[count_row];
            prevxpower = xpower;
            for (mvalue=1 ; mvalue <= loc_fitpar->mdegree; mvalue++) {
                count_col++;
                (loc_fitpar->deriv)[count_row][count_col] =
                                (loc_fitpar->deriv)[count_row][mvalue]*xpower;
            }
        }

        for (fvalue=1; fvalue<=loc_fitpar->nfibres; fvalue++) {
            count_col++;
            (loc_fitpar->deriv)[count_row][count_col] =
                            (fvalue == loc_fitpar->seq_order[count_row][2]) ? 1 : 0;
        }
    }
    return(0);
}

static int dcompare(const double *pos1, const double *pos2)
{
    if (*pos1 < *pos2) {
        return(-1);
    }
    else if (*pos1 > *pos2) {
        return(1);
    }
    else {
        return(0);
    }
}

static int 
dcompare_qsort(const void *pos1, const void *pos2)
{
    return dcompare((const double *)pos1, (const double *)pos2);
}


/*---------------------------------------------------------------------------
  Implementation
  ---------------------------------------------------------------------------*/
/**
   @name  flames_checksize()  
   @short  check input frame size descriptors
   @author I. Porceddu -  ITAL_FLAMES Consortium. Ported to CPL by A. Modigliani


   @param HALFIBREWIDTH
   @param MAXFIBRES
   @param LIMIT
   @param NUMBER
   @param FIBRENUMBERS
   @param FIBRESHIFTS
   @param OTAB
   @param DEFPOL Degrees of the x,m parameters
   @param IN_B
   @param REFSTART
   @param REFSTEP
   @param REFNPIX
   @param CHIPCHOICE


   @return success or failure code (and the table with fitted order traces)

   DRS Functions called:          
   none                                         

   Pseudocode:                                                             
   check several MIDAS descriptors                                     

   @note
 */


int 
flames_fitting(
                const double *HALFIBREWIDTH,
                const int *MAXFIBRES,
                const int *LIMIT,
                const int *NUMBER,
                const int *FIBRENUMBERS,
                const double *FIBRESHIFTS,
                const char *OTAB,
                const int *DEFPOL,
                const char *IN_B,
                const double *REFSTART,
                const double *REFSTEP,
                const int *REFNPIX,
                const char *CHIPCHOICE)
{
    double dbuf=0;
    int status=0, selstatus=0;
    char dummytab[CATREC_LEN+4];
    char inptab[CATREC_LEN+1];
    char outptab[CATREC_LEN+1];
    char message[CATREC_LEN+15];
    char coeffc[21];
    int refitted=0;
    int cliploop=0;
    int maxiters=10;




    int i=0;
    int j=0;
    double kthreshold=2;
    double accepted=.9;
    int defpol[2]={0,0};
    int upointer=0;
    int null=0;
    int null2=0;
    int order_col=0;
    int fibre_col=0;
    int ordfib_col=0;
    char orderfibre[20];
    int order=0;
    int minorder=0;
    int fibre=0;
    int xcol=0;
    //  int ycol=0;
    double ycol=0;
    int actvals=0;
    int x_col=0;
    int y_col=0;
    int nbcol=0;
    int nbrow=0;

    /* Table ids' of data table */
    int tid=0;
    int yfit_col;
    int delta_col;
    int seq=0;
    int dim_pcoeffd=0;
    int row=0;
    int col=0;
    int start_median=0;
    int fibreson=0;
    double *fibrepos=0;
    double mean_par=0;
    double topx=0;
    double lowx=0;
    double xrange=0;
    double xrangepower=0;
    double xmax=0;
    double xmin=0;
    double **pcoeffd=0;
    double yvalue=0;
    double *compu_y=0;
    double *delta_y=0;
    double start[2]={0,0};
    double step[2]={0,0};
    char form[CATREC_LEN+1];
    int unit=0;
    int maxfibres=0;
    int nflats=0;
    double halfibrewidth=0;
    int npix[2]={0,0};
    int orderlim[2]={0,0};
    int coeffi[7]={0,0,0,0,0,0,0};
    float coeffr[5]={0,0,0,0,0};
    char chipchoice=0;
    int *fibremask=0;
    int *fibrenumbers=0;
    char trouble=0;
    int ifibre=0;
    char output[200];
    double *realfibrepos=0;
    int firstfibre=0;
    int lastfibre=0;
    double previousfibrepos=0;
    int previousfibre=0;
    double *fibreshifts=0;
    double *ordershifts=0;
    int medshiftindex=0;
    double medshift=0;

    fit_info *fitpar=0, *fitparsel1=0, *fitparsel2=0, *fitparsel3=0;
    orderpos *orderstruct=0;

    orderstruct = cpl_calloc(1, sizeof(orderpos));
    fitpar = cpl_calloc(1, sizeof(fit_info));
    fitparsel1 = cpl_calloc(1, sizeof(fit_info));
    fitparsel2 = cpl_calloc(1, sizeof(fit_info));
    memset(dummytab, 0, CATREC_LEN+4);
    memset(inptab, 0, CATREC_LEN+1);
    memset(outptab, 0, CATREC_LEN+1);
    memset(message, 0, CATREC_LEN+15);
    memset(form, 0, CATREC_LEN+1);

    memset(output, 0, 200);

    /* === interface to MIDAS =============================================== */

    /*Let's initialize the MIDAS environment */
    SCSPRO("flames_fitting");

    /* keyword HALFIBREWIDTH stores the half fibre width */
    SCKRDD(HALFIBREWIDTH, 1, 1, &actvals, &halfibrewidth, &upointer, &null);
    if (halfibrewidth <= 0) {
        SCTPUT("Error: HALFIBREWIDTH must be > 0!");
        return flames_midas_fail();
    }

    /* keyword MAXFIBRES stores the maximum number of fibres */
    SCKRDI(MAXFIBRES, 1, 1, &actvals, &maxfibres, &upointer, &null);
    if (maxfibres <= 0) {
        SCTPUT("Error: MAXFIBRES must be > 0!");
        return flames_midas_fail();
    }

    /* keyword limit stores the number of frames + 2 */
    /* jmlarsen: Changed limit to uppercase */
    SCKRDI(LIMIT, 1, 1, &actvals, &nflats, &upointer, &null);
    nflats -= 2;

    /* Reading NUMBER keyword (no. of detected fibers) */
    SCKRDI(NUMBER, 1, 1,&actvals, &fibreson, &upointer,&null);

    fibrenumbers = ivector(0, (nflats*maxfibres)-1);
    /* keyword FIBRENUMBERS contains the assignments between traces and fibres */
    SCKRDI(FIBRENUMBERS, 1, nflats*maxfibres, &actvals,
           fibrenumbers, &upointer, &null);

    fibreshifts = dvector(0, (nflats*maxfibres)-1);
    /* keyword FIBRESHIFTS contains the y shifts measured between the orders
     in the primeval guess order table used to get here and the positions
     of fibres as measured here */
    SCKRDD(FIBRESHIFTS, 1, nflats*maxfibres, &actvals,
           fibreshifts, &upointer, &null);

    /* initialise the fibremask, first clear it... */
    fibremask = ivector(0,(int32_t)maxfibres-1);
    for (i=0; i<=maxfibres-1; i++) {
        fibremask[i] = FALSE;
    }
    /* ... then match fibres which were found with the appropriate numbers,
     making sure that no fibre is repeated in more than one frame */
    trouble = FALSE;
    for (i=0; i<=fibreson-1; i++) {
        ifibre = fibrenumbers[i]-1;
        if (fibremask[ifibre] == TRUE) {
            /* trouble here: this fibre was assigned already! */
            sprintf(output, "Error: fibre %d is marked lit in more than one frame!",
                            ifibre+1);
            trouble = TRUE;
        }
        else fibremask[ifibre] = TRUE;
    }
    if (trouble==TRUE) {
        /* complain, clean up and exit with an error */
        SCTPUT(output);
        SCTPUT("Each fibre to be found and subsequently processed");
        SCTPUT("must be present in one and one only fibre/order");
        SCTPUT("positioning frame. Either (some of) the FIBREMASK(s)");
        SCTPUT("was set incorrectly or a fibre is indeed present in more");
        SCTPUT("than one frame. Please correct this problem and try again");
        free_ivector(fibremask, 0, maxfibres-1);
        free_ivector(fibrenumbers, 0, (nflats*maxfibres)-1);
        free_dvector(fibreshifts, 0, (nflats*maxfibres)-1);
        free(orderstruct);
        free(fitpar);
        free(fitparsel1);
        free(fitparsel2);
        return MAREMMA;
    }

    /* The OTAB  keyword refers to the table containing the ORDER, FIBRE
     X and Y values as computed from flames_tracing module   */
    //  SCKGETC(OTAB,1,60,&actvals,dummytab);
    //  if (striptblext(dummytab, inptab) != NOERR) {
    //    SCTPUT("Error stripping extension from input order table file name");
    //    return flames_midas_fail();
    //  }
    //jmlarsen: just do like this, we should use full filenames
    SCKGETC(OTAB,1,60,&actvals,inptab);

    /* Let's open the order table file (output table) */
    TCTOPN(inptab,F_IO_MODE,&tid);
    uves_msg_debug("input table:%s dummy tab: %s",inptab,dummytab);
    /* We read the parameters fully describing the inptab table */
    TCIGET (tid, &nbcol, &nbrow);

    /* We allocate the memory for the array which will host the order sequence */
    fitpar->seq_order = imatrix(1,(int32_t)nbrow,1,3);
    fitparsel1->seq_order = imatrix(1,(int32_t)nbrow,1,3);
    fitparsel2->seq_order = imatrix(1,(int32_t)nbrow,1,3);

    /* The DEFPOL descriptor does contain the actual value for the x and
     m degree values */
    SCKRDI(DEFPOL,1,2,&actvals,defpol,&upointer,&null);

    /* Filling up the x- and m-degree variables with the actual values */
    fitpar->xdegree = defpol[0];
    fitpar->mdegree = defpol[1];
    fitparsel1->xdegree = defpol[0];
    fitparsel1->mdegree = defpol[1];
    fitparsel2->xdegree = defpol[0];
    fitparsel2->mdegree = defpol[1];
    orderstruct->mdegree = defpol[1];
    orderstruct->xdegree = defpol[0];

    /* The IN_B keyword refers to the table middummr, which contains
     ORDER, FIBRE, X and Y values as computed from flames_tracing module */
    SCKGETC(IN_B,1,60,&actvals,dummytab);
    if (striptblext(dummytab, outptab) != NOERR) {
        SCTPUT("Error stripping extension from input order table file name");
        return flames_midas_fail();
    }

    fitpar->nfibres = fibreson;
    fitparsel1->nfibres = fibreson;
    fitparsel2->nfibres = fibreson;

    /* We allocate the memory for the array which will host the derivatives */
    fitpar->deriv = dmatrix(1,(int32_t)nbrow,1,
                    (int32_t)((defpol[1]+1)*(defpol[0]+1)+fibreson-1));
    fitparsel1->deriv = dmatrix(1,(int32_t)nbrow,1,
                    (int32_t)((defpol[1]+1)*(defpol[0]+1)+
                                    fibreson-1));
    fitparsel2->deriv = dmatrix(1,(int32_t)nbrow,1,
                    (int32_t)((defpol[1]+1)*(defpol[0]+1)+
                                    fibreson-1));

    /* We allocate the memory for the array which will host coeffd data */
    pcoeffd = dmatrix(0,fitpar->mdegree,0,fitpar->xdegree);
    for (i=0; i<=(fitpar->mdegree+1)*(fitpar->xdegree+1)-1; i++)
        *(pcoeffd[0]+i) = 0;

    /* We allocate the memory for the array which will host fibers'
     shifts data */
    fibrepos = dvector(0,(int32_t)maxfibres-1);
    for (i=0; i<=fibreson-1; i++) fibrepos[i] = 0;

    /* We allocate the memory for the vector which will host the y data to
     be fitted */
    fitpar->f_yvalue = dvector(1,(int32_t)nbrow);
    fitparsel1->f_yvalue = dvector(1,(int32_t)nbrow);
    fitparsel2->f_yvalue = dvector(1,(int32_t)nbrow);

    /* We allocate the memory for the vector which will host the x data to
     be fitted */
    fitpar->f_xvalue = dvector(1,(int32_t)nbrow);
    fitparsel1->f_xvalue = dvector(1,(int32_t)nbrow);
    fitparsel2->f_xvalue = dvector(1,(int32_t)nbrow);

    /* We allocate the memory for the vector which will host the variance
     for y data */
    fitpar->f_sigma = dvector(1,(int32_t)nbrow);
    fitparsel1->f_sigma = dvector(1,(int32_t)nbrow);
    fitparsel2->f_sigma = dvector(1,(int32_t)nbrow);

    /* Further filling the structure fitpar with number of parameters to be
     fitted ....*/
    fitpar->n_par = ((defpol[1]+1)*(defpol[0]+1)+fibreson - 1);
    fitparsel1->n_par = ((defpol[1]+1)*(defpol[0]+1)+fibreson - 1);
    fitparsel2->n_par = ((defpol[1]+1)*(defpol[0]+1)+fibreson - 1);

    /* We allocate the memory for the vector which will host the fitted
     parameters */
    fitpar->par = dvector(1,fitpar->n_par);
    fitparsel1->par = dvector(1,fitpar->n_par);
    fitparsel2->par = dvector(1,fitpar->n_par);

    /* ..... and number of points (equations) to be used  */
    fitpar->n_xy = (int32_t) nbrow;
    fitparsel1->n_xy = 0;
    fitparsel2->n_xy = 0;

    for (i=1; i<=nbrow; i++) {
        for (j=1; j<=3; j++) {
            fitpar->seq_order[i][j] = 0;
            fitparsel1->seq_order[i][j] = 0;
            fitparsel2->seq_order[i][j] = 0;
        }
        for (j=1; j<=((defpol[1]+1)*(defpol[0]+1) + fibreson-1); j++){
            fitpar->deriv[i][j] = 0;
            fitparsel1->deriv[i][j] = 0;
            fitparsel2->deriv[i][j] = 0;
        }
        fitpar->f_yvalue[i] = 0;
        fitparsel1->f_yvalue[i] = 0;
        fitparsel2->f_yvalue[i] = 0;
        fitpar->f_xvalue[i] = 0;
        fitparsel1->f_xvalue[i] = 0;
        fitparsel2->f_xvalue[i] = 0;
        fitpar->f_sigma[i] = 0;
        fitparsel1->f_sigma[i] = 0;
        fitparsel2->f_sigma[i] = 0;
    }
    for (i=0; i<=fitpar->n_par; i++){
        fitpar->par[i] = 0;
        fitparsel1->par[i] = 0;
        fitparsel2->par[i] = 0;
    }

    /* We allocate the memory for the vector which will host the y data to
     be fitted */
    compu_y = dvector(0,(int32_t)nbrow-1);
    delta_y = dvector(0,(int32_t)nbrow-1);
    for (i=0; i<=nbrow-1; i++) {
        compu_y[i] = 0;
        delta_y[i] = 0;
    }

    /* Matching the name of the column with its number */
    TCCSER (tid, "ORDER", &order_col);
    TCCSER (tid, "X", &x_col);
    TCCSER (tid, "Y", &y_col);
    TCCSER (tid, "YFIT", &yfit_col);
    TCCSER (tid, "RESIDUAL", &delta_col);
    TCCSER (tid, "FIBRE", &fibre_col);
    TCCSER (tid, "ORDERFIB", &ordfib_col);

    null=1;
    for (row=1; (row<=nbrow) && (null==1); row++)
        TCERDI (tid, row, order_col, &order,  &null);
    if (null==0) {
        minorder=order;
        uves_msg_debug("order=%d min=%d\n",order,minorder);
        for (row++; row<=nbrow; row++) {
            TCERDI (tid, row, order_col, &order,  &null);
            if (minorder>order) minorder=order;
        }
        minorder--;
        uves_msg_debug("order=%d min=%d\n",order,minorder);
    }
    else minorder=0;
    SCDWRI(tid, "TAB_IN_OUT_OSHIFT", &minorder, 1, 1, &null);

    for (row=1; row<=nbrow; row++) {
        /* read all values in fitpar, only selected ones in fitparsel1 */
        TCSGET (tid, row, &selstatus);
        TCERDI (tid, row, fibre_col, &fibre,  &null);
        TCERDI (tid, row, order_col, &order,  &null2);
        if ((null==0) && (null2==0)) {
            order -= minorder;
            TCEWRI(tid, row, order_col, &order);
            memset(orderfibre, 0, 20);
            sprintf(orderfibre, "%d,%d",order,fibrenumbers[fibre-1]);
            TCEWRC (tid, row, ordfib_col, orderfibre);
        }
        //jmlarsen: Hard to understand how this could ever work in the MIDAS pipeline
        //   'Y' contains floating point values, therefore should not be truncated.
        //   To avoid this read as floating point value, not integer.
        TCERDI (tid, row, x_col, &xcol,  &null);
        //    TCERDI (tid, row, y_col, &ycol,  &null);
        TCERDD (tid, row, y_col, &ycol,  &null);

        fitpar->f_xvalue[row] = (double)xcol; /* x value */
        //    fitpar->f_yvalue[row] =  (double)ycol; /* y value */
        fitpar->f_yvalue[row] =  ycol;

        fitpar->seq_order[row][1] = order; /* associated order to x and y values */
        fitpar->seq_order[row][2] = fibre; /* associated fibre to x and y values */
        fitpar->seq_order[row][3] = row; /* associated row to x and y values */
        fitpar->f_sigma[row] = 1.;
        if (selstatus == TRUE) {
            fitparsel1->n_xy++;
            fitparsel1->f_xvalue[fitparsel1->n_xy] = fitpar->f_xvalue[row];
            /* x value */
            fitparsel1->f_yvalue[fitparsel1->n_xy] = fitpar->f_yvalue[row];
            /* y value */
            fitparsel1->seq_order[fitparsel1->n_xy][1] = order;
            /* associated order
	 to x and y values */
            fitparsel1->seq_order[fitparsel1->n_xy][2] = fibre;
            /* associated fibre
	 to x and y values */
            fitparsel1->seq_order[fitparsel1->n_xy][3] = row;
            /* associated row to x
	 and y values */
            fitparsel1->f_sigma[fitparsel1->n_xy] = 1.;
        }
    }
    /* make sure we found some selected rows */
    if (fitparsel1->n_xy == 0) {
        sprintf(message, "Error: no selected entries in the %s table", inptab);
        SCTPUT(message);
        /* deallocate allocated memory here (TODO), then return an error*/
        return(MAREMMA);
    }

    /* Normalize the x-s  vector        BEGIN  */
    xmax = xmin = fitparsel1->f_xvalue[1];
    for(seq=2; seq <= fitparsel1->n_xy; seq++) {
        if(fitparsel1->f_xvalue[seq] > xmax)  xmax = fitparsel1->f_xvalue[seq];
        if(fitparsel1->f_xvalue[seq] < xmin)  xmin = fitparsel1->f_xvalue[seq];
    }
    for (seq = 1; seq <= fitparsel1->n_xy; seq++) {
        (fitparsel1->f_xvalue)[seq] = (fitparsel1->f_xvalue)[seq]/ (xmax-xmin);
    }
    topx = xmax;
    lowx = xmin;
    /* Normalize the x-s  vector        END  */

    /* Derivatives are computed here */
    if ((status=Derivatives(fitparsel1)) != NOERR) return(status);

    /* Actual fitting is done here */
    if ((status=mvfit(fitparsel1)) != NOERR) return(status);

    /*                           */
    /* Parameters are now written into FLAMESTABLE table's descriptors.
     We have been also back scaling the x values by using xmax */

    /* Computing the median value from the fibres constant values */
    start_median = (defpol[1]+1)*(defpol[0]+1) ;

    xmax = xmin = fitparsel1->par[start_median];
    for (row = start_median; row <= start_median + fibreson-1; row++) {
        if(fitparsel1->par[row] > xmax)  xmax=fitparsel1->par[row];
        if(fitparsel1->par[row] < xmin)  xmin=fitparsel1->par[row];
    }
    mean_par = (xmax + xmin)/2;

    /* Filling the pcoeffd array */
    pcoeffd[0][0] = mean_par;
    seq=0;
    for (col=1; col <= fitparsel1->mdegree; col++) {
        seq++;
        pcoeffd[col][0] = fitparsel1->par[seq];
    }
    xrange = topx-lowx;
    xrangepower = 1;
    for (row=1; row <= fitparsel1->xdegree; row++) {
        xrangepower *= xrange;
        for (col=0; col <= fitparsel1->mdegree; col++) {
            seq++;
            /*pcoeffd[col][row] = fitpar->par[seq];*/
            pcoeffd[col][row] = fitparsel1->par[seq]/xrangepower;
        }
    }
    for(row = 0; row <= fibreson-1; row++) {
        seq++;
        fibrepos[row] = fitparsel1->par[seq] - pcoeffd[0][0];
        uves_msg_debug("Fibre %d position: %f",
                       row, fibrepos[row]);
    }
    orderstruct->orderpol = pcoeffd;

    /* compute interpolated positions for all x's (selected and unselected) */
    for ( row=1; row <= nbrow; row++) {
        get_ordpos(orderstruct,
                        (double)fitpar->seq_order[row][1],
                        (double)fitpar->f_xvalue[row],&yvalue);
        compu_y[row-1] = yvalue+fibrepos[fitpar->seq_order[row][2]-1];
        delta_y[row-1] = fitpar->f_yvalue[row]-compu_y[row-1];
    }

    /* Here starts the sigma-clipping loop */
    refitted = TRUE;
    for (cliploop=1; (cliploop<=maxiters)&&(refitted==TRUE); cliploop++) {
        int oldorder = fitparsel1->seq_order[1][1];
        int oldfibre = fitparsel1->seq_order[1][2];
        int oldfibreend2 = fitparsel2->n_xy = 0;
        int oldfibrestart1 = 1;
        selstatus = FALSE;
        /* run over the previously selected values */
        for (row=1; row<=fitparsel1->n_xy; row++) {
            if (fitparsel1->seq_order[row][1] != oldorder ||
                            fitparsel1->seq_order[row][2] != oldfibre) {
                /* the fibre/order changed, check if the old one is to be wholly
	   discarded */
                /* are there enough pixels left to save the fibre/order? */
                if (fitparsel2->n_xy - oldfibreend2 <
                                (int) (ceil (accepted*(double) (row-oldfibrestart1)))) {
                    /* no, not enough, discard the fibre/order */
                    for (i=oldfibrestart1; i<=row-1; i++) {
                        /* deselect the whole fibre/order in the table */
                        TCSPUT(tid, fitparsel1->seq_order[i][3], &selstatus);
                    }
                    /* take away the fibre/order from the new selection */
                    fitparsel2->n_xy = oldfibreend2;
                }
                else {
                    oldfibreend2 = fitparsel2->n_xy;
                }
                oldfibrestart1 = row;
                oldorder = fitparsel1->seq_order[row][1];
                oldfibre = fitparsel1->seq_order[row][2];
            }
            /* is this pixel good? */
            if (fabs(delta_y[fitparsel1->seq_order[row][3]]) < kthreshold) {
                /* copy this pixel and all associated data to fitparsel2 */
                fitparsel2->n_xy++;
                fitparsel2->f_yvalue[fitparsel2->n_xy] = fitparsel1->f_yvalue[row];
                fitparsel2->f_xvalue[fitparsel2->n_xy] = fitparsel1->f_xvalue[row];
                fitparsel2->f_sigma[fitparsel2->n_xy] = fitparsel1->f_sigma[row];
                for (i=1; i<=fitparsel1->n_par; i++) {
                    fitparsel2->deriv[fitparsel2->n_xy][i] = fitparsel1->deriv[row][i];
                }
                for (i=1; i<=3; i++) {
                    fitparsel2->seq_order[fitparsel2->n_xy][i] =
                                    fitparsel1->seq_order[row][i];
                }
            }
        }
        /* check if I have to discard last fibre/order */
        if (fitparsel2->n_xy - oldfibreend2 <
                        (int) (ceil (0.8 * (double) (fitparsel1->n_xy+1-oldfibrestart1)))) {
            /* no, not enough, discard the fibre/order */
            for (i=oldfibrestart1; i<=fitparsel1->n_xy; i++) {
                /* deselect the whole fibre/order in the table */
                TCSPUT(tid, fitparsel1->seq_order[i][3], &selstatus);
            }
            /* take away the fibre/order from the new selection */
            fitparsel2->n_xy = oldfibreend2;
        }

        /* ok, sigma clipping complete, do I need to redo the fit? */
        if (fitparsel2->n_xy < fitparsel1->n_xy) {
            /* yes, I do, swap fitparsel2 and fitparsel1 */
            fitparsel3 = fitparsel1;
            fitparsel1 = fitparsel2;
            fitparsel2 = fitparsel3;
            /* redo the fit */
            if ((status=mvfit(fitparsel1)) != NOERR) return(status);
            /* Computing the median value from the fibres constant values */
            xmax = xmin = fitparsel1->par[start_median];
            for (row = start_median; row <= start_median + fibreson-1; row++) {
                if(fitparsel1->par[row] > xmax)  xmax=fitparsel1->par[row];
                if(fitparsel1->par[row] < xmin)  xmin=fitparsel1->par[row];
            }
            mean_par = (xmax + xmin)/2;
            /* Filling the pcoeffd array */
            pcoeffd[0][0] = mean_par;
            seq=0;
            for (col=1; col <= fitparsel1->mdegree; col++) {
                seq++;
                pcoeffd[col][0] = fitparsel1->par[seq];
            }
            xrangepower=1;
            for ( row=1; row <= fitparsel1->xdegree; row++) {
                xrangepower *= xrange;
                for (col=0; col <= fitparsel1->mdegree; col++) {
                    seq++;
                    pcoeffd[col][row] = fitparsel1->par[seq]/xrangepower;
                }
            }
            for( row = 0; row <= fibreson-1; row++ ) {
                seq++;
                fibrepos[row] = fitparsel1->par[seq] - pcoeffd[0][0];
                uves_msg_debug("Clipping loop %d: Fibre %d position: %f",
                               cliploop, row, fibrepos[row]);

            }
            /* compute interpolated positions for all x's
	 (selected and unselected) */
            for ( row=1; row <= nbrow; row++) {
                get_ordpos(orderstruct,
                                (double)fitpar->seq_order[row][1],
                                (double)fitpar->f_xvalue[row],&yvalue);
                compu_y[row-1] = yvalue+fibrepos[fitpar->seq_order[row][2]-1];
                delta_y[row-1] = fitpar->f_yvalue[row]-compu_y[row-1];
            }
            /* mark refitted as TRUE, so that it will retry clipping at next run */
            refitted = TRUE;
        }
        else {
            /* mark refitted as FALSE, the fit converged */
            refitted = FALSE;
        }
    } /* for cliploop */

    /* find first and last order on the table */
    orderlim[0] = orderlim[1] = fitpar->seq_order[1][1];
    for (row=1; row <= nbrow; row++) {
        if (fitpar->seq_order[row][1] < orderlim[0]) {
            orderlim[0] = fitpar->seq_order[row][1];
        }
        if (fitpar->seq_order[row][1] > orderlim[1]) {
            orderlim[1] = fitpar->seq_order[row][1];
        }
    }

    //jmlarsen: not used
    //nbord = (int)(orderlim[1]-orderlim[0]+1);
    //SCKWRI(ECHORD, &nbord, 1, 1, &null);

    /* Writing table values for COEFFI and COEFFD descriptors */
    memset(coeffc, 0, 21);
    strncpy(coeffc, inptab, 16);
    strcat(coeffc,"MULT");
    SCDWRC(tid, "COEFFC", 1, coeffc, 1, 20, &null);

    coeffi[0] = (int)fitparsel1->n_xy;
    coeffi[1] = 2;
    coeffi[2] = (int)y_col;
    /* coeffi[3] and coeffi[4] are expected by MIDAS to contain the columns
     of the two independent variables of the fit; in this case the 
     independent variables are actually three, therefore we are slightly 
     cheating here, but it is irrelevant anyway */
    coeffi[3] = (int)x_col;
    coeffi[4] = (int)order_col;
    coeffi[5] = (int)fitpar->xdegree;
    coeffi[6] = (int)fitpar->mdegree;
    SCDWRI(tid, "COEFFI", coeffi, 1, 7, &null);

    coeffr[0] = (float)lowx;
    coeffr[1] = (float)topx;
    coeffr[2] = (float)orderlim[0];
    coeffr[3] = (float)orderlim[1];
    coeffr[4] = 0;
    SCDWRR(tid, "COEFFR", coeffr, 1, 5, &null);

    dim_pcoeffd = (fitpar->mdegree+1)*(fitpar->xdegree+1);
    SCDWRD(tid, "COEFFD", pcoeffd[0],1, dim_pcoeffd,&null);

    /* compute order shift from fibre shifts */
    ordershifts = dvector(0, (nflats*maxfibres)-1);
    for (i=0; i<=fibreson-1; i++) ordershifts[i] = fibreshifts[i]-fibrepos[i];
    /* now find the median of the order shifts */
    qsort(ordershifts, (size_t) fibreson, sizeof(double), dcompare_qsort);
    /* now that the pixels are sorted, pick the median */
    if (2*(fibreson/2) != fibreson) {
        /* fibreson is odd, just pick the right index */
        medshiftindex = (fibreson-1)/2;
        medshift = ordershifts[medshiftindex];
    }
    else {
        /*  fibreson is even, compute the average of the two
	median values */
        medshiftindex = fibreson/2;
        medshift = (ordershifts[medshiftindex-1]+ordershifts[medshiftindex])/2;
    }
    free_dvector(ordershifts, 0, (nflats*maxfibres)-1);

    /* write the value we found to the TAB_IN_OUT_YSHIFT descriptor */
    SCDWRD(tid, "TAB_IN_OUT_YSHIFT", &medshift, 1, 1, &null);

    /* put fibrepos in proper order, according to provided fibremasks */
    realfibrepos = dvector(0,maxfibres-1);
    for (i=0; i<=maxfibres-1; i++) realfibrepos[i] = 0;
    firstfibre = lastfibre = fibrenumbers[0]-1;
    for (i=0; i<=fibreson-1; i++) {
        ifibre = fibrenumbers[i]-1;
        realfibrepos[ifibre] = fibrepos[i];
        if (ifibre<firstfibre) firstfibre=ifibre;
        if (ifibre>lastfibre) lastfibre=ifibre;
    }
    /* now check that they are indeed ordered (they ought to) */
    previousfibrepos = realfibrepos[firstfibre];
    previousfibre = firstfibre;
    for (i=firstfibre+1; i<=lastfibre; i++) {
        if (fibremask[i] == TRUE) {
            if (realfibrepos[i]>previousfibrepos){
                previousfibrepos = realfibrepos[i];
                previousfibre = i;
            }
            else {

                sprintf(output, "Error: fibre %d was found at position %g, \n \
while fibre %d was found at position %g", previousfibre, previousfibrepos, 
i, realfibrepos[i]);
                SCTPUT(output);
                SCTPUT("To produce this outcome, either some FIBREMASK was incorrect");
                SCTPUT("or some major instrument setup change occurred (earthquake?)");
                SCTPUT("or there is a software problem. Check out the first two");
                SCTPUT("possibilities first.");
                SCTPUT("AMO: 220602 Skipped the MAREMMA section");

                /*
	  free_dvector(fitpar->f_xvalue,1,(int32_t)nbrow);
	  free_dvector(fitpar->f_yvalue,1,(int32_t)nbrow);
	  free_dvector(fitpar->f_sigma,1,(int32_t)nbrow);
	  free_dvector(fitpar->par, 1, fitpar->n_par);
	  free_imatrix(fitpar->seq_order,1,(int32_t)nbrow,1,3);
	  free_dmatrix(fitpar->deriv,1,(int32_t)nbrow,1,
	  (int32_t)((defpol[1]+1)*(defpol[0]+1) + fibreson-1));
	  free_dvector(fitparsel1->f_xvalue,1,(int32_t)nbrow);
	  free_dvector(fitparsel1->f_yvalue,1,(int32_t)nbrow);
	  free_dvector(fitparsel1->f_sigma,1,(int32_t)nbrow);
	  free_dvector(fitparsel1->par, 1, fitpar->n_par);
	  free_imatrix(fitparsel1->seq_order,1,(int32_t)nbrow,1,3);
	  free_dmatrix(fitparsel1->deriv,1,(int32_t)nbrow,1,
	  (int32_t)((defpol[1]+1)*(defpol[0]+1) + fibreson-1));
	  free_dvector(fitparsel2->f_xvalue,1,(int32_t)nbrow);
	  free_dvector(fitparsel2->f_yvalue,1,(int32_t)nbrow);
	  free_dvector(fitparsel2->f_sigma,1,(int32_t)nbrow);
	  free_dvector(fitparsel2->par, 1, fitpar->n_par);
	  free_imatrix(fitparsel2->seq_order,1,(int32_t)nbrow,1,3);
	  free_dmatrix(fitparsel2->deriv,1,(int32_t)nbrow,1,
	  (int32_t)((defpol[1]+1)*(defpol[0]+1) + fibreson-1));
	  free_dvector(fibrepos, 0, maxfibres-1);
	  free_dvector(realfibrepos, 0, maxfibres-1);
	  free_ivector(fibrenumbers, 0, (nflats*maxfibres)-1);
	  free_dvector(fibreshifts, 0, (nflats*maxfibres)-1);
	  free_dvector(compu_y, 0,(int32_t) nbrow-1);
	  free_dvector(delta_y, 0,(int32_t) nbrow-1);
	  free_dmatrix(orderstruct->orderpol, 0, fitpar->mdegree, 
	  0, fitpar->xdegree);
	  free(orderstruct);
	  free(fitpar);
	  free(fitparsel1);
	  free(fitparsel2);
	  return MAREMMA;
                 */
            }
        }
    }

    /* update fibre numbers in the table with their reordered numbers */
    for (row=1; row<=nbrow; row++) {
        TCERDI (tid, row, fibre_col, &fibre,  &null);
        if (null==0) {
            fibre = fibrenumbers[fibre-1];
            TCEWRI(tid, row, fibre_col, &fibre);
        }
    }

    /* Writing table values for FIBREPOS descriptor */
    SCDWRD(tid, "FIBREPOS", realfibrepos,1, maxfibres, &null);

    /* write the fibremask */
    SCDWRI(tid, "FIBREMASK", fibremask, 1, maxfibres, &null);
    free_ivector(fibremask, 0, maxfibres-1);

    /* some more useful information to be stored as descriptors */
    SCDWRI(tid, "MAXFIBRES", &maxfibres, 1, 1, &null);
    SCDWRI(tid, "FIBRESON", &fibreson, 1, 1, &null);
    SCDWRD(tid, "HALFIBREWIDTH", &halfibrewidth, 1, 1, &null);
    SCDWRI(tid, "ORDERLIM", orderlim, 1, 2, &null);
    SCDWRC(tid, "CORRECTED", 1, "f", 1, 1, &null);

    /* create some descriptors which will be initialised in later steps of
     the DRS, right now they get dummy values */
    dbuf=0;
    SCDWRD(tid, "YCORRECTION", &dbuf, 1, 1, &null);
    SCDWRD(tid, "GAUSSFIBRESIGMA", &dbuf, 1, 1, &null);
    SCDWRD(tid, "GAUSSHALFWIDTH", &dbuf, 1, 1, &null);
    for (i=0; i<=(maxfibres-1); i++) fibrepos[i]=0;
    SCDWRD(tid, "GAUSSSELFSHIFT", fibrepos, 1, maxfibres, &null);

    SCKRDD(REFSTART, 1, 2, &actvals, start, &unit, &null);
    SCDWRD(tid, "REFSTART", start, 1, 2, &null);
    SCKRDD(REFSTEP, 1, 2, &actvals, step, &unit, &null);
    SCDWRD(tid, "REFSTEP", step, 1, 2, &null);
    SCKRDI(REFNPIX, 1, 2, &actvals, npix, &unit, &null);
    SCDWRI(tid, "REFNPIX", npix, 1, 2, &null);
    SCKRDC(CHIPCHOICE, 1, 1, 1, &actvals, &chipchoice, &unit, &null);
    SCDWRC(tid, "CHIPCHOICE", 1, &chipchoice, 1, 1, &null);
    npix[0]=0;
    SCDWRI(tid, "NAXIS", npix, 1, 1, &null);

    /* Writing table values for traced y values */
    for (row =1; row <= nbrow; row++) {
        TCEWRD(tid, row, yfit_col, &compu_y[row-1]);
        TCEWRD(tid, row, delta_col, &delta_y[row-1]);
    }

    TCTCLO(tid);

    free_dvector(fitpar->f_xvalue,1,(int32_t)nbrow);
    free_dvector(fitpar->f_yvalue,1,(int32_t)nbrow);
    free_dvector(fitpar->f_sigma,1,(int32_t)nbrow);
    free_dvector(fitpar->par, 1, fitpar->n_par);
    free_imatrix(fitpar->seq_order,1,(int32_t)nbrow,1,3);
    free_dmatrix(fitpar->deriv,1,(int32_t)nbrow,1,
                 (int32_t)((defpol[1]+1)*(defpol[0]+1) + fibreson-1));
    free_dvector(fitparsel1->f_xvalue,1,(int32_t)nbrow);
    free_dvector(fitparsel1->f_yvalue,1,(int32_t)nbrow);
    free_dvector(fitparsel1->f_sigma,1,(int32_t)nbrow);
    free_dvector(fitparsel1->par, 1, fitpar->n_par);
    free_imatrix(fitparsel1->seq_order,1,(int32_t)nbrow,1,3);
    free_dmatrix(fitparsel1->deriv,1,(int32_t)nbrow,1,
                 (int32_t)((defpol[1]+1)*(defpol[0]+1) + fibreson-1));
    free_dvector(fitparsel2->f_xvalue,1,(int32_t)nbrow);
    free_dvector(fitparsel2->f_yvalue,1,(int32_t)nbrow);
    free_dvector(fitparsel2->f_sigma,1,(int32_t)nbrow);
    free_dvector(fitparsel2->par, 1, fitpar->n_par);
    free_imatrix(fitparsel2->seq_order,1,(int32_t)nbrow,1,3);
    free_dmatrix(fitparsel2->deriv,1,(int32_t)nbrow,1,
                 (int32_t)((defpol[1]+1)*(defpol[0]+1) + fibreson-1));
    free_dvector(fibrepos, 0,(int32_t)fibreson-1);
    free_dvector(realfibrepos, 0, maxfibres-1);
    free_ivector(fibrenumbers, 0, (nflats*maxfibres)-1);
    free_dvector(fibreshifts, 0, (nflats*maxfibres)-1);
    free_dvector(compu_y, 0,(int32_t) nbrow-1);
    free_dvector(delta_y, 0,(int32_t) nbrow-1);
    free_dmatrix(orderstruct->orderpol, 0,fitpar->mdegree,0,fitpar->xdegree);
    free(orderstruct);
    free(fitpar);
    free(fitparsel1);
    free(fitparsel2);
    return SCSEPI();
}

/**@}*/