// bytran -|- spectral calculations for portable devices using the HITRAN database
//
// Developed by Dzianis Pliutau (Denis Plutov) and Kostya Roslyakov
// www.bytran.org | bytran@bytran.org
// Copyright (c) 2016 Dzianis Pliutau
// Distributed under the MIT license
//
// The instrument function calculation algorithms are converted from the Python code of the HITRAN Application
// Programming Interface (HAPI) distributed under the Academic Free License 3.0, http://hitran.org/hapi/
// R.V. Kochanov, I.E. Gordon, L.S. Rothman, P. Wcisło, C. Hill, J.S. Wilzewski, "HITRAN Application Programming Interface (HAPI):
// A comprehensive approach to working with spectroscopic data", JQSRT, Volume 177, July 2016, Pages 15–30; (hitran.org/hapi/)
//
// The MIT License (MIT)
// Copyright (c) 2016 Dzianis Pliutau
//
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy,
// modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software
// is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
// LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR
// IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.



// Object containing modules to carry out instrument function related calculations
// This module is based on the Python code in the HITRAN API programming interface
#include "instrumentfunction.h"
#include "auxiliary.h"
#include "const.h"
#include "math.h"
#include <QMessageBox>

InstrumentFunction::InstrumentFunction(QObject *parent) : QObject(parent)
{
    this->instrFunctIndex = 3;
    this->instrResol = 3;  // The derault is 0.1 cm-1 in HAPI
    this->instrFunctWing = 10;  // The default is 10 cm-1
}

// Methods to get and set the instrument function to use
quint8 InstrumentFunction::getInstrFunctIndex() { return this->instrFunctIndex; }
void InstrumentFunction::setInstrFunctIndex(quint8 newInstrFunctIndex) { this->instrFunctIndex = newInstrFunctIndex; }

// Geting and seting the instrument function resolution
double InstrumentFunction::getInstrResol() { return this->instrResol; }
void InstrumentFunction::setInstrResol(double newInstrResol) { this->instrResol = newInstrResol; }

// Geting and seting the instrument function wing
double InstrumentFunction::getInstrFunctWing() { return this->instrFunctWing; }
void InstrumentFunction::setInstrFunctWing(double newInstrFunctWing) { this->instrFunctWing = newInstrFunctWing; }

// Method to set the wavelength object
void InstrumentFunction::setWavelength(Wavelength *newWavelength) { this->wavelength = newWavelength; }

// Method to calculate the convolustion of the spectrum with the instrument funation
// Convolution of the spectrum is based in part on: http://www.cs.cornell.edu/courses/cs1114/2013sp/sections/s06_convolution.pdf
void InstrumentFunction::convolveSpectrum(Spectrum *spectrum, Spectrum *convolvedSpec)
{
    // Instrumental function
    // 0 - None
    // 1 - Rectangular (Boxcar)
    // 2 - Triangular
    // 3 - Gaussian
    // 4 - Dispersion (Lorentz)
    // 5 - Diffraction
    // 6 - Michelson

    QVector<double> instrFunct;  // Instrument function array
    QVector<double> x;   // Array with the wavenumber values for the slit function
    double step;  // Wavenumber step
    double sum;   // sum of the values in the slit function array
    double wn;    // Auxiliary wavenumber variable
    int i;  // Loop variable


    // Handle the case when the instrument function wing is wider than the spectral range specified or close to it (maybe let user know).

    if (wavelength->getWns()->count() > 2)
    {
        // Forming the wavenumber array for the slit function calculation
        x.resize(0);
        step = wavelength->getWnsAt(1) - wavelength->getWnsAt(0);

        // Check to make sure the wavenumber step is smaller than the instrument resolution
        if (step < this->instrResol)
        {
            i = 0;
            wn = (-1) * this->instrFunctWing;
            while (wn <= (this->instrFunctWing + step))
            {
                wn = (-1) * this->instrFunctWing + step * i;
                x.append(wn);
                i++;
            }

            switch (this->instrFunctIndex)
            {
                case 0:
                    // Copy properly, convolved spectrum is the same as the original spectrum convolvedSpec = spectrum;
                    Auxiliary::copySpectrum(spectrum, convolvedSpec);
                    break;
                case 1:  // 1 - Rectangular (Boxcar)
                    slitRectangular(&x, this->instrResol, &instrFunct);  // Rectangular slit instrument function
                    break;
                case 2:  // 2 - Triangular
                    slitTriangular(&x, this->instrResol, &instrFunct);   // Triangular slit function
                    break;
                case 3:  // 3 - Gaussian
                    slitGaussian(&x, this->instrResol, &instrFunct);     // Gaussian slit function
                    break;
                case 4:  // 4 - Dispersion (Lorentz)
                    slitDispersion(&x, this->instrResol, &instrFunct);   // Dispersion slit function
                    break;
                case 5:  // 5 - Diffraction
                    slitDiffraction(&x, this->instrResol, &instrFunct);  // Diffraction slit function
                    break;
                case 6:  // 6 - Michelson
                    slitMichelson(&x, this->instrResol, &instrFunct);    // Michelson slit function
                    break;
                case 7:  // 7 - Cosine
                    slitCosine(&x, this->instrResol, &instrFunct);       // Cosine instrument function
                default:
                    break;
            }

            // Normalization of the slit function
            // slit /= sum(slit) * step # simple normalization
            sum = 0;
            for (i = 0; i < instrFunct.count(); i++) { sum = sum + instrFunct[i]; }
            for (i = 0; i < instrFunct.count(); i++)
            {
                instrFunct[i] = (instrFunct[i] / sum);  // In the HAPI there is also the "step" term for some reason
            }

            // Convolving the spectrum
            convolveSpectrum(spectrum->getY(), &(instrFunct), convolvedSpec->getY());
        }
        else
        {
            emit openOkDialog("The line-by-line calculation spectral\nresolution has to be smaller than the\ninstrument resolution.\n\nInstrument function calculation\nwill not be performed.");
        }
    }
}


// Method to calculate the convolustion of the spectrum with the instrument function
// Convolution of the spectrum is based in part on: http://www.cs.cornell.edu/courses/cs1114/2013sp/sections/s06_convolution.pdf
void InstrumentFunction::convolveSpectrum(QVector<double> *spectrum, QVector<double> *instrFunct, QVector<double> *convolvedSpec)
{
    int i;  // Loop variable
    int j;  // Loop variable
    int halfLength;  // Half length of the instrument function
    double sum;  // Sum of values along the kermel (instrument function bin)


    // Determining the half-distance from the center
    halfLength = static_cast<int> (floor(instrFunct->count() * 0.5));

    // Resizing the convolved spectrum
    convolvedSpec->resize(spectrum->count());

    // Processing the convolution if the length of the array larger than that of the edge regions
    if (instrFunct->length() < spectrum->count())
    {
        // Writing zero values into the edje regions of the convolved spectrum
        for (i = 0; i <= (halfLength + 1); i++)
        {
            (*convolvedSpec)[i] = 0;
            (*convolvedSpec)[convolvedSpec->count() - 1 - i] = 0;
        }

        // Forming the convolved spectrum
        for (i = halfLength; i < spectrum->count() - halfLength; i++)
        {
            sum = 0;
            for (j = 0; j < instrFunct->count(); j++)
            {
                sum = sum + (*spectrum)[i - halfLength + j] * (*instrFunct)[j];
            }
            // Saving the value of the convolved spectrum
            (*convolvedSpec)[i] = sum;
        }
    }
}

// Rectangular slit instrument function
void InstrumentFunction::slitRectangular(QVector<double> *x, double g, QVector<double> *instrFunct)
{
    // Instrumental (slit) function.
    // B(x) = 1/γ , if |x| ≤ γ/2 & B(x) = 0, if |x| > γ/2,
    // where γ is a slit width or the instrumental resolution.
    // """
    // index_inner = abs(x) <= g/2
    // index_outer = ~index_inner
    // y = zeros(len(x))
    // y[index_inner] = 1/g
    // y[index_outer] = 0
    // return y

    double halfg;  // Half of the resolution
    int i;  // Loop variable


    // Calculating the instrument function
    halfg = g * 0.5;
    instrFunct->resize(x->count());
    for (i = 0; i < instrFunct->count(); i++)
    {
        if (fabs((*x)[i]) <= halfg) { (*instrFunct)[i] = 1 / g; }
    }
}

// Triangular slit function
void InstrumentFunction::slitTriangular(QVector<double> *x, double g, QVector<double> *instrFunct)
{
    // Instrumental (slit) function.
    // B(x) = 1/γ*(1-|x|/γ), if |x| ≤ γ & B(x) = 0, if |x| > γ,
    // where γ is the line width equal to the half base of the triangle.
    // """
    // index_inner = abs(x) <= g
    // index_outer = ~index_inner
    // y = zeros(len(x))
    // y[index_inner] = 1/g * (1 - abs(x[index_inner])/g)
    // y[index_outer] = 0
    // return y

    int i;  // Loop variable

    // Calculating the instrument function
    instrFunct->resize(x->count());
    for (i = 0; i < instrFunct->count(); i++)
    {
        if (fabs((*x)[i]) <= g) { (*instrFunct)[i] = 1 / g * (1 - fabs((*x)[i]) / g); }
    }
}

// Gaussian slit function
void InstrumentFunction::slitGaussian(QVector<double> *x, double g, QVector<double> *instrFunct)
{
    // Instrumental (slit) function.
    // B(x) = sqrt(ln(2)/pi)/γ*exp(-ln(2)*(x/γ)**2),
    // where γ/2 is a Gaussian half-width at half-maximum.
    // """
    // g /= 2
    // return sqrt(log(2))/(sqrt(pi)*g)*exp(-log(2)*(x/g)**2)

    int i;  // Loop variable

    // Calculating the instrument function
    g = g * 0.5;
    instrFunct->resize(x->count());
    for (i = 0; i < instrFunct->count(); i++)
    {
        (*instrFunct)[i] = sqrt(log(2)) / (sqrt(Const::pi)*g) * exp(-log(2) * ((*x)[i]/g) * ((*x)[i]/g));
    }
}

// Dispersion slit function
void InstrumentFunction::slitDispersion(QVector<double> *x, double g, QVector<double> *instrFunct)
{
    // Instrumental (slit) function.
    // B(x) = γ/pi/(x**2+γ**2),
    // where γ/2 is a Lorentzian half-width at half-maximum.
    // """
    // g /= 2
    // return g/pi/(x**2+g**2)

    int i;  // Loop variable

    // Calculating the instrument function
    g = g * 0.5;
    instrFunct->resize(x->count());
    for (i = 0; i < instrFunct->count(); i++)
    {
        (*instrFunct)[i] = g / Const::pi / ((*x)[i] * (*x)[i] + g * g);
    }
}

// Cosine slit function
void InstrumentFunction::slitCosine(QVector<double> *x, double g, QVector<double> *instrFunct)
{
    // This is from the HAPI's Python code
    // return (cos(pi/g*x)+1)/(2*g);

    int i;  // Loop variable

    instrFunct->resize(x->count());
    for (i = 0; i < instrFunct->count(); i++)
    {
        (*instrFunct)[i] = (cos(Const::pi / (g * (*x)[i])) + 1) / (2 * g);
    }
}

// Diffraction slit function
void InstrumentFunction::slitDiffraction(QVector<double> *x, double g, QVector<double> *instrFunct)
{
    // Instrumental (slit) function.
    // """
    // y = zeros(len(x))
    // index_zero = x==0
    // index_nonzero = ~index_zero
    // dk_ = pi/g
    // x_ = dk_*x[index_nonzero]
    // w_ = sin(x_)
    // r_ = w_**2/x_**2
    // y[index_zero] = 1
    // y[index_nonzero] = r_/g
    // return y


    // Auxiliary calculation variables
    double dk_;
    double x_;
    double w_;
    double r_;
    int midIndex;
    int i;  // Loop variable


    // Resizing the instrument function
    instrFunct->resize(x->count());

    // Setting the center value to be equal to 1
    midIndex = static_cast<int> (ceil(instrFunct->count() * 0.5));
    (*instrFunct)[midIndex] = 1;

    // Calculating the instrument function
    dk_ = Const::pi / g;
    for (i = 0; i < midIndex; i++)
    {
        x_ = dk_ * (*x)[i];
        w_ = sin(x_);
        r_ = (w_* w_) / (x_ * x_);
        (*instrFunct)[i] = r_ / g;
    }
    for (i = midIndex + 1; i < (*instrFunct).count(); i++)
    {
        x_ = dk_ * (*x)[i];
        w_ = sin(x_);
        r_ = (w_* w_) / (x_ * x_);
        (*instrFunct)[i] = r_ / g;
    }
}

// Michelson slit function
void InstrumentFunction::slitMichelson(QVector<double> *x, double g, QVector<double> *instrFunct)
{
    // Instrumental (slit) function.
    // B(x) = 2/γ*sin(2pi*x/γ)/(2pi*x/γ) if x!=0 else 1,
    // where 1/γ is the maximum optical path difference.
    // """
    // y = zeros(len(x))
    // index_zero = x==0
    // index_nonzero = ~index_zero
    // dk_ = 2*pi/g
    // x_ = dk_*x[index_nonzero]
    // y[index_zero] = 1
    // y[index_nonzero] = 2/g*sin(x_)/x_
    // return y

    // Auxiliary calculation variables
    double dk_;
    double x_;
    int midIndex;
    int i;  // Loop variable


    // Resizing the instrument function
    instrFunct->resize(x->count());

    // Setting the center value to be equal to 1
    midIndex = static_cast<int> (floor(instrFunct->count() * 0.5));
    (*instrFunct)[midIndex] = 1;

    // Calculating the instrument function
    dk_ = 2 * Const::pi / g;
    for (i = 0; i < midIndex; i++)
    {
        x_ = dk_ * (*x)[i];
        (*instrFunct)[i] = 2 / g * sin(x_) / x_;
    }
    for (i = midIndex + 1; i < (*instrFunct).count(); i++)
    {
        x_ = dk_ * (*x)[i];
        (*instrFunct)[i] = 2 / g * sin(x_) / x_;
    }
}


// HAPI example
// # spectral convolution with an apparatus (slit) function
// def convolveSpectrum(Omega,CrossSection,Resolution=0.1,AF_wing=10.,SlitFunction=SLIT_RECTANGULAR):
//    """
//    INPUT PARAMETERS:
//        Omega:         wavenumber grid                           (required)
//        CrossSection:  high-res cross section calculated on grid (required)
//        Resolution:    instrumental resolution γ                 (optional)
//        AF_wing:       instrumental function wing                (optional)
//        SlitFunction:  instrumental function for low-res spectra calculation (optional)
//    OUTPUT PARAMETERS:
//        Omega: wavenumber grid
//        CrossSection: low-res cross section calculated on grid
//        i1: lower index in Omega input
//        i2: higher index in Omega input
//        slit: slit function calculated over grid [-AF_wing; AF_wing]
//                with the step equal to instrumental resolution.
//    ---
//    DESCRIPTION:
//        Produce a simulation of experimental spectrum via the convolution
//        of a “dry” spectrum with an instrumental function.
//        Instrumental function is provided as a parameter and
//        is calculated in a grid with the width=AF_wing and step=Resolution.
//    ---
//    EXAMPLE OF USAGE:
//        nu_,radi_,i,j,slit = convolveSpectrum(nu,radi,Resolution=2.0,AF_wing=10.0,
//                                                SlitFunction=SLIT_MICHELSON)
//    ---
//    """
//    step = Omega[1]-Omega[0]
//    if step>=Resolution: raise Exception('step must be less than resolution')
//    x = arange(-AF_wing,AF_wing+step,step)
//    slit = SlitFunction(x,Resolution)
//    # FIXING THE BUG: normalize slit function
//    slit /= sum(slit)*step # simple normalization
//    left_bnd = len(slit)/2
//    right_bnd = len(Omega) - len(slit)/2
//    #CrossSectionLowRes = convolve(CrossSection,slit,mode='valid')*step
//    CrossSectionLowRes = convolve(CrossSection,slit,mode='same')*step
//    #return Omega[left_bnd:right_bnd],CrossSectionLowRes,left_bnd,right_bnd,slit
//    return Omega[left_bnd:right_bnd],CrossSectionLowRes[left_bnd:right_bnd],left_bnd,right_bnd,slit
