/*                                                                              *
 *   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       *
 *                                                                              */

/*
 * $Author: amodigli $
 * $Date: 2013-10-11 12:43:10 $
 * $Revision: 1.17 $
 * $Name: not supported by cvs2svn $
 * $Log: not supported by cvs2svn $
 * Revision 1.15  2012/03/02 16:49:48  amodigli
 * fixed warning related to upgrade to CPL6
 *
 * Revision 1.14  2011/12/08 14:05:48  amodigli
 * Fix warnings with CPL6
 *
 * Revision 1.13  2010/09/24 09:32:07  amodigli
 * put back QFITS dependency to fix problem spot by NRI on FIBER mode (with MIDAS calibs) data
 *
 * Revision 1.11  2007/06/06 08:17:33  amodigli
 * replace tab with 4 spaces
 *
 * Revision 1.10  2007/04/24 12:50:29  jmlarsen
 * Replaced cpl_propertylist -> uves_propertylist which is much faster
 *
 * Revision 1.9  2007/04/10 07:09:37  jmlarsen
 * Changed interface of uves_spline_hermite()
 *
 * Revision 1.8  2006/11/06 15:19:41  jmlarsen
 * Removed unused include directives
 *
 * Revision 1.7  2006/08/17 13:56:53  jmlarsen
 * Reduced max line length
 *
 * Revision 1.6  2006/05/12 15:10:07  jmlarsen
 * Shortened lines
 *
 * Revision 1.5  2006/04/10 12:38:13  jmlarsen
 * Bugfix: don't read uninitialized memory (caused atmospheric extinction step 
 * to be randomly disabled)
 *
 * Revision 1.4  2006/04/06 08:49:23  jmlarsen
 * Propagate errors when normalizing spectrum
 *
 * Revision 1.3  2005/12/19 16:17:56  jmlarsen
 * Replaced bool -> int
 *
 * Revision 1.2  2005/12/16 14:22:23  jmlarsen
 * Removed midas test data; Added sof files
 *
 * Revision 1.1  2005/11/11 13:18:54  jmlarsen
 * Reorganized code, renamed source files
 *
 */

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

/*----------------------------------------------------------------------------*/
/**
 * @addtogroup uves_reduce
 */
/*----------------------------------------------------------------------------*/
/**@{*/

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

#include <uves_reduce_utils.h>
#include <math.h>
#include <uves_pfits.h>
#include <uves_utils.h>
#include <uves_utils_wrappers.h>
#include <uves_error.h>
#include <uves_msg.h>

#include <cpl.h>

/*-----------------------------------------------------------------------------
                            Functions prototypes
 -----------------------------------------------------------------------------*/

/*-----------------------------------------------------------------------------
                            Implementation
 -----------------------------------------------------------------------------*/
static double
uves_get_airmass(const double z) {
    /* http://en.wikipedia.org/wiki/Air_mass_%28astronomy%29 */

      /* X= sec z * ( 1-0.0012( (sec z)^2 -1 ) ) */
    double x=z/180.*CPL_MATH_PI;

    double sec=1./cos(x);
    return ( sec * ( 1.-0.0012 * ( sec*sec-1. ) ) );

    /* Yung 1994
    double A= 1.002432*cos(x)*cos(x)+0.148386*cos(z)+0.0096467;
    double B= cos(x)*cos(x)*cos(x)+0.149864*cos(x)*cos(x)+0.0102963*cos(x)+0.000303978;
    return A/B;
    */


}

/*----------------------------------------------------------------------------*/
/**
   @brief    Normalize a spectrum
   @param    spectrum           The 1d (merged) or 2d (non-merged or 2d
                                extracted+merged) spectrum to be normalized
   @param    spectrum_error     Error (1 sigma) of @em spectrum, or NULL.
   @param    spectrum_header    Header describing the geometry of the input spectrum
   @param    raw_header         Header of the raw frame
   @param    n_traces           Number of spatial bins (1 unless 2d extracted)
   @param    chip               CCD chip
   @param    atm_extinction     The table of extinction coefficients
   @param    correct_binning    Flag indicating whether or not to divide by
                                the x-binning factor
   @param    scaled_error       (output) If non-NULL, error of output spectrum
   @return   The normalized spectrum

   The spectrum is divided by exposure time, gain and (optionally) binning. 
   Also corrected for atmospheric extinction using the provided table of
   extinction coefficients.

   Bad pixels are propagated.
*/
/*----------------------------------------------------------------------------*/
cpl_image *
uves_normalize_spectrum(const cpl_image *spectrum, const cpl_image *spectrum_error,
            const uves_propertylist *spectrum_header,
            const uves_propertylist *raw_header,
            int n_traces,
            enum uves_chip chip,
            const cpl_table *atm_extinction,
            bool correct_binning,
            cpl_image **scaled_error)
{
    cpl_image *scaled = NULL;
    double exptime, gain;
    int binx;
    int norders, ny, nx;
    
    assure_nomsg( spectrum != NULL, CPL_ERROR_NULL_INPUT);
    assure_nomsg( scaled_error == NULL || spectrum_error != NULL, CPL_ERROR_NULL_INPUT);
    assure_nomsg( spectrum_header != NULL, CPL_ERROR_NULL_INPUT);

    nx = cpl_image_get_size_x(spectrum);
    ny = cpl_image_get_size_y(spectrum);

    if (spectrum_error != NULL)
    {
        assure( nx == cpl_image_get_size_x(spectrum_error) &&
            ny == cpl_image_get_size_y(spectrum_error), CPL_ERROR_INCOMPATIBLE_INPUT,
            "Error spectrum geometry differs from spectrum: %" CPL_SIZE_FORMAT "x%" CPL_SIZE_FORMAT " vs. %d x %d",
            cpl_image_get_size_x(spectrum_error),
            cpl_image_get_size_y(spectrum_error),
            nx, ny);
    }
    
    assure( ny % n_traces == 0, CPL_ERROR_INCOMPATIBLE_INPUT,
        "Spectrum image height (%d) is not a multiple of "
        "the number of traces (%d). Confused, bailing out",
        ny, n_traces);
    
    norders = ny / n_traces;

    /*
     * Correct for exposure time, gain, bin 
     */
    check( exptime = uves_pfits_get_exptime(raw_header),
           "Could not read exposure time");
    
    check( gain = uves_pfits_get_gain(raw_header, chip),
       "Could not read gain factor");
    
    if (correct_binning)
    {
        /* x-binning of rotated image is y-binning of raw image */
        check( binx  = uves_pfits_get_biny(raw_header),
           "Could not read binning");
    }
    else
    {
        uves_msg("Spectrum will not be normalized to unit binning");
        binx = 1;
    }
    
    assure( exptime > 0, CPL_ERROR_ILLEGAL_INPUT, "Non-positive exposure time: %f s", exptime);
    assure( gain    > 0, CPL_ERROR_ILLEGAL_INPUT, "Non-positive gain: %f", gain);
    assure( binx    > 0, CPL_ERROR_ILLEGAL_INPUT, "Illegal binning: %d", binx);
    
    uves_msg("Correcting for exposure time = %f s, gain = %f, binx = %d", exptime, gain, binx);
    
    check( scaled        = cpl_image_divide_scalar_create(spectrum, exptime * gain * binx),
       "Error correcting spectrum for gain, exposure time, binning");
    
    if (scaled_error != NULL)
    {
        check( *scaled_error = cpl_image_divide_scalar_create(spectrum_error, 
                                  exptime * gain * binx),
           "Error correcting rebinned spectrum for gain, exposure time, binning");
    }
    
    /* 
     * Correct for atmospheric extinction 
     */
    {
    double airmass;
    double dlambda, lambda_start;
    int order;
    int miss_airm_start=0;
    int miss_airm_end=0;
    {
        double airmass_start, airmass_end;
        if(cpl_propertylist_has(raw_header,UVES_AIRMASS_START)) {
            check( airmass_start = uves_pfits_get_airmass_start(raw_header),
                            "Error reading airmass start");
        } else {
            double alt = uves_pfits_get_tel_alt_start(raw_header);
            double z=90.-alt;
            miss_airm_start=1;

            airmass_start=uves_get_airmass(z);
            uves_msg_warning("Missing %s approximate it", UVES_AIRMASS_START);

        }
        if(cpl_propertylist_has(raw_header,UVES_AIRMASS_END)) {
            check( airmass_end = uves_pfits_get_airmass_end(raw_header),
                            "Error reading airmass end");
        } else {
            double alt = uves_pfits_get_tel_alt_start(raw_header);
            double z=90.-alt;
            miss_airm_end=1;
            airmass_end=uves_get_airmass(z);
            uves_msg_warning("Missing %s approximate it", UVES_AIRMASS_END);

        }
        if(miss_airm_start && miss_airm_end) {
            uves_msg_error("Both airmass start and end aere missing. Stop");
            cpl_error_set(cpl_func,CPL_ERROR_ILLEGAL_INPUT);
        }  else {
            /* Use arithmetic mean of airmass start/end */
            airmass = (airmass_start + airmass_end) / 2;
        }

    }

    uves_msg("Correcting for extinction through airmass %f", airmass);
        
    check( dlambda = uves_pfits_get_cdelt1(spectrum_header),
           "Error reading bin width from header");
    
    for (order = 1; order <= norders; order++)
        {
        int trace;

        /* If spectrum was already merged, then read crval1,
         * otherwise read wstart for each order
         */
        
        if (norders == 1)
            {
            check( lambda_start = uves_pfits_get_crval1(spectrum_header),
                   "Error reading start wavelength from header");    
            }
        else
            {
            check( lambda_start = uves_pfits_get_wstart(spectrum_header, order),
                   "Error reading start wavelength from header");    
            }

        for (trace = 1; trace <= n_traces; trace++)
            {
            int spectrum_row = (order - 1)*n_traces + trace;
            int x;
            
            for (x = 1; x <= nx; x++)
                {
                int pis_rejected1;
                int pis_rejected2;
                double flux;
                double dflux = 0;
                double extinction;
                double lambda;

                lambda = lambda_start + (x-1) * dlambda;
                
                flux  = cpl_image_get(scaled, x, spectrum_row, &pis_rejected1);
                if (scaled_error != NULL)
                    {
                    dflux = cpl_image_get(*scaled_error, x, 
                                  spectrum_row, &pis_rejected2);
                    }

                if (!pis_rejected1 && (scaled_error == NULL || !pis_rejected2))
                    {
                                        int istart = 0;

                    /* Read extinction (units: magnitude per airmass) */
                    check( extinction = 
                           uves_spline_hermite_table(
                           lambda, atm_extinction,
                           "LAMBDA", "LA_SILLA", &istart),
                           "Error interpolating extinction coefficient");
                    
                    /* Correct for extinction using
                     * the magnitude/flux relation
                     * m = -2.5 log_10 F
                     *  => 
                     * F = 10^(-m*0.4)
                     *
                     * m_top-of-atmosphere = m - ext.coeff*airmass
                     * F_top-of-atmosphere = F * 10^(0.4 * ext.coeff*airmass)
                     */

                    cpl_image_set(
                        scaled, x, spectrum_row,
                        flux * pow(10, 0.4 * extinction * airmass));
                    if (scaled_error != NULL)
                        {
                        cpl_image_set(
                            *scaled_error, x, spectrum_row,
                            dflux * pow(10, 0.4 * extinction * airmass));
                        }
                    }
                else
                    {
                    cpl_image_reject(scaled, x, spectrum_row);
                    if (scaled_error != NULL)
                        {
                        cpl_image_reject(*scaled_error, x, spectrum_row);
                        }
                    }
                } /* for each x */

            } /* for each (possibly only 1) trace */

        } /* for each (possibly only 1) order */
    }
    
  cleanup:
    if (cpl_error_get_code() != CPL_ERROR_NONE)
    {
        uves_free_image(&scaled);
        if (scaled_error != NULL)
        {
            uves_free_image(scaled_error);
        }
    }
    
    return scaled;
}

/**@}*/