/* @(#)create_ordertable.c */
/*===========================================================================
  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_backfit   Substep: fit the background 
 *
 */
/*-------------------------------------------------------------------------*/

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

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


#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <flames_freeordpos.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 <math.h>
#include <flames_midas_def.h>
#include <flames_uves.h>
#include <flames_getordpos.h>
#include <flames_newmatrix.h>

#include <uves_msg.h>
#include <flames_create_backtable.h>
#include <flames_readordpos.h>

/**@{*/

/*---------------------------------------------------------------------------
  Implementation
  ---------------------------------------------------------------------------*/
/**
   @name  flames_create_backtable()  
   @short Create background table                                             
   @author G. Mulas  -  ITAL_FLAMES Consortium. Ported to CPL by A. Modigliani

   @param MYORDER         User-defined name of the input order table          
   @param MYBACK        User-defined name of the created background table  
   @param WINSPERORDER    maximum number of background windows per interorder 
   @param MYWIDTHS          x and y maximum widths of each window             

   @return initialized frame

   DRS Functions called:                                                
   none                                                                 

   Pseudocode:                                                          
   Missing                                                               


   @note
 */


int 
flames_create_backtable(
                const char *MYORDER,
                const char *MYBACK,
                const int *WINSPERORDER,
                const double *MYWIDTHS)
{
    char ordtab[CATREC_LEN+4];
    char backtab[CATREC_LEN+4];
    char output[160];
    int winsperorder=0;
    double winstep=0;
    double xminwidth=0;
    double halfxminwidth=0;
    double firstx=0;

    double *xcentres=0;
    double *xstarts=0;
    double *xends=0;
    double yminwidth=0;
    double halfyminwidth=0;
    double ymax=0;
    double ycentre=0;
    double ystart=0;
    double yend=0;
    double yswap=0;
    double dorder=0;
    int32_t ixcount=0;
    int32_t iorder=0;
    int tid=0;
    int ordcol=0;
    int xcol=0;
    int xstartcol=0;
    int xendcol=0;
    int ycol=0;
    int ystartcol=0;
    int yendcol=0;
    int windowcounter=0;
    int unit=0;
    int null=0;
    int actvals=0;
    orderpos *ordpos=0;

    memset(ordtab, 0, CATREC_LEN+4);
    memset(backtab, 0, CATREC_LEN+4);
    memset(output, 0, 160);

    SCSPRO("create_backtable");

    /* read procedure input parameters */
    if (SCKGETC(MYORDER, 1, 160, &actvals, ordtab) != 0) {
        SCTPUT("Error reading MYORDER keyword");
        return flames_midas_fail();
    }

    if (SCKGETC(MYBACK, 1, 160, &actvals, backtab) != 0) {
        SCTPUT("Error reading MYORDER keyword");
        return flames_midas_fail();
    }

    if (SCKRDI(WINSPERORDER, 1, 1, &actvals, &winsperorder, &unit, &null)!=0){
        SCTPUT("Error reading WINSPERORDER keyword");
        return flames_midas_fail();
    }

    if (SCKRDD(MYWIDTHS, 1, 1, &actvals, &xminwidth, &unit, &null)!=0 ||
                    SCKRDD(MYWIDTHS, 2, 1, &actvals, &yminwidth, &unit, &null)!=0) {
        SCTPUT("Error reading MYWIDTHS keyword");
        return flames_midas_fail();
    }


    /* allocate order table */
    if ((ordpos=calloc(1, sizeof(orderpos))) == NULL) {
        SCTPUT("Error allocating ordpos: out of memory?");
        return flames_midas_fail();
    }

    /* read order table */
    if (readordpos(ordtab, ordpos) != NOERR) {
        sprintf(output, "Error reading %s order table", ordtab);
        SCTPUT(output);
        return flames_midas_fail();
    }


    /* some consistency/sanity checks on the input parameters */
    if (winsperorder*2>ordpos->npix[0]) {
        winsperorder = ordpos->npix[0]/2;
        if (winsperorder<1) winsperorder=1;
        sprintf(output, "Warning: too many windows per order requested, \
falling back to %d", winsperorder);
        SCTPUT(output);
    }
    if (winsperorder>1) {
        if ((xminwidth/ordpos->step[0])<2) {
            xminwidth=2*ordpos->step[0];
            sprintf(output, "Warning: requested window x size is too small, \
increasing to %f", xminwidth);
            SCTPUT(output);
        }
    }
    else {
        xminwidth=ordpos->step[0];
        sprintf(output, "Warning: setting window x size to %f", xminwidth);
        SCTPUT(output);
    }
    if (xminwidth>((double)ordpos->npix[0]*ordpos->step[0]/(double)winsperorder)){
        xminwidth = (double)ordpos->npix[0]*ordpos->step[0]/(double)winsperorder;
        sprintf(output, "Warning: requested window x size is too large, reducing to %f",
                xminwidth);
        SCTPUT(output);
    }
    if ((yminwidth/ordpos->step[1])<1) {
        yminwidth=ordpos->step[1];
        sprintf(output, "Warning: requested window y size is too small, \
increasing to %f", yminwidth);
        SCTPUT(output);
    }

    /* setting up some parameters */
    halfyminwidth = yminwidth/2;
    ymax = ordpos->start[1]+(double)ordpos->npix[1]*ordpos->step[1];
    halfxminwidth = xminwidth/2;
    firstx = ordpos->start[0]+halfxminwidth;
    winstep = ordpos->step[0]*(double)ordpos->npix[0]/(double)winsperorder;
    if (winstep<xminwidth) winstep=xminwidth;
    xcentres = dvector(0, winsperorder);
    xstarts = dvector(0, winsperorder);
    xends = dvector(0, winsperorder);
    /* compute the x boundaries once and for all, they do not change */
    double xcentre=0;
    for (ixcount=0; ixcount<=winsperorder-1; ixcount++) {
        xcentres[ixcount] = xcentre = firstx+(double)ixcount*winstep;
        xstarts[ixcount] = xcentre-halfxminwidth;
        xends[ixcount] = xcentre+halfxminwidth;
    }

    /* ok, now we can proceed to create the background table */

    TCTINI(backtab, F_O_MODE, (ordpos->lastorder+2-ordpos->firstorder)*winsperorder,
           &tid);
    TCCINI(tid, D_R8_FORMAT, 1, "G14.7", "", "ORDER", &ordcol);
    TCCINI(tid, D_R8_FORMAT, 1, "G14.7", "", "X", &xcol);
    TCCINI(tid, D_R8_FORMAT, 1, "G14.7", "", "XSTA", &xstartcol);
    TCCINI(tid, D_R8_FORMAT, 1, "G14.7", "", "XEND", &xendcol);
    TCCINI(tid, D_R8_FORMAT, 1, "G14.7", "", "YBKG", &ycol);
    TCCINI(tid, D_R8_FORMAT, 1, "G14.7", "", "YSTA", &ystartcol);
    TCCINI(tid, D_R8_FORMAT, 1, "G14.7", "", "YEND", &yendcol);

    windowcounter = 0;
    for (iorder=ordpos->firstorder-1; iorder<=ordpos->lastorder; iorder++) {
        dorder = iorder+.5;
        for (ixcount=0; ixcount<=winsperorder-1; ixcount++) {
            if (get_ordpos(ordpos, dorder, xstarts[ixcount], &ystart) !=NOERR) {
                SCTPUT("Error computing interorder position");
                return flames_midas_fail();
            }
            if (get_ordpos(ordpos, dorder, xcentres[ixcount], &ycentre) !=NOERR) {
                SCTPUT("Error computing interorder position");
                return flames_midas_fail();
            }
            if (get_ordpos(ordpos, dorder, xends[ixcount], &yend) !=NOERR) {
                SCTPUT("Error computing interorder position");
                return flames_midas_fail();
            }
            if (ystart<=yend) {
                ystart -= halfyminwidth;
                yend += halfyminwidth;
            }
            else {
                yswap = yend-halfyminwidth;
                yend = ystart+halfyminwidth;
                ystart = yswap;
            }
            /* do we need to trim the window? */
            if (ystart<ordpos->start[1]) ystart=ordpos->start[1];
            if (yend>ymax) yend=ymax;
            /* now, before writing a given background window to the back table,
	 check that it falls indeed into the frame */
            if (((yend-ystart)/ordpos->step[1]>1) && (ycentre>ystart) &&
                            (ycentre<yend)) {
                /* do write this window */
                windowcounter++;
                TCEWRD(tid, windowcounter, ordcol, &dorder);
                TCEWRD(tid, windowcounter, xcol, &xcentres[ixcount]);
                TCEWRD(tid, windowcounter, xstartcol, &xstarts[ixcount]);
                TCEWRD(tid, windowcounter, xendcol, &xends[ixcount]);
                TCEWRD(tid, windowcounter, ycol, &ycentre);
                TCEWRD(tid, windowcounter, ystartcol, &ystart);
                TCEWRD(tid, windowcounter, yendcol, &yend);
            }
        }
    }

    TCTCLO(tid);

    /* free dynamic arrays */
    free_dvector(xcentres, 0, winsperorder);
    free_dvector(xstarts, 0, winsperorder);
    free_dvector(xends, 0, winsperorder);
    if (freeordpos(ordpos)!=NOERR) {
        SCTPUT("Warning: error freeing the orderpos structure");
    }
    free(ordpos);

    if (windowcounter==0) {
        SCTPUT("Error: no background windows written");
        return flames_midas_fail();
    }

    return SCSEPI();
}

/**@}*/