/*===========================================================================
  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
  ===========================================================================*/
/* Program  : prepslitff.c                                                 */
/* Author   : G. Mulas  -  ITAL_FLAMES Consortium                          */
/* Date     :                                                              */
/*                                                                         */
/* Purpose  : Missing                                                      */
/*                                                                         */
/*                                                                         */
/* Input:  see interface                                                   */ 
/*                                                                      */
/* Output:                                                              */
/*                                                                         */
/* DRS Functions called:                                                   */
/* none                                                                    */ 
/*                                                                         */ 
/* Pseudocode:                                                             */
/* Missing                                                                 */ 
/*                                                                         */ 
/* Version  :                                                              */
/* Last modification date: 2002/08/05                                      */
/* Who     When        Why                Where                            */
/* AMo     02-08-05   Add header         header                            */
/*-------------------------------------------------------------------------*/

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

#include <flames_readslit0.h>
#include <flames_readslit.h>
#include <flames_prepslitff.h>
#include <flames_getordpos.h>
#include <flames_readordpos.h>
#include <flames_freeordpos.h>
#include <flames_stripfitsext.h>

#include <flames_midas_def.h>
#include <flames_uves.h>
#include <flames_newmatrix.h>
#include <flames_writeslitff.h>
#include <flames_freeslitflats.h>
#include <uves_msg.h>

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <ctype.h>

#define ULMAX(a,b) ((uint32_t)(a) > (uint32_t)(b) ? (uint32_t)(a) : (uint32_t)(b))
#define ULMIN(a,b) ((uint32_t)(a) < (uint32_t)(b) ? (uint32_t)(a) : (uint32_t)(b))

static int 
yshiftcompare(const void *slit1, const void *slit2);

static flames_err 
prepslitff(allslitflats *slitflats, orderpos *ordpos, 
           double decentSNR);




/**
 @brief prepares and normalise slit flat field data

 @param SLITCAT input slit FF data set
 @param OUTCAT  ouput slit FF data set
 @param MYORDER fibre order table
 @param BASENAME prefix to be used for slit FF data filenames
 @param DECENTSNR How large must the SNR on a fibre be in a calibration frame, at a given
                  order and x, for that slice to be considered "good"?

 @doc  -read Input frames and parameters
       -sort the FF frames in order of increasing YSHIFT
       -allocate and fill structures with data and finally calls actual prepslitff()
       -do slit FF cross-normalisation, then write output to disk
       -free memory
 */

int flames_prepslitff(const cpl_frameset *SLITCAT, cpl_frameset **OUTCAT,
                      const char *MYORDER, const char *BASENAME, const double *DECENTSNR) {
    int fileid = 0;
    int actvals = 0;
    int status = 0;
    int entrynum = 0;
    int unit = 0;
    int null = 0;
    //char incat[CATREC_LEN+1];
    const cpl_frameset *incat;
    //char outcat[CATREC_LEN+1];
    cpl_frameset **outcat;
    char basename[CATREC_LEN + 1];char
    filename[CATREC_LEN + 1];char
    ordername[CATREC_LEN + 1];char
    identifier[CATREC_LEN + 1];int32_t
    iframe = 0;
    double decentsnr = 0;

    orderpos *ordpos = 0;
    allslitflats *slitflats = 0;

    //memset(incat, 0, CATREC_LEN+1);
    //memset(outcat, 0, CATREC_LEN+1);
    memset(basename, 0, CATREC_LEN + 1);
    memset(filename, 0, CATREC_LEN + 1);
    memset(ordername, 0, CATREC_LEN + 1);
    memset(identifier, 0, CATREC_LEN + 1);

    /* allocate memory for the structures */
    ordpos = (orderpos *) calloc(1, sizeof(orderpos));
    slitflats = (allslitflats *) calloc(1, sizeof(allslitflats));

    /* enter the MIDAS environment */
    SCSPRO("prepslitff");

    /* read the SLITCAT keyword to know the name of the catalog file
   containing the list of int32_t slit FF frames */
    if ((status = SCKGETC_fs(SLITCAT, 1, 79, &actvals, &incat)) != 0) {
        /* the keyword seems undefined, protest... */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* get the outcat keyword */
    if ((status = SCKGETC_fsp(OUTCAT, 1, CATREC_LEN, &actvals, &outcat)) != 0) {
        /* I could not get the outcat keyword: complain... */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* get the myorder keyword */
    if ((status = SCKGETC(MYORDER, 1, CATREC_LEN, &actvals, ordername)) != 0) {
        /* I could not get the outcat keyword: complain... */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* get the basename for the files to be written to disk */
    if ((status = SCKGETC(BASENAME, 1, CATREC_LEN, &actvals, filename)) != 0) {
        /* I could not get the basename keyword: complain... */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }
    /* strip filename of the .fits extension, if it has one */
    if ((status = stripfitsext(filename, basename)) != NOERR) {
        /* error stripping extension */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* initialise DECENTSNR from keyword */
    if (SCKRDD(DECENTSNR, 1, 1, &actvals, &decentsnr, &unit, &null) != 0) {
        /* problems reading DECENTSNR */
        SCTPUT("Error reading the minimum acceptable SNR");
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* first open the input catalog and count the FF frames present */
    slitflats->nflats = 0;
    entrynum = 0;
    do {
        if ((status = SCCGET(incat,0,filename,identifier,&entrynum)) != 0) {
            /* error getting catalog entry */
            free(ordpos);
            free(slitflats);
            return flames_midas_fail();
        }
        /* did I get a valid catalog entry? */
        if (filename[0] != ' ') {
            slitflats->nflats++;
        }
    } while (filename[0] != ' ');

    /* check that nflats>0, otherwise, well... */
    if (slitflats->nflats == 0) {
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* in order to initialise firstorder and lastorder in slitflats, I
   need to read in the dummy table, therefore get that first */
    entrynum = 0;
    if ((status = SCCGET(incat,0,filename,identifier, &entrynum)) != 0) {
        /* error getting catalog entry */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* try to open the first frame */
    if ((status = SCFOPN(filename, FLAMESDATATYPE,
                    0, F_IMA_TYPE, &fileid)) != 0) {
        /* I could not open the frame */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    if ((status = SCDRDC(fileid, "CHIPCHOICE", 1, 1, 1, &actvals,
                    &slitflats->chipchoice, &unit, &null)) != 0) {
        /* something went wrong in SCDRDI */
        free(ordpos);
        free(slitflats);
        return (status);
    }

    /* close the first frame */
    if ((status = SCFCLO(fileid)) != 0) {
        /* problems closing the sigma frame */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* go for ordpos */
    /* initialise the ordpos structure from the dummy table descriptors */
    if ((status = readordpos(ordername, ordpos)) != NOERR) {
        /* something went wrong in the initialisation */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* check whether frames and order chip choices match */
    if (ordpos->chipchoice != slitflats->chipchoice) {
        /* no, they don't match */
        SCTPUT("Error: chip mismatch between frames and order table");
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* copy firstorder, lastorder and tab_io_oshift from ordpos */
    slitflats->firstorder = ordpos->firstorder;
    slitflats->lastorder = ordpos->lastorder;
    slitflats->tab_io_oshift = ordpos->tab_io_oshift;

    uves_msg_debug(
                    "Shifts = %d %d %d\n", ordpos->firstorder, ordpos->lastorder, ordpos->tab_io_oshift);

    /* initialise slitflats from the first frame */
    if ((status = readslit0(slitflats, 0, filename)) != NOERR) {
        /* error reading frame */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* now run a loop to fetch the other FF frames */
    for (iframe = 1; iframe <= slitflats->nflats - 1; iframe++) {
        entrynum = (int) iframe;
        /* get filename */
        if ((status = SCCGET(incat,0,filename,identifier,&entrynum)) != 0) {
            /* error getting catalog entry */
            return flames_midas_fail();
        }
        /* read frame */
        if ((status = readslit(slitflats, iframe, filename)) != NOERR) {
            /* error reading frame */
            return flames_midas_fail();
        }
    }
    /* good, all FF frames have been read in */
    /* take care to sort the FF frames in order of increasing YSHIFT */
    qsort(slitflats->slit, (size_t) slitflats->nflats, sizeof(slitFF),
          yshiftcompare);

    /* now that the structures have been allocated and filled with data,
   call the actual prepslitff */
    uves_msg("snr=%f", decentsnr);
    //writeslitff(slitflats, basename, outcat);

    if ((status = prepslitff(slitflats, ordpos, decentsnr)) != NOERR) {
        /* something went wrong with cross normalisation */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* the actual cross-normalisation has been done, here begins
   the phase of writing it down to disk */

    /* write it all to disk */
    if ((status = writeslitff(slitflats, basename, outcat)) != NOERR) {
        /* something went wrong writing this stuff to disk */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }

    /* free the memory allocated for the structures */
    if ((status = freeordpos(ordpos)) != NOERR) {
        /* error freeing ordpos */
        free(ordpos);
        free(slitflats);
        return flames_midas_fail();
    }
    free(ordpos);
    if ((status = freeslitflats(slitflats)) != NOERR) {
        /* error freeing slitflats */
        return flames_midas_fail();
    }
    free(slitflats);

    return (SCSEPI());

}

static flames_err prepslitff(allslitflats *slitflats, orderpos *ordpos,
                             double decentSNR) {

    flames_err status;
    int32_t iframe, iorder, ix, iy, goodpixels;
    double order, x, ordercentre, framecentre, plow, phigh;
    frame_data numerator, denominator, numerator2, denominator2;
    frame_data numsigma, densigma, pixel, sigma;
    int *framelist, *toclear;
    int32_t i, current, topmost, nframes, oldnframes;
    double currcentre, prevcentre;
    frame_data normfactor;
    slitFF *newslit;
    frame_data **dataframe, **sigmaframe;
    frame_mask **badframe;
    frame_data decenthreshold, decentSNR2;
    slitFF *myslit = 0;
    slitFF *currslit = 0;
    slitFF *topslit = 0;

    slitFF *firstslit = 0;


    frame_data fdbuf1 = 0;
    frame_data slitpixel = 0;
    frame_data slitpixel2 = 0;
    frame_data *fdvecbuf1 = 0;
    frame_data *fdvecbuf2 = 0;
    frame_data *fdvecbuf3 = 0;
    frame_data *fdvecbuf4 = 0;
    frame_data *fdvecbuf5 = 0;
    frame_data *fdvecbuf6 = 0;
    frame_data *fdvecbuf7 = 0;
    frame_mask *fmvecbuf1 = 0;
    frame_mask *fmvecbuf2 = 0;
    frame_mask *fmvecbuf3 = 0;
    frame_mask *fmvecbuf4 = 0;
    int32_t *lvecbuf1 = 0;
    int32_t *lvecbuf2 = 0;
    int32_t *lvecbuf3 = 0;
    int32_t *lvecbuf4 = 0;
    int32_t *lvecbuf5 = 0;
    int32_t *lvecbuf6 = 0;
    int32_t *lvecbuf7 = 0;
    int32_t *lvecbuf8 = 0;
    int32_t iyixuplimit = 0;
    int32_t iyixindex = 0;
    int32_t iorderixoffset = 0;
    int32_t iorderixindex = 0;
    int32_t commbottom = 0;
    int32_t commtop = 0;

    /* since it happens to have almost null slit values on unmasked bad
   columns, try to avoid them */
    decenthreshold = (frame_data) decentSNR
                    * pow(slitflats->gain * slitflats->ron, .5);
    uves_msg("thresh=%f", decenthreshold);
    decentSNR2 = (frame_data) (decentSNR * decentSNR);

    iyixuplimit = (slitflats->subrows * slitflats->subcols) - 1;

    lvecbuf1 = slitflats->lowbound[0];
    lvecbuf2 = slitflats->highbound[0];
    lvecbuf3 = slitflats->slit[0].lowbound[0];
    lvecbuf4 = slitflats->slit[0].highbound[0];
    fmvecbuf1 = slitflats->goodx[0];
    fdvecbuf1 = slitflats->normfactor[0];
    fdvecbuf2 = slitflats->slit[0].data[0];
    fdvecbuf3 = slitflats->slit[0].sigma[0];
    fmvecbuf2 = slitflats->slit[0].badpixel[0];

    /* before beginning the actual work, check that the flat part the orders
   in each frame is at least 1 pixel wide */
    for (iframe = 0; iframe <= slitflats->nflats - 1; iframe++) {
        if (2 * slitflats->slit[iframe].halfwidth / slitflats->substepy < 1) {
            return (MAREMMA);
        }
        myslit = slitflats->slit + iframe;
        fdvecbuf4 = myslit->data[0];
        fdvecbuf5 = myslit->sigma[0];
        fmvecbuf3 = myslit->badpixel[0];
        for (iyixindex = 0; iyixindex <= iyixuplimit; iyixindex++) {
            if ((fmvecbuf3[iyixindex] == 0)
                            && (((fdbuf1 = fdvecbuf4[iyixindex]) < decenthreshold)
                                            || ((fdbuf1 * fdbuf1 / fdvecbuf5[iyixindex]) < decentSNR2))) {
                fmvecbuf3[iyixindex] = 1;
                /*
         uves_msg("ok1 rat1=%f check1=%f rat2=%f check2=%f",
         fdbuf1,decenthreshold,
         fdbuf1*fdbuf1/fdvecbuf5[iyixindex],decentSNR2);
                 */
            }
        }
    }

    /* begin looping over orders at first */
    for (iorder = 0; iorder <= ((ordpos->lastorder) - (ordpos->firstorder));
                    iorder++) {
        iorderixoffset = iorder * slitflats->subcols;
        order = (double) (iorder + (ordpos->firstorder));
        /* loop over x... */
        for (ix = 0; ix <= (slitflats->subcols - 1); ix++) {
            iorderixindex = iorderixoffset + ix;
            /* convert the ix pixel coordinate to the x world coordinate */
            x = slitflats->substartx + (slitflats->substepx) * ((double) ix);
            /* find the unshifted central position and slope of this
       order at this x */
            /* bail out if the function call return an error status */
            if ((status = get_ordpos(ordpos, order, x, &ordercentre)) != NOERR) {
                return (status);
            }
            /* loop over FF frames */
            for (iframe = 0; iframe <= slitflats->nflats - 1; iframe++) {
                /* compute the order boundaries for each frame in this order in
         pixel coordinates*/
                myslit = slitflats->slit + iframe;
                lvecbuf5 = myslit->lowbound[0];
                lvecbuf6 = myslit->highbound[0];
                framecentre = ordercentre + myslit->yshift;
                /* remember that each pixel has a finite size, its position is
         the position of its centre and we want to be conservative in
         setting these boundaries */
                plow = ceil(
                                (framecentre - myslit->halfwidth) / slitflats->substepy + 0.5);
                phigh = floor(
                                (framecentre + myslit->halfwidth) / slitflats->substepy - 0.5);
                /* check boundaries and truncate to integer appropriately */
                /* is the whole interval out of the boundaries? */
                if ((phigh < 0) || (plow > (double) (slitflats->subrows) - 1)) {
                    /* the y loop must be skipped */
                    lvecbuf5[iorderixindex] = 1;
                    lvecbuf6[iorderixindex] = 0;
                } else {
                    /* is the upper limit above the upper boundary? */
                    if (phigh >= (double) (slitflats->subrows - 1)) {
                        /* yes it is */
                        lvecbuf6[iorderixindex] = slitflats->subrows - 1;
                    } else {
                        /* no it is not, truncate the upper limit to integer */
                        lvecbuf6[iorderixindex] = (int32_t) phigh;
                    }
                    /* is the lower limit below the lower boundary? */
                    if (plow <= 0) {
                        /* yes it is */
                        lvecbuf5[iorderixindex] = 0;
                    } else {
                        /* no it is not, truncate the lower limit to integer */
                        lvecbuf5[iorderixindex] = (int32_t) plow;
                    }
                }
            }
            /* set the goodx mask to good by default */
            fmvecbuf1[iorderixindex] = 0;
            /* set the overall boundaries equal to the boudaries for the first
       frame, at first */
            lvecbuf1[iorderixindex] = lvecbuf3[iorderixindex];
            lvecbuf2[iorderixindex] = lvecbuf4[iorderixindex];
            /* now compute the overall illumination factors; we do this using
       the first frame, since all subsequent frames will be scaled to
       match this */
            goodpixels = 0;
            numerator = 0;
            for (iy = lvecbuf1[iorderixindex]; iy <= lvecbuf2[iorderixindex]; iy++) {
                iyixindex = (iy * slitflats->subcols) + ix;
                /* take out insanely valued pixels before they cause real damage */
                if (fdvecbuf2[iyixindex] < decenthreshold) {
                    fmvecbuf2[iyixindex] = 1;
                    /*
           uves_msg("ok2 rat1=%f check1=%f",
           fdvecbuf2[iyixindex],decenthreshold);
                     */
                    //uves_msg("ok2");
                }
                /* is this pixel good? */
                if (fmvecbuf2[iyixindex] == 0) {
                    /* yes, it is, add its contribution to the overall
           illumination factor */
                    goodpixels++;
                    numerator += fdvecbuf2[iyixindex];
                }
            }
            /* did I find any good pixels? */
            if (goodpixels > 0
                            && (normfactor = numerator / ((frame_data) goodpixels))
                            > decenthreshold)
                fdvecbuf1[iorderixindex] = normfactor;
            else
                fdvecbuf1[iorderixindex] = 0;
            /* now compute the relative normalisation factors, for each frame
       relative to the previous one; of course this is skipped in case
       we have just one int32_t slit FF frame */
            myslit = slitflats->slit;
            lvecbuf5 = myslit->lowbound[0];
            lvecbuf6 = myslit->highbound[0];
            fdvecbuf4 = myslit->data[0];
            fdvecbuf5 = myslit->sigma[0];
            fmvecbuf3 = myslit->badpixel[0];
            for (iframe = 1; iframe <= slitflats->nflats - 1; iframe++) {
                /* compare the overall boundaries so far with the boundaries of
         this frame */
                lvecbuf7 = lvecbuf5;
                lvecbuf8 = lvecbuf6;
                fdvecbuf6 = fdvecbuf4;
                fdvecbuf7 = fdvecbuf5;
                fmvecbuf4 = fmvecbuf3;
                myslit = slitflats->slit + iframe;
                fdvecbuf4 = myslit->data[0];
                fdvecbuf5 = myslit->sigma[0];
                fmvecbuf3 = myslit->badpixel[0];
                lvecbuf5 = myslit->lowbound[0];
                lvecbuf6 = myslit->highbound[0];
                if (lvecbuf1[iorderixindex] > lvecbuf5[iorderixindex])
                    lvecbuf1[iorderixindex] = lvecbuf5[iorderixindex];
                if (lvecbuf2[iorderixindex] < lvecbuf6[iorderixindex])
                    lvecbuf2[iorderixindex] = lvecbuf6[iorderixindex];
                numerator = 0;
                numsigma = 0;
                denominator = 0;
                densigma = 0;
                goodpixels = 0;
                /* loop over the common part of the slit in iframe and the
         preceding one */
                if (lvecbuf5[iorderixindex] > lvecbuf7[iorderixindex])
                    commbottom = lvecbuf5[iorderixindex];
                else
                    commbottom = lvecbuf7[iorderixindex];
                if (lvecbuf6[iorderixindex] < lvecbuf8[iorderixindex])
                    commtop = lvecbuf6[iorderixindex];
                else
                    commtop = lvecbuf8[iorderixindex];
                for (iy = commbottom; iy <= commtop; iy++) {
                    iyixindex = (iy * slitflats->subcols) + ix;
                    /* take out insanely valued pixels before they cause real damage */
                    if (fdvecbuf4[iyixindex] < decenthreshold)
                        fmvecbuf3[iyixindex] = 1;
                    /*
           uves_msg("ok3 rat1=%f check1=%f",
           fdvecbuf3[iyixindex],decenthreshold);
                     */
                    //uves_msg("ok3");
                    /* is this pixel good in both frames? */
                    if ((fmvecbuf3[iyixindex] == 0) && (fmvecbuf4[iyixindex] == 0)) {
                        /* yes it is, add its contribution to all normalisation factors */
                        goodpixels++;
                        numerator += fdvecbuf6[iyixindex];
                        numsigma += fdvecbuf7[iyixindex];
                        denominator += fdvecbuf4[iyixindex];
                        densigma += fdvecbuf5[iyixindex];
                    }
                }
                /* were any overlapping good pixels found? */
                if ((goodpixels != 0) && (denominator > FDEPSILON)) {
                    /* yes, therefore do normalise */
                    for (iy = lvecbuf5[iorderixindex]; iy <= lvecbuf6[iorderixindex];
                                    iy++) {
                        iyixindex = (iy * slitflats->subcols) + ix;
                        /* first compute the variance of the normalised pixel */
                        numerator2 = numerator * numerator;
                        denominator2 = denominator * denominator;
                        slitpixel = fdvecbuf4[iyixindex];
                        slitpixel2 = slitpixel * slitpixel;
                        fdvecbuf5[iyixindex] = fdvecbuf5[iyixindex] * numerator2
                                        / denominator2 + numsigma * slitpixel2 / denominator2
                                        + densigma * slitpixel2 * numerator2
                                        / (denominator2 * denominator2);
                        /* now compute the scaled pixel value */
                        fdvecbuf4[iyixindex] = numerator * fdvecbuf4[iyixindex]
                                                                     / denominator;
                    }
                } else {
                    /* no overlapping good pixels found, mark this goodx as bad */
                    fmvecbuf1[iorderixindex] = 1;
                }
            }
            /* now run over the overall interval and compute each pixel as a
       weighted average of the available values for that position */
            if (fmvecbuf1[iorderixindex] == 0) {
                for (iy = lvecbuf1[iorderixindex]; iy <= lvecbuf2[iorderixindex];
                                iy++) {
                    iyixindex = (iy * slitflats->subcols) + ix;
                    denominator = 0;
                    numerator = 0;
                    goodpixels = 0;
                    /* cycle through frames to see which ones include this pixel */
                    for (iframe = 0; iframe <= slitflats->nflats - 1; iframe++) {
                        myslit = slitflats->slit + iframe;
                        fdvecbuf4 = myslit->data[0];
                        fdvecbuf5 = myslit->sigma[0];
                        fmvecbuf3 = myslit->badpixel[0];
                        lvecbuf5 = myslit->lowbound[0];
                        lvecbuf6 = myslit->highbound[0];
                        if ((iy >= lvecbuf5[iorderixindex])
                                        && (iy <= lvecbuf6[iorderixindex])) {
                            /* if it is a good pixel in this frame, add it to the average */
                            if (fmvecbuf3[iyixindex] == 0) {
                                goodpixels++;
                                numerator += fdvecbuf4[iyixindex] / fdvecbuf5[iyixindex];
                                denominator += 1 / fdvecbuf5[iyixindex];
                            }
                        }
                    }
                    /* if an average could be computed, put it back in the frames */
                    if (goodpixels != 0) {
                        /* average pixel value */
                        pixel = (frame_data) (numerator / denominator);
                        /* sigma of the average pixel value */
                        sigma = (frame_data) (1 / denominator);
                        /* cycle through frames to see which ones include this pixel */
                        for (iframe = 0; iframe <= slitflats->nflats - 1; iframe++) {
                            myslit = slitflats->slit + iframe;
                            fdvecbuf4 = myslit->data[0];
                            fdvecbuf5 = myslit->sigma[0];
                            fmvecbuf3 = myslit->badpixel[0];
                            lvecbuf5 = myslit->lowbound[0];
                            lvecbuf6 = myslit->highbound[0];
                            if ((iy >= lvecbuf5[iorderixindex])
                                            && (iy <= lvecbuf6[iorderixindex])) {
                                /* this is a good pixel now */
                                fmvecbuf3[iyixindex] = 0;
                                fdvecbuf4[iyixindex] = pixel;
                                fdvecbuf5[iyixindex] = sigma;
                            }
                        }
                    }
                }
            }
        }
    }

    /* find the minimum set of int32_t slit frames that covers the whole interval */
    framelist = ivector(0, slitflats->nflats - 1);
    /* find the lowest starting */
    current = 0;
    currslit = slitflats->slit;
    for (i = 1; i <= slitflats->nflats - 1; i++) {
        myslit = slitflats->slit + i;
        if ((myslit->yshift - myslit->halfwidth)
                        < (currslit->yshift - currslit->halfwidth)) {
            current = i;
            currslit = myslit;
        }
    }
    framelist[0] = topmost = current;
    topslit = currslit;
    nframes = 1;
    do {
        oldnframes = nframes;
        /* find the overlapping frame with the highest upper limit */
        for (i = 0; i <= slitflats->nflats - 1; i++) {
            myslit = slitflats->slit + i;
            /* is the lower limit of the i frame sufficiently lower than the
       upper limit of the current frame (i.e. can it overlap the current
       frame?)? */
            if ((topslit->yshift + topslit->halfwidth)
                            - (myslit->yshift - myslit->halfwidth) > slitflats->substepy) {
                /* is the upper limit of the i frame higher than the upper limit of
         the current frame? */
                if ((myslit->yshift + myslit->halfwidth)
                                > (currslit->yshift + currslit->halfwidth)) {
                    current = i;
                    currslit = myslit;
                }
            }
        }
        /* did I find another frame? */
        if (current != topmost) {
            /* add it to the list */
            framelist[nframes] = topmost = current;
            topslit = currslit;
            nframes++;
        }
    } while (nframes != oldnframes);

    /* ok, now I have the list of necessary frames, prepare a new
   slit array and put just the necessary frames in it */
    newslit = (slitFF *) calloc((size_t) nframes, sizeof(slitFF));
    toclear = ivector(0, slitflats->nflats - 1);
    for (i = 0; i <= slitflats->nflats - 1; i++) {
        toclear[i] = 0;
    }
    /* fill newslit with the necessary frames */
    for (i = 0; i <= nframes - 1; i++) {
        slitFF* mynewslit = newslit + i;
        myslit = slitflats->slit + framelist[i];
        mynewslit->data = myslit->data;
        mynewslit->sigma = myslit->sigma;
        mynewslit->badpixel = myslit->badpixel;
        mynewslit->yshift = myslit->yshift;
        mynewslit->framename = myslit->framename;
        mynewslit->sigmaname = myslit->sigmaname;
        mynewslit->badname = myslit->badname;
        mynewslit->boundname = myslit->boundname;
        mynewslit->halfwidth = myslit->halfwidth;
        mynewslit->lowbound = myslit->lowbound;
        mynewslit->highbound = myslit->highbound;
        toclear[framelist[i]] = 1;
    }
    /* free unused frames */
    for (i = 0; i <= slitflats->nflats - 1; i++) {
        /* is this slit unused? */
        if (toclear[i] == 0) {
            myslit = slitflats->slit + i;
            free_fdmatrix(myslit->data, 0, slitflats->subrows, 0, slitflats->subcols);
            free_fdmatrix(myslit->sigma, 0, slitflats->subrows, 0,
                          slitflats->subcols);
            free_fmmatrix(myslit->badpixel, 0, slitflats->subrows, 0,
                          slitflats->subcols);
            free_lmatrix(myslit->lowbound, 0, ordpos->lastorder - ordpos->firstorder,
                         0, slitflats->subcols);
            free_lmatrix(myslit->highbound, 0, ordpos->lastorder - ordpos->firstorder,
                         0, slitflats->subcols);
            free_cvector(myslit->framename, 0, CATREC_LEN);
            free_cvector(myslit->sigmaname, 0, CATREC_LEN);
            free_cvector(myslit->badname, 0, CATREC_LEN);
            free_cvector(myslit->boundname, 0, CATREC_LEN);
        }
    }
    /* free the old slit pointer */
    free(slitflats->slit);
    /* put the new pointer in slit */
    slitflats->slit = newslit;
    /* update the number of flats in the structure */
    slitflats->nflats = nframes;
    /* free temporary variables */
    free_ivector(framelist, 0, slitflats->nflats - 1);
    free_ivector(toclear, 0, slitflats->nflats - 1);

    /* to be even stricter, check whether the orders actually overlap;
   if not, I will simply put everything in one frame */
    /* of course, just do this if there is more than one frame left in the
   set of slit FF frames... */
    if (slitflats->nflats > 1) {
        /* run a loop over orders (from the second to the last) and x and find
     the absolute minimum of order separation */
        double fcurrmin = fabs((double) (slitflats->subrows + 1) * slitflats->substepy);
        for (iorder = 1; iorder <= ordpos->lastorder - ordpos->firstorder;
                        iorder++) {
            order = (double) (iorder + ordpos->firstorder);
            for (ix = 0; ix <= (slitflats->subcols - 1); ix++) {
                /* bail out if the function call return an error status */
                x = (double) (ix + slitflats->substartx);
                if ((status = get_ordpos(ordpos, order, x, &currcentre)) != NOERR) {
                    return (status);
                }
                double prevorder = (double) (iorder - 1 + ordpos->firstorder);
                if ((status = get_ordpos(ordpos, prevorder, x, &prevcentre)) != NOERR) {
                    return (status);
                }
                if (fabs(currcentre - prevcentre) < fcurrmin) {
                    fcurrmin = fabs(currcentre - prevcentre);
                }
            }
        }
        /* now I know the lowest order separation, is there overlap? */
        firstslit = slitflats->slit;
        slitFF* lastslit = firstslit + slitflats->nflats - 1;
        if (firstslit->yshift - firstslit->halfwidth + fcurrmin
                        - (lastslit->yshift + lastslit->halfwidth) > slitflats->substepy) {
            /* no overlap, therefore put all pixels in the first frame! */
            dataframe = fdmatrix(0, slitflats->subrows - 1, 0,
                            slitflats->subcols - 1);
            fdvecbuf4 = dataframe[0];
            memset(fdvecbuf4, 0,
                   slitflats->subrows * slitflats->subcols * sizeof(frame_data));
            sigmaframe = fdmatrix(0, slitflats->subrows - 1, 0,
                            slitflats->subcols - 1);
            fdvecbuf5 = sigmaframe[0];
            memset(fdvecbuf5, 0,
                   slitflats->subrows * slitflats->subcols * sizeof(frame_data));
            badframe = fmmatrix(0, slitflats->subrows - 1, 0, slitflats->subcols - 1);
            fmvecbuf3 = badframe[0];
            for (iyixindex = 0; iyixindex <= iyixuplimit; iyixindex++) {
                fmvecbuf3[iyixindex] = 1;
            }
            for (iorder = 0; iorder <= ordpos->lastorder - ordpos->firstorder;
                            iorder++) {
                iorderixoffset = iorder * slitflats->subcols;
                for (ix = 0; ix <= (slitflats->subcols - 1); ix++) {
                    iorderixindex = iorderixoffset + ix;
                    lvecbuf5 = lvecbuf3;
                    for (iframe = 0; iframe <= slitflats->nflats - 1; iframe++) {
                        myslit = slitflats->slit + iframe;
                        fdvecbuf6 = myslit->data[0];
                        fdvecbuf7 = myslit->sigma[0];
                        fmvecbuf4 = myslit->badpixel[0];
                        lvecbuf6 = myslit->highbound[0];
                        for (iy = lvecbuf5[iorderixindex]; iy <= lvecbuf6[iorderixindex];
                                        iy++) {
                            iyixindex = (iy * slitflats->subcols) + ix;
                            fdvecbuf4[iyixindex] = fdvecbuf6[iyixindex];
                            fdvecbuf5[iyixindex] = fdvecbuf7[iyixindex];
                            fmvecbuf3[iyixindex] = fmvecbuf4[iyixindex];
                        }
                        lvecbuf5 = lvecbuf6;
                    }
                    lvecbuf4[iorderixindex] = lvecbuf2[iorderixindex];
                }
            }
            /* now the ugly memory handling */
            /* copy the pointers to the new slit */
            newslit = (slitFF *) calloc(1, sizeof(slitFF));
            newslit->data = dataframe;
            newslit->sigma = sigmaframe;
            newslit->badpixel = badframe;
            newslit->framename = firstslit->framename;
            newslit->sigmaname = firstslit->sigmaname;
            newslit->badname = firstslit->badname;
            newslit->boundname = firstslit->boundname;
            newslit->yshift = (firstslit->yshift - firstslit->halfwidth
                            + lastslit->yshift + lastslit->halfwidth) / 2;
            newslit->halfwidth = (lastslit->yshift + lastslit->halfwidth
                            - firstslit->yshift + firstslit->halfwidth) / 2;
            newslit->lowbound = firstslit->lowbound;
            newslit->highbound = firstslit->highbound;
            /* free unneeded arrays */
            free_fdmatrix(firstslit->data, 0, slitflats->subrows, 0,
                          slitflats->subcols);
            free_fdmatrix(firstslit->sigma, 0, slitflats->subrows, 0,
                          slitflats->subcols);
            free_fmmatrix(firstslit->badpixel, 0, slitflats->subrows, 0,
                          slitflats->subcols);
            for (iframe = 1; iframe <= slitflats->nflats - 1; iframe++) {
                myslit = firstslit + iframe;
                free_fdmatrix(myslit->data, 0, slitflats->subrows, 0,
                              slitflats->subcols);
                free_fdmatrix(myslit->sigma, 0, slitflats->subrows, 0,
                              slitflats->subcols);
                free_fmmatrix(myslit->badpixel, 0, slitflats->subrows, 0,
                              slitflats->subcols);
                free_lmatrix(myslit->lowbound, 0,
                             ordpos->lastorder - ordpos->firstorder, 0, slitflats->subcols);
                free_lmatrix(myslit->highbound, 0,
                             ordpos->lastorder - ordpos->firstorder, 0, slitflats->subcols);
                free_cvector(myslit->framename, 0, CATREC_LEN);
                free_cvector(myslit->sigmaname, 0, CATREC_LEN);
                free_cvector(myslit->badname, 0, CATREC_LEN);
                free_cvector(myslit->boundname, 0, CATREC_LEN);
            }
            /* free the old slit pointer */
            free(slitflats->slit);
            /* replace it with the new pointer */
            slitflats->slit = newslit;
            /* set the number of flats in the structure to 1 */
            slitflats->nflats = 1;
        }
    }

    /* time to write to disk, isn't it? */

    return NOERR;

}

static int yshiftcompare(const void* s1, const void *s2) {
    const slitFF *slit1 = (const slitFF *) s1;
    const slitFF *slit2 = (const slitFF *) s2;

    if (slit1->yshift < slit2->yshift) {
        return (-1);
    } else if (slit1->yshift > slit2->yshift) {
        return (1);
    } else {
        return (0);
    }
}