// bytran -|- spectral calculations for portable devices using the HITRAN database
//
// Developed by Dzianis Pliutau (Denis Plutov) and Kostya Roslyakov
// www.bytran.org | contact@bytran.org
// Copyright (c) 2016 Dzianis Pliutau
// Distributed under the MIT license
//
// Parts of the line-by-line calculation algorithm are based on the Python code from 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/)
//
// Appendix 1.1 of the theoretical section from the HITRAN-PC 4.0 manual was also consulted, www.cas.usf.edu/lidarlab/hitranpc/Hitran-PC_4.0_Manual.pdf
// Killinger D.K., Wilcos W.E., Pliutau D., \"HITRAN=PC: 25 years of academic development and commercialization of laser
// atmospheric transmission software for environmental remote sensing\", Technology & Innovation, Vol. 14, No. 3-4, pp. 303-327
//
// The Voigt calculation routine uses the algorithm described in: Abrarov, S., and Quine, B. (2015).
// "A rational approximation for efficient computation of the Voigt function in quantitative spectroscopy".
// Journal of Mathematics Research, 7(2), p 163. doi:http://dx.doi.org/10.5539/jmr.v7n2.
//
// 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 with all line-by-line calculation elements
#include <math.h>
#include <complex>
#include <QVector>
#include <QMessageBox>
#include <QApplication>
#include <QFile>
#include <QTextStream>
#include <QApplication>
#include "lbl.h"
#include "auxiliary.h"
#include "hitranonerec.h"
#include "const.h"


Lbl::Lbl(QObject *parent, HitranOnlineDownloader *newHitran) : QObject(parent)
{
    // Passing in a reference to the HITRAN download object
    this->hitran = newHitran;

    // The following coefficents used in the the Voigt profile implementation are taken from the publication reported in:
    // Abrarov, S., & Quine, B. (2015). A Rational Approximation for Efficient Computation of the Voigt Function
    // in Quantitative Spectroscopy. Journal Of Mathematics Research, 7(2), p163. doi:http://dx.doi.org/10.5539/jmr.v7n2
    coeff1_12.resize(12);
    coeff2_12.resize(12);
    coeff3_12.resize(12);
    coeff1_12[0] = 2.307372754308023e-001; coeff2_12[0] = 4.989787261063716e-002; coeff3_12[0] = 1.464495070025765e+000;
    coeff1_12[1] = 7.760531995854886e-001; coeff2_12[1] = 4.490808534957343e-001; coeff3_12[1] = -3.230894193031240e-001;
    coeff1_12[2] = 4.235506885098250e-002; coeff2_12[2] = 1.247446815265929e+000; coeff3_12[2] = -5.397724160374686e-001;
    coeff1_12[3] = -2.340509255269456e-001; coeff2_12[3] = 2.444995757921221e+000; coeff3_12[3] = -6.547649406082363e-002;
    coeff1_12[4] = -4.557204758971222e-002; coeff2_12[4] = 4.041727681461610e+000; coeff3_12[4] = 2.411056013969393e-002;
    coeff1_12[5] = 5.043797125559205e-003; coeff2_12[5] = 6.037642585887094e+000; coeff3_12[5] = 4.001198804719684e-003;
    coeff1_12[6] = 1.180179737805654e-003; coeff2_12[6] = 8.432740471197681e+000; coeff3_12[6] = -5.387428751666454e-005;
    coeff1_12[7] = 1.754770213650354e-005; coeff2_12[7] = 1.122702133739336e+001; coeff3_12[7] = -2.451992671326258e-005;
    coeff1_12[8] = -3.325020499631893e-006; coeff2_12[8] = 1.442048518447414e+001; coeff3_12[8] = -5.400164289522879e-007;
    coeff1_12[9] = -9.375402319079375e-008; coeff2_12[9] = 1.801313201244001e+001; coeff3_12[9] = 1.771556420016014e-008;
    coeff1_12[10] = 8.034651067438904e-010; coeff2_12[10] = 2.200496182129099e+001; coeff3_12[10] = 4.940360170163906e-010;
    coeff1_12[11] = 3.355455275373310e-011; coeff2_12[11] = 2.639597461102705e+001; coeff3_12[11] = 5.674096644030151e-014;

    coeff1_16.resize(16);
    coeff2_16.resize(16);
    coeff3_16.resize(16);
    coeff1_16[0] = 1.608290174437121e-001; coeff2_16[0] = 3.855314219175531e-002; coeff3_16[0] = 1.366578214428949e+000;
    coeff1_16[1] = 6.885967427017463e-001; coeff2_16[1] = 3.469782797257978e-001; coeff3_16[1] = -5.742919588559361e-002;
    coeff1_16[2] = 2.651151642675390e-001; coeff2_16[2] = 9.638285547938826e-001; coeff3_16[2] = -5.709602545656873e-001;
    coeff1_16[3] = -2.050008245317253e-001; coeff2_16[3] = 1.889103967396010e+000; coeff3_16[3] = -2.011075414803758e-001;
    coeff1_16[4] = -1.274551644219086e-001; coeff2_16[4] = 3.122804517532180e+000; coeff3_16[4] = 1.069871368716704e-002;
    coeff1_16[5] = -1.134971805306579e-002; coeff2_16[5] = 4.664930205202391e+000; coeff3_16[5] = 1.468639542320982e-002;
    coeff1_16[6] = 4.201921570328543e-003; coeff2_16[6] = 6.515481030406647e+000; coeff3_16[6] = 1.816268776500938e-003;
    coeff1_16[7] = 8.084740485193432e-004; coeff2_16[7] = 8.674456993144942e+000; coeff3_16[7] = -6.875907999947567e-005;
    coeff1_16[8] = 1.946391440605860e-005; coeff2_16[8] = 1.114185809341728e+001; coeff3_16[8] = -2.327910355924500e-005;
    coeff1_16[9] = -4.132639863292073e-006; coeff2_16[9] = 1.391768433122366e+001; coeff3_16[9] = -1.004011418729134e-006;
    coeff1_16[10] = -2.656262492217795e-007; coeff2_16[10] = 1.700193570656409e+001; coeff3_16[10] = 2.304990232059197e-008;
    coeff1_16[11] = -1.524188131553777e-009; coeff2_16[11] = 2.039461221943855e+001; coeff3_16[11] = 2.275276345355270e-009;
    coeff1_16[12] = 2.239681784892829e-010; coeff2_16[12] = 2.409571386984707e+001; coeff3_16[12] = 3.383885053101652e-011;
    coeff1_16[13] = 4.939143128687883e-012; coeff2_16[13] = 2.810524065778962e+001; coeff3_16[13] = -4.398940326332977e-013;
    coeff1_16[14] = 4.692078138494072e-015; coeff2_16[14] = 3.242319258326621e+001; coeff3_16[14] = -1.405511706545786e-014;
    coeff1_16[15] = -2.512454984032184e-016; coeff2_16[15] = 3.704956964627684e+001; coeff3_16[15] = -3.954682293307548e-016;

    // Method to create actions
    createActions();
}

// Method to create connections
void Lbl::createActions()
{
    // Need to expand this for other profiles
    connect(this->hitran, SIGNAL(dataAvailable()), this, SLOT(calcWnSegment()));
    connect(this->hitran, SIGNAL(noMoreDataAvailable()), this, SLOT(postCalculation()));
}

// Method to get and set a pointer to the object with the calculation parameters
Params *Lbl::getParams() { return this->par; }
void Lbl::setParams(Params *newParams) { this->par = newParams; }

// Method to get and set an object with wavelengths
Wavelength *Lbl::getWavelength() { return this->wavelength; }
void Lbl::setWavelength(Wavelength *newWavelength) { this->wavelength = newWavelength; }

// Methods to get and set the HITRAN records to run calculations for
QVector<HitranOneRec> *Lbl::getHitranRecs() { return this->hitranRecs; }
void Lbl::setHitranRecs(QVector<HitranOneRec> *newHitranRecs) { this->hitranRecs = newHitranRecs; }

// Methods to get and set a pointer to the spectra array (this object is used to save the calculated spectra to)
OneLayerSpectra *Lbl::getHorizSpec() { return this->horizSpectra; }
void Lbl::setHorizSpec(OneLayerSpectra *newHorizSpectra) { this->horizSpectra = newHorizSpectra; }

// Method to pass the HITRAN downloader class reference to the object (used for segment reading of the data)
void Lbl::setHitranDownloader(HitranOnlineDownloader *newHitran) { this->hitran = newHitran; }

// Method is used to set a variable used to cancel the calculation
void Lbl::setCancelCalculation(bool newCancelCalculation)
{
    this->cancelCalculation = newCancelCalculation;
    QApplication::processEvents();
}

// Method to initialize the spectra (resize them according to the molparam.txt file information)
void Lbl::initHorizSpectra()
{
    int i;  // Loop variable
    int j;  // Loop variable
    int k;  // Loop variable


    // Initializing the horizontal path spectra
    horizSpectra->getMol()->resize(par->getMol()->count());
    for (i = 0; i < par->getMol()->count(); i++)
    {
        // Resizing the spectra according to the number of isotopes in the set of spectra
        horizSpectra->getMolAt(i)->getSpec()->getY()->resize(wavelength->getWns()->count());
        horizSpectra->getMolAt(i)->getSpec()->getX()->resize(0);
        horizSpectra->getMolAt(i)->getIso()->resize(par->getMolAt(i)->getIso()->count());
        for (j = 0; j < horizSpectra->getMolAt(i)->getIso()->count(); j++)
        {
            horizSpectra->getMolAt(i)->getIsoAt(j)->getSpec()->getY()->resize(wavelength->getWns()->count());
            horizSpectra->getMolAt(i)->getIsoAt(j)->getSpec()->getX()->resize(0);
        }
    }

    // Filling the arrays with 0 for the composite spectrum
    horizSpectra->getComp()->getY()->resize(wavelength->getWns()->count());
    horizSpectra->getComp()->getX()->resize(0);
    for (i = 0; i < wavelength->getWns()->count(); i++)
    {
        horizSpectra->getComp()->setYAt(i, 0);
    }

    // Filling arrays with 0 for the molecular and isotope spectra
    for (i = 0; i < par->getMol()->count(); i++)
    {
        // Filling the arrays with 0 for the molecular spectra
        for (k = 0; k < wavelength->getWns()->count(); k++)
        {
            horizSpectra->getMolAt(i)->getSpec()->setYAt(k, 0);
        }

        // Here the same should be done for the isotope spectra later on (if used)
    }
}

// Method to find maximum values in the spectra
void Lbl::findMinMaxValuesInSpectra()
{
    int i;  // Loop variable


    // Filling the arrays with 0 for the composite spectrum
    horizSpectra->getComp()->setMaxValue(Auxiliary::maximum(horizSpectra->getComp()->getY()));
    // horizSpectra->getComp()->setMinValue(Auxiliary::minimum(horizSpectra->getComp()->getY()));

    // Filling arrays with 0 for the molecular and isotope spectra
    for (i = 0; i < par->getMol()->count(); i++)
    {
        horizSpectra->getMolAt(i)->getSpec()->setMaxValue(Auxiliary::maximum(horizSpectra->getMolAt(i)->getSpec()->getY()));
        // horizSpectra->getMolAt(i)->getSpec()->setMinValue(Auxiliary::minimum(horizSpectra->getMolAt(i)->getSpec()->getY()));
    }
}

// Method to calculate spectra for a horizontal path
void Lbl::calculateHorizPath()
{
    // Initizlizing the wavenumber grid
    wavelength->initWavenumberGrid();

    // Initializing the spectral arrays
    initHorizSpectra();

    // Running calculation for a horizontal path layer depending on the lineshape profiles used
    performLblCalculation();

    // Finding the maximum spectra
    findMinMaxValuesInSpectra();
}

// Method to calculate spectra for a slant path
void Lbl::calculateSlantPath()
{
    ;
}

// Method to calculate spectra for a single layer using the Lorentz profile
// This implementation is based in part on the HAPI implementation, the thoretical section of the freely available manual for HITRAN-PC was also consulted.
void Lbl::performLblCalculation()
{
    // If the cancel button has been pressed, then the calculation is canceled
    if (this->cancelCalculation == true)
    {
        postCalculation();
    }
    else
    {
        // Initial request to get the data (had to implement this through a signal slog because option 3 with download every time needs to time to download HITRAN data from the server)
        currRecNum = 0;
        hitran->getNextSubrange();
    }
    QApplication::processEvents();   // Updating the processed plot
}

// Slot is executed when the calculation is over
void Lbl::postCalculation()
{
    // Hiding the notification
    emit notify("", 0, false, true);
    QApplication::processEvents();

    // Forming the composite spectrum is done in the spectraresultsdataobject.cpp to enable combining of the spectra after the calculation

    // Finding the maximum spectra
    findMinMaxValuesInSpectra();

    // Emitting a signal to update the plotted spectra
    emit updateSpectra(true, true);
    QApplication::processEvents();

    // Emiting a signal stating that the calculation is complete
    emit calculationComplete();
}

// Slot to calculate a segment (for different lineshape profiles)
void Lbl::calcWnSegment()
{
    // Running calculation for a horizontal path layer depending on the lineshape profiles used
    switch (par->getLineshapeProfile())
    {
        case 0:
            calcDopplerWnSegment(); break;
        case 1:
            calcLorentzWnSegment(); break;
        case 2:
            calcVoigtWnSegment(); break;
        case 3:
            calcVoigtWnSegmentHighRes(); break;
        default:
            calcVoigtWnSegment(); break;
    }
}

// Slot to calculate on e Lorentz wavenumber segment
void Lbl::calcLorentzWnSegment()
{
    int i;  // Loop variable
    int j;  // Loop variable
    int molIndex;      // Molecule index
    int isoIndex;      // Isotope index
    int totalNumRecs;  // Total number of records, used to display the calculation progress
    QVector<HitranOneRec> *recs;  // Array with the read in HITRAN data
    double pressRatio; // Pressure ratio (current to reference of 1 atm) used to save a precalculated value
    double tempRatio;  // Temperature ratio (current to reference of 296K)
    Spectrum *spec;    // Spectrum to contain the calculated result for a given record in the HITRAN database
    double optDepth;   // Optical depth value
    double m1;         // Multiplication factor, used to save a pre-calculated value
    double m2;         // Multiplication constant used on an intermediate stage of calculations
    double factor;     // Volume concentration multiplication factor
    double shift0;     // Shift 0
    double wnWingF;    // Lineshape wing contribution parameter
    double ch;         // Line intensity recalculation exponent containing  nominator
    double zn;         // Line intensity recalculation exponent containing  denominator
    double sigmaT;     // TIPS at the current temperature
    double sigmaTRef;  // TIPS at the reference temperature
    double gamma2;        // Gamma 0 squared, Lorentz broadening coefficient squared
    double gamma0 = 1;    // Gamma 0, Lorentz broadening coefficient
    int boundIndexLower;  // Lower bounding index
    int boundIndexUpper;  // Upper bounding index
    double lineIntens;    // Calculated line intensity
    double wnShifted;     // Wavenumber with the shift applied



    // Initializng an object to contain the spectral calculation results
    spec = new Spectrum();

    // Calculating the total number of records and the number of records after which the plot is to be updated
    totalNumRecs = hitran->getTotalNumRecsToProcess();

    // Reading the next set of indices
    if (hitran->getHitranFileUsage() == 2) currRecNum = 0;

    // Getting a pointer to the list of records read in
    recs = hitran->getHitranRecs();

    // Calculating the temperature and pressure ratio
    pressRatio = par->getHorizParams()->getTotalPress() / Const::totalPressRef;
    tempRatio = Const::tempRef / par->getHorizParams()->getTemp();

    // Calculating the multiplication constant, 100 is to convert to centimeters and 0.000001 is to convert from parts per million (ppmv)
    m1 = 100 * 0.000001 * 296.0 * Const::Loch / Const::pi;

    // This was slowing dow the calculation considerably on the desktop, mobile doesn't seem to matter
    // Incrementing the sub-segment HITRAN record number
    QApplication::processEvents();
    emit notify("Found " + QString::number(totalNumRecs) + " HITRAN records.\nRunning the line-by-line calculation.\nProcessing segment " + QString::number(currRecNum + 1) + "-" + QString::number(recs->count() + currRecNum + 1) + ".", 0, true, true);
    QApplication::processEvents();

    // Running calculations for all available spectral lines
    for (i = 0; i < recs->count(); i++)
    {
        // Incrementing the current record number
        currRecNum++;

        // Saving the molecule and isotope indices
        molIndex = ((*recs)[i]).getMolNum() - 1;
        isoIndex = ((*recs)[i]).getIsoNum() - 1;

        // Calculating for selected molecules only unless all molecules are to be included in calculations
        if ((par->getMolAt(molIndex)->getSelected() == true) || par->getCalcNotSelectedMols())
        {
            if ( ((*recs)[i]).getIsoNum() == 0 ) { isoIndex = 10; }  // Special case of the CO2 molecule

            // If the abundance is unavailable, then skipping the current calculation
            if (par->getMolAt(molIndex)->getIsoAt(isoIndex)->getTerrAbund() == 0) { continue; }

            // Depending on the HITRAN database units taking into account the multiplication factor
            if (par->getHitranUnits() == true) // Multiplication factor is equal to 1.0
                { factor = 1.0; }
            else  // Multiplication factor calculated from the pressure and temperature
                { factor = par->getHorizParams()->getTotalPress() / (9.869233e-7 * Const::BoltsCgs * par->getHorizParams()->getTemp()); }  // CGS units used, this is from HAPI

            // Pressure shift, taking parameters from the HITRAN database if the line shift is used
            if (par->getUseLineShift()) { shift0 = ((*recs)[i]).getAirPShift() * pressRatio; }
            else { shift0 = 0; }

            // Getting all environment dependences from Voigt parameters, line intensity, used to be a separate procedure environmentDependencyIntensity
            ch = exp(-Const::c2 * ((*recs)[i]).getLowE() / par->getHorizParams()->getTemp()) * (1 - exp(-Const::c2 * ((*recs)[i]).getWn() / par->getHorizParams()->getTemp()));
            zn = exp(-Const::c2 * ((*recs)[i]).getLowE() / Const::tempRef) * (1 - exp(-Const::c2 * ((*recs)[i]).getWn() / Const::tempRef));
            // Partition functions for T and Tref, need to be optimized
            sigmaTRef = par->getMolAt(molIndex)->getIsoAt(isoIndex)->getQ296K();
            sigmaT = par->getMolAt(molIndex)->getIsoAt(isoIndex)->getTipsAtIndex(round(par->getHorizParams()->getTemp()) - par->getTipsTempValueAtIndex(0));  // July 3, 2017 - There was a change in the TIPS files, the TIPS values now start from 1 Kelvin as opposed to 70 in the past

            // Calculating the line intensity taking into account the TIPS values
            lineIntens = ((*recs)[i]).getIntens() * (sigmaTRef / sigmaT) * (ch / zn);

            // Filtering by LineIntensity: compare it with IntencityThreshold. Need to apply wing narrowing instead of filtering, this would be more appropriate
            if (lineIntens < par->getIntensThresh()) { continue; }  // Skip this line

            // Lorentz broadening coefficient
            if (par->getLineBroadening() == 0)  // Air broadening only
                { gamma0 = ((*recs)[i]).getAirHw() * pressRatio * (1 - par->getHorizParams()->getPartPress(molIndex) * 0.000001) * pow((double) tempRatio, ((*recs)[i]).getNAir()); }
            else if (par->getLineBroadening() == 1)  // Self
                { gamma0 = ((*recs)[i]).getSelfHw() * pressRatio * par->getHorizParams()->getPartPress(molIndex) * 0.000001 * pow((double) tempRatio, ((*recs)[i]).getNAir()); }
            else if (par->getLineBroadening() == 2)  // Both air and self
                { gamma0 = pressRatio * ( (((*recs)[i]).getAirHw() * (1 - par->getHorizParams()->getPartPress(molIndex) * 0.000001) + ((*recs)[i]).getSelfHw() * par->getHorizParams()->getPartPress(molIndex) * 0.000001) * pow((double) tempRatio, ((*recs)[i]).getNAir()) ); }

            // Get final wing of the line according to Gamma0, OmegaWingHW and OmegaWing, the maximum of the 2 values is used
            if ((par->getWingContribAbs() > (par->getWingContribSatHw() * gamma0))) { wnWingF = par->getWingContribAbs(); }
            else { wnWingF = par->getWingContribSatHw() * gamma0; }

            // Determining the lower and upper boundary indexes for a single spectral line based on the wing contribution for a single line
            wnShifted = ((*recs)[i]).getWn() + shift0;
            // boundIndexCenter = wavelength->getWns()->count() * ((wnWithShift - wavelength->getStartWn()) / (wavelength->getEndWn() - wavelength->getStartWn()));
            boundIndexLower = wavelength->getWns()->count() * ((wnShifted - wnWingF - wavelength->getStartWn()) / (wavelength->getEndWn() - wavelength->getStartWn()));
            // boundIndexUpper = 2 * boundIndexCenter - boundIndexLower;
            boundIndexUpper = wavelength->getWns()->count() * ((wnShifted + wnWingF - wavelength->getStartWn()) / (wavelength->getEndWn() - wavelength->getStartWn()));
            spec->setStartIndex(boundIndexLower);
            spec->setEndIndex(boundIndexUpper);
            spec->getY()->resize(boundIndexUpper - boundIndexLower + 1);

            // Calculating Lorentz lineshape profile, note from HAPI
            // def PROFILE_LORENTZ(sg0,Gam0,sg):
            //
            // Input parameters:
            // sg0: Unperturbed line position in cm-1 (Input).
            // Gam0: Speed-averaged line-width in cm-1 (Input).
            // sg: Current WaveNumber of the Computation in cm-1 (Input).
            //
            // return Gam0/(pi*(Gam0**2+(sg-sg0)**2))
            //
            // Optical depth determination:
            // OptDepth = xSect * Const::Loch * (296 / Temp) * ppmv * pathLength

            // Calculating m2 multiplication factor
            m2 = m1 * gamma0 * factor * lineIntens * par->getMolAt(molIndex)->getIsoAt(isoIndex)->getAbund() / (par->getMolAt(molIndex)->getIsoAt(isoIndex)->getTerrAbund());
            m2 = (m2 / par->getHorizParams()->getTemp()) * par->getHorizParams()->getPartPress(molIndex) * par->getHorizParams()->getPathLength() * pressRatio;
            gamma2 = gamma0 * gamma0;

            // Limiting the calculated region to the current spectral range
            if (spec->getStartIndex() < 0) { spec->setStartIndex(0); }
            if (spec->getEndIndex() > (wavelength->getWns()->count() - 1)) { spec->setEndIndex(wavelength->getWns()->count() - 1); }
            // Later take into account symmetry in the lineshape profile to make the calculations 2 times faster (only calculate half of the spectral line profile)
            // Later on will probably change this to make it a separate isotopes addition
            // Maybe it makes sense to later include the sides in the caluclated single spectral lines
            for (j = spec->getStartIndex(); j <= spec->getEndIndex(); j++)
            {
                // Calculating the value for the current wavenumber position
                optDepth =  m2 / (gamma2 + (wavelength->getWnsAt(j) - wnShifted) * (wavelength->getWnsAt(j) - wnShifted));
                // spec->setYAt(j - spec->getStartIndex(), optDepth);

                // Appending the spectrum to the total molecules spectrum
                optDepth = optDepth + horizSpectra->getMolAt(molIndex)->getSpec()->getYAt(j);
                horizSpectra->getMolAt(molIndex)->getSpec()->setYAt(j, optDepth);
            }

            // Alternative implementation where spectral lines at the ends of the spectral regions are used completely
            //        if (oneLineSuccess)
            //        {
            //            if (spec->getStartIndex() < 0) { startIndex = 0; }
            //            if (spec->getEndIndex() > (wavelength->getWns()->count() - 1)) { endIndex = wavelength->getWns()->count() - 1; }
            //            // Appending the spectrum to the total spectrum for a given molecule
            //            // (the composite is calculated later on by adding all molecular spectra together)
            //            for (j = startIndex; j <= endIndex; j++)
            //            {
            //                molIndex = (*recs)[i].getMolNum() - 1;
            //                value = (*layerSpectra).getMolAt(molIndex)->getSpec()->getYAt(j) + spec->getYAt(j - spec->getStartIndex());
            //                (*layerSpectra).getMolAt(molIndex)->getSpec()->setYAt(j, value);
            //            }
            //            // Separate isotopes and separate lines spectra are not in the picture right now but will be added later on
            //        }

        }
    }

    // Updating the plot
    findMinMaxValuesInSpectra();     // Finding the maximum spectra
    emit updateSpectra(true, true);  // Emitting a signal to update the plotted spectra
    QApplication::processEvents();   // Updating the processed plot

    // If the cancel button has been pressed, then the calculation is canceled
    if (this->cancelCalculation == true)
        { postCalculation(); }
    else // Requesting the next segment
        { hitran->getNextSubrange(); }

    // Updating the processed plot
    QApplication::processEvents();
}

// Slot to calculate on e Lorentz wavenumber segment
void Lbl::calcDopplerWnSegment()
{
    int i;  // Loop variable
    int j;  // Loop variable
    int molIndex;      // Molecule index
    int isoIndex;      // Isotope index
    int totalNumRecs;  // Total number of records, used to display the calculation progress
    QVector<HitranOneRec> *recs;  // Array with the read in HITRAN data
    double pressRatio; // Pressure ratio (current to reference of 1 atm) used to save a precalculated value
    Spectrum *spec;    // Spectrum to contain the calculated result for a given record in the HITRAN database
    double optDepth;   // Optical depth value
    double m1;         // Multiplication factor, used to save a pre-calculated value
    double m2;         // Multiplication constant used on an intermediate stage of calculations
    double factor;     // Volume concentration multiplication factor
    double shift0;     // Line center pressure shift coefficent
    double wnWingF;    // Lineshape wing contribution parameter
    double ch;         // Line intensity recalculation exponent containing  nominator
    double zn;         // Line intensity recalculation exponent containing  denominator
    double sigmaT;     // TIPS at the current temperature
    double sigmaTRef;  // TIPS at the reference temperature
    double gammaD;     // GammaD, Doppler broadening coefficient
    int boundIndexLower;  // Lower bounding index
    int boundIndexUpper;  // Upper bounding index
    double lineIntens;    // Calculated line intensity
    double wnShifted;     // Wavenumber with the shift applied
    double sqrtMass;      // Square root of the molar mass



    // Initializng an object to contain the spectral calculation results
    spec = new Spectrum();

    // Calculating the total number of records and the number of records after which the plot is to be updated
    totalNumRecs = hitran->getTotalNumRecsToProcess();

    // Reading the next set of indices
    if (hitran->getHitranFileUsage() == 2) currRecNum = 0;

    // Getting a pointer to the list of records read in
    recs = hitran->getHitranRecs();

    // Calculating the temperature and pressure ratio
    pressRatio = par->getHorizParams()->getTotalPress() / Const::totalPressRef;

    // Calculating the multiplication constant
    m1 = 100 * 0.000001 * 296.0 * Const::Loch;  // 100 is to convert to centimeters and 0.000001 is to convert from parts per million (ppmv), embedded into the m1

    // This was slowing dow the calculation considerably on the desktop, mobile doesn't seem to matter
    // Incrementing the sub-segment HITRAN record number
    QApplication::processEvents();
    emit notify("Found " + QString::number(totalNumRecs) + " HITRAN records.\nRunning the line-by-line calculation.\nProcessing segment " + QString::number(currRecNum + 1) + "-" + QString::number(recs->count() + currRecNum + 1) + ".", 0, true, true);
    QApplication::processEvents();

    // Running calculations for all available spectral lines
    for (i = 0; i < recs->count(); i++)
    {
        // Incrementing the current record number
        currRecNum++;

        // Saving the molecule and isotope indices
        molIndex = ((*recs)[i]).getMolNum() - 1;
        isoIndex = ((*recs)[i]).getIsoNum() - 1;

        // Calculating for selected molecules only unless all molecules are to be included in calculations
        if ((par->getMolAt(molIndex)->getSelected() == true) || par->getCalcNotSelectedMols())
        {
            if ( ((*recs)[i]).getIsoNum() == 0 ) { isoIndex = 10; }  // Special case of the CO2 molecule

            // If the abundance is unavailable, then skipping the current calculation
            if (par->getMolAt(molIndex)->getIsoAt(isoIndex)->getTerrAbund() == 0) { continue; }

            // Precalculation of volume concentration
            // Depending on the HITRAN database units taking into account the multiplication factor
            if (par->getHitranUnits() == true)
                { factor = 1.0; } // Multiplication factor is equal to 1.0
            else // Multiplication factor calculated from the pressure and temperature
                { factor = par->getHorizParams()->getTotalPress() / (9.869233e-7 * Const::BoltsCgs * par->getHorizParams()->getTemp()); }   // CGS units used in this formula

            // Pressure shift, taking parameters from the HITRAN database if the line shift is used
            if (par->getUseLineShift()) { shift0 = ((*recs)[i]).getAirPShift() * pressRatio; }
            else { shift0 = 0; }

            // Partition functions for T and Tref, need to be optimized
            sigmaT = par->getMolAt(molIndex)->getIsoAt(isoIndex)->getTips()->at(round(par->getHorizParams()->getTemp()) - par->getTipsTempValueAtIndex(0));
            sigmaTRef = par->getMolAt(molIndex)->getIsoAt(isoIndex)->getQ296K();

            // Getting all environment dependences from Voigt parameters, line intensity, used to be a separate procedure environmentDependencyIntensity
            ch = exp(-Const::c2 * ((*recs)[i]).getLowE() / par->getHorizParams()->getTemp()) * (1 - exp(-Const::c2 * ((*recs)[i]).getWn() / par->getHorizParams()->getTemp()));
            zn = exp(-Const::c2 * ((*recs)[i]).getLowE() / Const::tempRef) * (1 - exp(-Const::c2 * ((*recs)[i]).getWn() / Const::tempRef));
            lineIntens = ((*recs)[i]).getIntens() * (sigmaTRef / sigmaT) * (ch / zn);

            // Filtering by LineIntensity: compare it with IntencityThreshold. Need to apply wing narrowing instead of filtering, this would be more appropriate
            if (lineIntens < par->getIntensThresh()) { continue; }  // Skip this line then, not included

            // GammaD calculation
            sqrtMass = sqrt(par->getMolAt(molIndex)->getIsoAt(isoIndex)->getMolarMass());
            gammaD = (Const::sqrt2Ln2 / Const::c) * sqrt(Const::Bolts / Const::cMassMol) * sqrt(par->getHorizParams()->getTemp()) * ((*recs)[i]).getWn() / sqrtMass;

            // Get final wing of the line according to Gamma0, OmegaWingHW and OmegaWing, the maximum of the 2 values is used
            if ((par->getWingContribAbs() > (par->getWingContribSatHw() * gammaD))) { wnWingF = par->getWingContribAbs(); }
            else { wnWingF = par->getWingContribSatHw() * gammaD; }

            // Determining the lower and upper boundary indexes for a single spectral line based on the wing contribution for a single line
            wnShifted = ((*recs)[i]).getWn() + shift0;
            // boundIndexCenter = wavelength->getWns()->count() * ((wnWithShift - wavelength->getStartWn()) / (wavelength->getEndWn() - wavelength->getStartWn()));
            boundIndexLower = wavelength->getWns()->count() * ((wnShifted - wnWingF - wavelength->getStartWn()) / (wavelength->getEndWn() - wavelength->getStartWn()));
            // boundIndexUpper = 2 * boundIndexCenter - boundIndexLower;
            boundIndexUpper = wavelength->getWns()->count() * ((wnShifted + wnWingF - wavelength->getStartWn()) / (wavelength->getEndWn() - wavelength->getStartWn()));
            spec->setStartIndex(boundIndexLower);
            spec->setEndIndex(boundIndexUpper);
            spec->getY()->resize(boundIndexUpper - boundIndexLower + 1);

            // Doppler profile
            //
            // Input parameters:
            // sg0: Unperturbed line position in cm-1 (Input).
            // GamD: Doppler HWHM in cm-1 (Input)
            // sg: Current WaveNumber of the Computation in cm-1 (Input).
            //
            // return cSqrtLn2divSqrtPi * exp(-cLn2*((sg-sg0) / GamD)**2) / GamD
            // Optical depth determination:
            // OptDepth = xSect * Const::Loch * (296 / Temp) * ppmv * pathLength

            // Calculating m2 multiplication factor
            m2 = (m1 / gammaD) * factor * lineIntens * par->getMolAt(molIndex)->getIsoAt(isoIndex)->getAbund() / (par->getMolAt(molIndex)->getIsoAt(isoIndex)->getTerrAbund());
            m2 = (m2 / par->getHorizParams()->getTemp()) * par->getHorizParams()->getPartPress(molIndex) * par->getHorizParams()->getPathLength() * pressRatio;

            // Limiting the calculated region to the current spectral range
            if (spec->getStartIndex() < 0) { spec->setStartIndex(0); }
            if (spec->getEndIndex() > (wavelength->getWns()->count() - 1)) { spec->setEndIndex(wavelength->getWns()->count() - 1); }
            // Later take into account symmetry in the lineshape profile in the expectation to make the calculations 2 times faster (only calculate half of the spectral line profile)
            // Later on will probably change this to make it a separate isotopes addition
            // Maybe it makes sense to later include the sides in the caluclated single spectral lines
            for (j = spec->getStartIndex(); j <= spec->getEndIndex(); j++)
            {
                // Calculating the value for the current wavenumber position
                // cSqrtLn2divSqrtPi * exp(-cLn2*((sg-sg0) / GamD)**2)
                optDepth = m2 * Const::sqrtLn2divSqrtPi * exp(-Const::ln2 * ((wavelength->getWnsAt(j) - wnShifted) / gammaD) * ((wavelength->getWnsAt(j) - wnShifted) / gammaD));
                // spec->setYAt(j - spec->getStartIndex(), optDepth);

                // Appending the spectrum to the total molecules spectrum
                optDepth = optDepth + horizSpectra->getMolAt(molIndex)->getSpec()->getYAt(j);
                horizSpectra->getMolAt(molIndex)->getSpec()->setYAt(j, optDepth);
            }

            // Alternative implementation where spectral lines at the ends of the spectral regions are used completely
            //        if (oneLineSuccess)
            //        {
            //            if (spec->getStartIndex() < 0) { startIndex = 0; }
            //            if (spec->getEndIndex() > (wavelength->getWns()->count() - 1)) { endIndex = wavelength->getWns()->count() - 1; }
            //            // Appending the spectrum to the total spectrum for a given molecule
            //            // (the composite is calculated later on by adding all molecular spectra together)
            //            for (j = startIndex; j <= endIndex; j++)
            //            {
            //                molIndex = (*recs)[i].getMolNum() - 1;
            //                value = (*layerSpectra).getMolAt(molIndex)->getSpec()->getYAt(j) + spec->getYAt(j - spec->getStartIndex());
            //                (*layerSpectra).getMolAt(molIndex)->getSpec()->setYAt(j, value);
            //            }
            //            // Separate isotopes and separate lines spectra are not in the picture right now but will be added later on
            //        }
        }
    }

    // Updating the plot
    findMinMaxValuesInSpectra();     // Finding the maximum spectra
    emit updateSpectra(true, true);  // Emitting a signal to update the plotted spectra
    QApplication::processEvents();   // Updating the processed plot

    // If the cancel button has been pressed, then the calculation is canceled
    if (this->cancelCalculation == true)
        { postCalculation(); }
    else // Requesting the next segment
        { hitran->getNextSubrange(); }

    // Updating the processed plot
    QApplication::processEvents();
}

// The following routine relies on the Voigt profile calculation algorithm and code described in the following publication:
// Abrarov, S., & Quine, B. (2015). A Rational Approximation for Efficient Computation of the Voigt Function
// in Quantitative Spectroscopy. Journal Of Mathematics Research, 7(2), p163. doi:http://dx.doi.org/10.5539/jmr.v7n2
// Slot to calculate one Voigt wavenumber segment (Low resolution implementation)
// Also: https://www.princeton.edu/cefrc/Files/2013%20Lecture%20Notes/Hanson/pLecture6.pdf
void Lbl::calcVoigtWnSegment()
{
    int i;  // Loop variable
    int j;  // Loop variable
    int t;  // Loop variable
    int molIndex;      // Molecule index
    int isoIndex;      // Isotope index
    int totalNumRecs;  // Total number of records, used to display the calculation progress
    QVector<HitranOneRec> *recs;  // Array with the read in HITRAN data
    double pressRatio; // Pressure ratio (current to reference of 1 atm) used to save a precalculated value
    double tempRatio;  // Temperature ratio
    Spectrum *spec;    // Spectrum to contain the calculated result for a given record in the HITRAN database
    double optDepth;   // Optical depth value
    double m1;         // Multiplication factor, used to save a pre-calculated value
    double m2;         // Multiplication constant used on an intermediate stage of calculations
    double factor;     // Volume concentration multiplication factor
    double shift0;     // Line center pressure shift coefficent
    double wnWingF;    // Lineshape wing contribution parameter
    double ch;         // Line intensity recalculation exponent containing  nominator
    double zn;         // Line intensity recalculation exponent containing  denominator
    double sigmaT;     // TIPS at the current temperature
    double sigmaTRef;  // TIPS at the reference temperature
    double gammaD;     // GammaD, Doppler broadening coefficient
    double gamma0 = 1;    // Gamma 0, Lorentz broadening coefficient
    int boundIndexLower;  // Lower bounding index
    int boundIndexUpper;  // Upper bounding index
    double lineIntens;    // Calculated line intensity
    double wnShifted;     // Wavenumber with the shift applied
    double sqrtMass;      // Square root of the molar mass

    // Auxiliaty variables for the calculation of the Voigt profile
    double mVoigt; // Multiplication term for the Voigt profile
    double vf;     // Value of the voigt function
    double xTerm;
    double x;
    double y;
    double varSigma = 2.75;
    double arr1;
    double arr2;
    double arr3;
    double ySquared;



    // Initializng an object to contain the spectral calculation results
    spec = new Spectrum();

    // Calculating the total number of records and the number of records after which the plot is to be updated
    totalNumRecs = hitran->getTotalNumRecsToProcess();

    // Reading the next set of indices
    if (hitran->getHitranFileUsage() == 2) currRecNum = 0;

    // Getting a pointer to the list of records read in
    recs = hitran->getHitranRecs();

    // Calculating the temperature and pressure ratio
    pressRatio = par->getHorizParams()->getTotalPress() / Const::totalPressRef;
    tempRatio = Const::tempRef / par->getHorizParams()->getTemp();

    // Calculating the multiplication constant
    m1 = 100 * 0.000001 * 296.0 * Const::Loch;  // Path is converted to cm, 0.000001 is to convert from parts per million

    // This was slowing dow the calculation considerably on the desktop, mobile doesn't seem to matter
    // Incrementing the sub-segment HITRAN record number
    QApplication::processEvents();
    emit notify("Found " + QString::number(totalNumRecs) + " HITRAN records.\nRunning the line-by-line calculation.\nProcessing segment " + QString::number(currRecNum + 1) + "-" + QString::number(recs->count() + currRecNum + 1) + ".", 0, true, true);
    QApplication::processEvents();

    // Running calculations for all available spectral lines
    for (i = 0; i < recs->count(); i++)
    {
        // Incrementing the current record number
        currRecNum++;

        // Saving the molecule and isotope indices
        molIndex = ((*recs)[i]).getMolNum() - 1;
        isoIndex = ((*recs)[i]).getIsoNum() - 1;

        // Calculating for selected molecules only unless all molecules are to be included in calculations
        if ((par->getMolAt(molIndex)->getSelected() == true) || par->getCalcNotSelectedMols())
        {
            if ( ((*recs)[i]).getIsoNum() == 0 ) { isoIndex = 10; }  // Special case of the CO2 molecule

            // If the abundance is unavailable, then skipping the current calculation
            if (par->getMolAt(molIndex)->getIsoAt(isoIndex)->getTerrAbund() == 0) { continue; }

            // Precalculation of volume concentration
            // Depending on the HITRAN database units taking into account the multiplication factor
            if (par->getHitranUnits() == true) // Multiplication factor is equal to 1.0
                { factor = 1.0; }
            else // Multiplication factor calculated from the pressure and temperature
                { factor = par->getHorizParams()->getTotalPress() / (9.869233e-7 * Const::BoltsCgs * par->getHorizParams()->getTemp()); }  // multiplication factor is in CGS units

            // Pressure shift, taking parameters from the HITRAN database if the line shift is used
            if (par->getUseLineShift()) { shift0 = ((*recs)[i]).getAirPShift() * pressRatio; }
            else { shift0 = 0; }

            // Partition functions for T and Tref, need to be optimized
            sigmaT = par->getMolAt(molIndex)->getIsoAt(isoIndex)->getTips()->at(round(par->getHorizParams()->getTemp()) - par->getTipsTempValueAtIndex(0));
            sigmaTRef = par->getMolAt(molIndex)->getIsoAt(isoIndex)->getQ296K();

            // Getting all environment dependences from Voigt parameters, line intensity, used to be a separate procedure environmentDependencyIntensity
            ch = exp(-Const::c2 * ((*recs)[i]).getLowE() / par->getHorizParams()->getTemp()) * (1 - exp(-Const::c2 * ((*recs)[i]).getWn() / par->getHorizParams()->getTemp()));
            zn = exp(-Const::c2 * ((*recs)[i]).getLowE() / Const::tempRef) * (1 - exp(-Const::c2 * ((*recs)[i]).getWn() / Const::tempRef));
            lineIntens = ((*recs)[i]).getIntens() * (sigmaTRef / sigmaT) * (ch / zn);

            // Filtering by LineIntensity: compare it with IntensityThreshold. Need to apply wing narrowing instead of filtering, this would be more appropriate
            if (lineIntens < par->getIntensThresh()) { continue; }  // Skip this line

            // GammaD calculation
            sqrtMass = sqrt(par->getMolAt(molIndex)->getIsoAt(isoIndex)->getMolarMass());
            gammaD = (Const::sqrt2Ln2 / Const::c) * sqrt(Const::Bolts / Const::cMassMol) * sqrt(par->getHorizParams()->getTemp()) * ((*recs)[i]).getWn() / sqrtMass;

            // Lorentz broadening coefficient, inlined method "environmentDependencyGamma0"
            // a^x -> exp(x * ln(a)) = a^x
            // Setting the broadening
            if (par->getLineBroadening() == 0)  // Air broadening only
                { gamma0 = ((*recs)[i]).getAirHw() * pressRatio * (1 - par->getHorizParams()->getPartPress(molIndex) * 0.000001) * pow((double) tempRatio, ((*recs)[i]).getNAir()); }
            else if (par->getLineBroadening() == 1)  // Self
                { gamma0 = ((*recs)[i]).getSelfHw() * pressRatio * par->getHorizParams()->getPartPress(molIndex) * 0.000001 * pow((double) tempRatio, ((*recs)[i]).getNAir()); }
            else if (par->getLineBroadening() == 2)  // Both air and self
                { gamma0 = pressRatio * ( (((*recs)[i]).getAirHw() * (1 - par->getHorizParams()->getPartPress(molIndex) * 0.000001) + ((*recs)[i]).getSelfHw() * par->getHorizParams()->getPartPress(molIndex) * 0.000001) * pow((double) tempRatio, ((*recs)[i]).getNAir()) ); }

            // Get final wing of the line according to Gamma0, OmegaWingHW and OmegaWing, the maximum of the 2 values is used
            // OmegaWingF = max(OmegaWing, OmegaWingHW * Gamma0, OmegaWingHW * GammaD);
            if ((par->getWingContribAbs() > (par->getWingContribSatHw() * gamma0)) && (par->getWingContribAbs() > (par->getWingContribSatHw() * gammaD))) { wnWingF = par->getWingContribAbs();; }
            if (((par->getWingContribSatHw() * gamma0) > par->getWingContribAbs()) && (gamma0 > gammaD)) { wnWingF = (par->getWingContribSatHw() * gamma0); }
            if (((par->getWingContribSatHw() * gammaD) > par->getWingContribAbs()) && (gammaD > gamma0)) { wnWingF = (par->getWingContribSatHw() * gammaD); }

            // Determining the lower and upper boundary indexes for a single spectral line based on the wing contribution for a single line
            wnShifted = ((*recs)[i]).getWn() + shift0;
            // boundIndexCenter = wavelength->getWns()->count() * ((wnWithShift - wavelength->getStartWn()) / (wavelength->getEndWn() - wavelength->getStartWn()));
            boundIndexLower = wavelength->getWns()->count() * ((wnShifted - wnWingF - wavelength->getStartWn()) / (wavelength->getEndWn() - wavelength->getStartWn()));
            // boundIndexUpper = 2 * boundIndexCenter - boundIndexLower;
            boundIndexUpper = wavelength->getWns()->count() * ((wnShifted + wnWingF - wavelength->getStartWn()) / (wavelength->getEndWn() - wavelength->getStartWn()));
            spec->setStartIndex(boundIndexLower);
            spec->setEndIndex(boundIndexUpper);
            spec->getY()->resize(boundIndexUpper - boundIndexLower + 1);

            // Auxiliary intermediate terms for the calculation for the Voigt profile
            xTerm = Const::sqrtLn2 / gammaD;  // 2 multiplication term is absent because Princeton is FWHM and here it is HWHM
            y = Const::sqrtLn2 * gamma0 / gammaD;
            y = fabs(y) + varSigma * 0.5;
            ySquared = y * y;

            // This term is according to the lecture from Princeton
            // https://www.princeton.edu/cefrc/Files/2013%20Lecture%20Notes/Hanson/pLecture6.pdf
            mVoigt = Const::sqrtLn2 / (sqrt(Const::pi) * gammaD);  // 2 multiplication term is absent because Princeton is FWHM and here it is HWHM

            // Optical depth determination:
            // OptDepth = xSect * Const::Loch * (296 / Temp) * ppmv * pathLength

            // sg0 - unperturbed line position
            //  sg - wavenumber grid (current wavenumbers to run calculations for)
            m2 = m1 * factor * lineIntens * par->getMolAt(molIndex)->getIsoAt(isoIndex)->getAbund() / (par->getMolAt(molIndex)->getIsoAt(isoIndex)->getTerrAbund());
            m2 = (m2 / par->getHorizParams()->getTemp()) * par->getHorizParams()->getPartPress(molIndex) * par->getHorizParams()->getPathLength() * pressRatio;

            // Limiting the calculated region to the current spectral range
            if (spec->getStartIndex() < 0) { spec->setStartIndex(0); }
            if (spec->getEndIndex() > (wavelength->getWns()->count() - 1)) { spec->setEndIndex(wavelength->getWns()->count() - 1); }
            // Later take into account symmetry in the lineshape profile to make the calculations 2 times faster (only calculate half of the spectral line profile)
            // Later on will probably change this to make it a separate isotopes addition
            // Maybe it makes sense to later include the sides in the caluclated single spectral lines
            for (j = spec->getStartIndex(); j < spec->getEndIndex(); j++)
            {
                // The following routine relies on the Voigt profile calculation algorithm and code described in the following publication:
                // Abrarov, S., & Quine, B. (2015). A Rational Approximation for Efficient Computation of the Voigt Function
                // in Quantitative Spectroscopy. Journal Of Mathematics Research, 7(2), p163. doi:http://dx.doi.org/10.5539/jmr.v7n2
                // Slot to calculate one Voigt wavenumber segment (Low resolution implementation)
                // Preparing the terms
                x = xTerm * (wavelength->getWnsAt(j) - wnShifted);
                arr1 = ySquared - x * x;
                arr2 = x * x + ySquared;
                arr3 = arr2 * arr2;
                // Calculating the Voigt profile value
                vf = 0;
                for (t = 0; t < coeff1_12.count(); t++)
                {
                    vf = vf + (coeff1_12[t] * (coeff2_12[t] + arr1) + coeff3_12[t] * y * (coeff2_12[t] + arr2)) / (coeff2_12[t] * coeff2_12[t] + 2 * coeff2_12[t] * arr1 + arr3) ;
                }
                optDepth = m2 * mVoigt * vf;
                optDepth = optDepth + horizSpectra->getMolAt(molIndex)->getSpec()->getYAt(j);
                horizSpectra->getMolAt(molIndex)->getSpec()->setYAt(j, optDepth);

                // Calculating the value for the current wavenumber position
                // spec->setYAt(j - spec->getStartIndex(), optDepth);
            }

            // Alternative implementation where spectral lines at the ends of the spectral regions are used completely
            //        if (oneLineSuccess)
            //        {
            //            if (spec->getStartIndex() < 0) { startIndex = 0; }
            //            if (spec->getEndIndex() > (wavelength->getWns()->count() - 1)) { endIndex = wavelength->getWns()->count() - 1; }
            //            // Appending the spectrum to the total spectrum for a given molecule
            //            // (the composite is calculated later on by adding all molecular spectra together)
            //            for (j = startIndex; j <= endIndex; j++)
            //            {
            //                molIndex = (*recs)[i].getMolNum() - 1;
            //                value = (*layerSpectra).getMolAt(molIndex)->getSpec()->getYAt(j) + spec->getYAt(j - spec->getStartIndex());
            //                (*layerSpectra).getMolAt(molIndex)->getSpec()->setYAt(j, value);
            //            }
            //            // Separate isotopes and separate lines spectra are not in the picture right now but will be added later on
            //        }
        }
    }

    // Updating the plot
    findMinMaxValuesInSpectra();     // Finding the maximum spectra
    emit updateSpectra(true, true);  // Emitting a signal to update the plotted spectra
    QApplication::processEvents();   // Updating the processed plot

    // If the cancel button has been pressed, then the calculation is canceled
    if (this->cancelCalculation == true)
        { postCalculation(); }
    else // Requesting the next segment
        { hitran->getNextSubrange(); }

    // Updating the processed plot
    QApplication::processEvents();
}

// The following routine relies on the Voigt profile calculation algorithm and code described in the following publication:
// Abrarov, S., & Quine, B. (2015). A Rational Approximation for Efficient Computation of the Voigt Function
// in Quantitative Spectroscopy. Journal Of Mathematics Research, 7(2), p163. doi:http://dx.doi.org/10.5539/jmr.v7n2
// This High resolution implementration relies on a 16-coefficient implementation
void Lbl::calcVoigtWnSegmentHighRes()
{
    int i;  // Loop variable
    int j;  // Loop variable
    int t;  // Loop variable
    int molIndex;      // Molecule index
    int isoIndex;      // Isotope index
    int totalNumRecs;  // Total number of records, used to display the calculation progress
    QVector<HitranOneRec> *recs;  // Array with the read in HITRAN data
    double pressRatio; // Pressure ratio (current to reference of 1 atm) used to save a precalculated value
    double tempRatio;  // Temperature ratio
    Spectrum *spec;    // Spectrum to contain the calculated result for a given record in the HITRAN database
    double optDepth;   // Optical depth value
    double m1;         // Multiplication factor, used to save a pre-calculated value
    double m2;         // Multiplication constant used on an intermediate stage of calculations
    double factor;     // Volume concentration multiplication factor
    double shift0;     // Line center pressure shift coefficent
    double wnWingF;    // Lineshape wing contribution parameter
    double ch;         // Line intensity recalculation exponent containing  nominator
    double zn;         // Line intensity recalculation exponent containing  denominator
    double sigmaT;     // TIPS at the current temperature
    double sigmaTRef;  // TIPS at the reference temperature
    double gammaD;     // GammaD, Doppler broadening coefficient
    double gamma0 = 1;    // Gamma 0, Lorentz broadening coefficient
    int boundIndexLower;  // Lower bounding index
    int boundIndexUpper;  // Upper bounding index
    double lineIntens;    // Calculated line intensity
    double wnShifted;     // Wavenumber with the shift applied
    double sqrtMass;      // Square root of the molar mass

    // Auxiliaty variables for the calculation of the Voigt profile
    double mVoigt; // Multiplication term for the Voigt profile
    double vf;     // Value of the voigt function
    double xTerm;
    double x;
    double y;
    double varSigma = 2.75;
    double arr1;
    double arr2;
    double arr3;
    double ySquared;



    // Initializng an object to contain the spectral calculation results
    spec = new Spectrum();

    // Calculating the total number of records and the number of records after which the plot is to be updated
    totalNumRecs = hitran->getTotalNumRecsToProcess();

    // Reading the next set of indices
    if (hitran->getHitranFileUsage() == 2) currRecNum = 0;

    // Getting a pointer to the list of records read in
    recs = hitran->getHitranRecs();

    // Calculating the temperature and pressure ratio
    pressRatio = par->getHorizParams()->getTotalPress() / Const::totalPressRef;
    tempRatio = Const::tempRef / par->getHorizParams()->getTemp();

    // Calculating the multiplication constant
    m1 = 100 * 0.000001 * 296.0 * Const::Loch;  // 100 is to convert to centimeters and 0.000001 is to convert from parts per million (ppmv)

    // This was slowing dow the calculation considerably on the desktop, mobile doesn't seem to matter
    // Incrementing the sub-segment HITRAN record number
    QApplication::processEvents();
    emit notify("Found " + QString::number(totalNumRecs) + " HITRAN records.\nRunning the line-by-line calculation.\nProcessing segment " + QString::number(currRecNum + 1) + "-" + QString::number(recs->count() + currRecNum + 1) + ".", 0, true, true);
    QApplication::processEvents();

    // Running calculations for all available spectral lines
    for (i = 0; i < recs->count(); i++)
    {
        // Incrementing the current record number
        currRecNum++;

        // Saving the molecule and isotope indices
        molIndex = ((*recs)[i]).getMolNum() - 1;
        isoIndex = ((*recs)[i]).getIsoNum() - 1;

        // Calculating for selected molecules only unless all molecules are to be included in calculations
        if ((par->getMolAt(molIndex)->getSelected() == true) || par->getCalcNotSelectedMols())
        {
            if ( ((*recs)[i]).getIsoNum() == 0 ) { isoIndex = 10; }  // Special case of the CO2 molecule

            // If the abundance is unavailable, then skipping the current calculation
            if (par->getMolAt(molIndex)->getIsoAt(isoIndex)->getTerrAbund() == 0) { continue; }

            // Precalculation of volume concentration
            // Depending on the HITRAN database units taking into account the multiplication factor
            if (par->getHitranUnits() == true)
                { factor = 1.0; }  // Multiplication factor is equal to 1.0
            else // Multiplication factor calculated from the pressure and temperature
                { factor = par->getHorizParams()->getTotalPress() / (9.869233e-7 * Const::BoltsCgs * par->getHorizParams()->getTemp()); }  // This multiplication factor is in CGS units

            // Pressure shift, taking parameters from the HITRAN database if the line shift is used
            if (par->getUseLineShift()) { shift0 = ((*recs)[i]).getAirPShift() * pressRatio; }
            else { shift0 = 0; }

            // Partition functions for T and Tref, need to be optimized
            sigmaT = par->getMolAt(molIndex)->getIsoAt(isoIndex)->getTips()->at(round(par->getHorizParams()->getTemp()) - par->getTipsTempValueAtIndex(0));   // partitionFunction(MoleculeNumberDB,IsoNumberDB,T)
            sigmaTRef = par->getMolAt(molIndex)->getIsoAt(isoIndex)->getQ296K();   // partitionFunction(MoleculeNumberDB, IsoNumberDB, Tref)

            // Getting all environment dependences from Voigt parameters, line intensity, used to be a separate procedure environmentDependencyIntensity
            ch = exp(-Const::c2 * ((*recs)[i]).getLowE() / par->getHorizParams()->getTemp()) * (1 - exp(-Const::c2 * ((*recs)[i]).getWn() / par->getHorizParams()->getTemp()));
            zn = exp(-Const::c2 * ((*recs)[i]).getLowE() / Const::tempRef) * (1 - exp(-Const::c2 * ((*recs)[i]).getWn() / Const::tempRef));
            lineIntens = ((*recs)[i]).getIntens() * (sigmaTRef / sigmaT) * (ch / zn);

            // Filtering by LineIntensity: compare it with IntensityThreshold. Need to apply wing narrowing instead of filtering, this would be more appropriate
            if (lineIntens < par->getIntensThresh()) { continue; }  // Skip this line then, not included

            // GammaD calculation
            sqrtMass = sqrt(par->getMolAt(molIndex)->getIsoAt(isoIndex)->getMolarMass());
            gammaD = (Const::sqrt2Ln2 / Const::c) * sqrt(Const::Bolts / Const::cMassMol) * sqrt(par->getHorizParams()->getTemp()) * ((*recs)[i]).getWn() / sqrtMass;

            // Lorentz broadening coefficient, inlined method "environmentDependencyGamma0"
            // a^x -> exp(x * ln(a)) = a^x
            // Setting the broadening
            if (par->getLineBroadening() == 0)  // Air broadening only
                { gamma0 = ((*recs)[i]).getAirHw() * pressRatio * (1 - par->getHorizParams()->getPartPress(molIndex) * 0.000001) * pow((double) tempRatio, ((*recs)[i]).getNAir()); }
            else if (par->getLineBroadening() == 1)  // Self
                { gamma0 = ((*recs)[i]).getSelfHw() * pressRatio * par->getHorizParams()->getPartPress(molIndex) * 0.000001 * pow((double) tempRatio, ((*recs)[i]).getNAir()); }
            else if (par->getLineBroadening() == 2)  // Both air and self
                { gamma0 = pressRatio * ( (((*recs)[i]).getAirHw() * (1 - par->getHorizParams()->getPartPress(molIndex) * 0.000001) + ((*recs)[i]).getSelfHw() * par->getHorizParams()->getPartPress(molIndex) * 0.000001) * pow((double) tempRatio, ((*recs)[i]).getNAir()) ); }

            // Get final wing of the line according to Gamma0, OmegaWingHW and OmegaWing, the maximum of the 2 values is used
            // OmegaWingF = max(OmegaWing, OmegaWingHW * Gamma0, OmegaWingHW * GammaD);
            if ((par->getWingContribAbs() > (par->getWingContribSatHw() * gamma0)) && (par->getWingContribAbs() > (par->getWingContribSatHw() * gammaD))) { wnWingF = par->getWingContribAbs();; }
            if (((par->getWingContribSatHw() * gamma0) > par->getWingContribAbs()) && (gamma0 > gammaD)) { wnWingF = (par->getWingContribSatHw() * gamma0); }
            if (((par->getWingContribSatHw() * gammaD) > par->getWingContribAbs()) && (gammaD > gamma0)) { wnWingF = (par->getWingContribSatHw() * gammaD); }

            // Determining the lower and upper boundary indexes for a single spectral line based on the wing contribution for a single line
            wnShifted = ((*recs)[i]).getWn() + shift0;
            // boundIndexCenter = wavelength->getWns()->count() * ((wnWithShift - wavelength->getStartWn()) / (wavelength->getEndWn() - wavelength->getStartWn()));
            boundIndexLower = wavelength->getWns()->count() * ((wnShifted - wnWingF - wavelength->getStartWn()) / (wavelength->getEndWn() - wavelength->getStartWn()));
            // boundIndexUpper = 2 * boundIndexCenter - boundIndexLower;
            boundIndexUpper = wavelength->getWns()->count() * ((wnShifted + wnWingF - wavelength->getStartWn()) / (wavelength->getEndWn() - wavelength->getStartWn()));
            spec->setStartIndex(boundIndexLower);
            spec->setEndIndex(boundIndexUpper);
            spec->getY()->resize(boundIndexUpper - boundIndexLower + 1);

            // Auxiliary intermediate terms for the calculation for the Voigt profile
            xTerm = Const::sqrtLn2 / gammaD;  // 2 multiplication term is absent because Princeton is FWHM and here it is HWHM
            y = Const::sqrtLn2 * gamma0 / gammaD;
            y = fabs(y) + varSigma * 0.5;
            ySquared = y * y;

            // This term is according to the lecture from Princeton
            // https://www.princeton.edu/cefrc/Files/2013%20Lecture%20Notes/Hanson/pLecture6.pdf
            mVoigt = Const::sqrtLn2 / (sqrt(Const::pi) * gammaD);  // 2 multiplication term is absent because Princeton is FWHM and here it is HWHM

            // Optical depth determination:
            // OptDepth = xSect * Const::Loch * (296 / Temp) * ppmv * pathLength

            // Calculating an intermediate multiplication factor
            m2 = m1 * factor * lineIntens * par->getMolAt(molIndex)->getIsoAt(isoIndex)->getAbund() / (par->getMolAt(molIndex)->getIsoAt(isoIndex)->getTerrAbund());
            m2 = (m2 / par->getHorizParams()->getTemp()) * par->getHorizParams()->getPartPress(molIndex) * par->getHorizParams()->getPathLength() * pressRatio;

            // Limiting the calculated region to the current spectral range
            if (spec->getStartIndex() < 0) { spec->setStartIndex(0); }
            if (spec->getEndIndex() > (wavelength->getWns()->count() - 1)) { spec->setEndIndex(wavelength->getWns()->count() - 1); }
            // Later take into account symmetry in the lineshape profile to make the calculations 2 times faster (only calculate half of the spectral line profile)
            // Later on will probably change this to make it a separate isotopes addition
            // Maybe it makes sense to later include the sides in the caluclated single spectral lines
            for (j = spec->getStartIndex(); j <= spec->getEndIndex(); j++)
            {
                // The following routine relies on the Voigt profile calculation algorithm and code described in the following publication:
                // Abrarov, S., & Quine, B. (2015). A Rational Approximation for Efficient Computation of the Voigt Function
                // in Quantitative Spectroscopy. Journal Of Mathematics Research, 7(2), p163. doi:http://dx.doi.org/10.5539/jmr.v7n2
                // Slot to calculate one Voigt wavenumber segment (Low resolution implementation)
                // Preparing the terms
                x = xTerm * (wavelength->getWnsAt(j) - wnShifted);
                arr1 = ySquared - x * x;
                arr2 = x * x + ySquared;
                arr3 = arr2 * arr2;
                // Calculating the Voigt profile value
                vf = 0;
                for (t = 0; t < coeff1_16.count(); t++)
                {
                    vf = vf + (coeff1_16[t] * (coeff2_16[t] + arr1) + coeff3_16[t] * y * (coeff2_16[t] + arr2)) / (coeff2_16[t] * coeff2_16[t] + 2 * coeff2_16[t] * arr1 + arr3) ;
                }
                optDepth = m2 * mVoigt * vf;

                // Calculating the value for the current wavenumber position
                // spec->setYAt(j - spec->getStartIndex(), optDepth);

                // Appending the spectrum to the total molecules spectrum
                optDepth = optDepth + horizSpectra->getMolAt(molIndex)->getSpec()->getYAt(j);
                horizSpectra->getMolAt(molIndex)->getSpec()->setYAt(j, optDepth);
            }

            // Alternative implementation where spectral lines at the ends of the spectral regions are used completely
            //        if (oneLineSuccess)
            //        {
            //            if (spec->getStartIndex() < 0) { startIndex = 0; }
            //            if (spec->getEndIndex() > (wavelength->getWns()->count() - 1)) { endIndex = wavelength->getWns()->count() - 1; }
            //            // Appending the spectrum to the total spectrum for a given molecule
            //            // (the composite is calculated later on by adding all molecular spectra together)
            //            for (j = startIndex; j <= endIndex; j++)
            //            {
            //                molIndex = (*recs)[i].getMolNum() - 1;
            //                value = (*layerSpectra).getMolAt(molIndex)->getSpec()->getYAt(j) + spec->getYAt(j - spec->getStartIndex());
            //                (*layerSpectra).getMolAt(molIndex)->getSpec()->setYAt(j, value);
            //            }
            //            // Separate isotopes and separate lines spectra are not in the picture right now but will be added later on
            //        }
        }
    }

    // Updating the plot
    findMinMaxValuesInSpectra();     // Finding the maximum spectra
    emit updateSpectra(true, true);  // Emitting a signal to update the plotted spectra
    QApplication::processEvents();   // Updating the processed plot

    // If the cancel button has been pressed, then the calculation is canceled
    if (this->cancelCalculation == true)
    {
        postCalculation();
    }
    else
    {
        // Requesting the next segment
        hitran->getNextSubrange();
    }
    QApplication::processEvents();   // Updating the processed plot
}
