/*
    Spectroscopic Toolkit version 1.93 by Pieter Suurmond, july 12, 2008.

    Copyright (c) 2000-2008 - Pieter Suurmond

    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.

    Any person wishing to distribute modifications to the Software is
    requested to send the modifications to the original developer so that
    they can be incorporated into the canonical version.

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

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

#include "ST_pconv.h"       /* For kST_bin_base */

/* On Ubuntu 8.04: strtold() prototype is not in <stdlib.h>!!??
   long double strtold(const char *nptr, char **endptr);
USING gcc-option "-std=c99" or "gnu99" helps. */


/*----------------------------------------------------------------------------*/
/*
    Reads a single line from the database. See headerfile ST_pconv.h for format.
*/
static int processLine(char*            textBuff,
                       pieterPacked*    PBIN,
                       unsigned long*   checksum)
{
    char    rem;
    int     n = 0;

    while (textBuff[n]) /* Same as n = strlen(textBuff); */
        *checksum += (unsigned long)((unsigned char)textBuff[n++]);

    if ((n < 155) || (n > 161))     /* Isotope shifts may be missing. */
        return 1;                   /* 1 = Wrong length. */
    if (textBuff[n-1] != '\n')
        return 2;                   /* 2 = Missing newline. */
    if ((textBuff[41] != ' ') ||
        (textBuff[69] != ' '))
        return 3;                   /* 3 = Missing blanks. */
    /*--------------------------------------- WLVAC: 11 */
    rem = textBuff[11];
    textBuff[11] = 0;
    PBIN->WLVAC = strtod(&textBuff[0], (char**)NULL);
    textBuff[11] = rem;
    /*--------------------------------------- LOGGF: 7  */
    rem = textBuff[18];
    textBuff[18] = 0;
    PBIN->LOGGF = strtod(&textBuff[11], (char**)NULL);
    textBuff[18] = rem;
    /*--------------------------------------- ELEM: 6  (to int) */
    rem = textBuff[24];
    textBuff[24] = 0;
    PBIN->ELEM = (short)(0.5 + (100.0 * strtod(&textBuff[18], (char**)NULL)));
    textBuff[24] = rem;
    /*--------------------------------------- E_LOW: 12 */
    rem = textBuff[36];
    textBuff[36] = 0;
    PBIN->LOW.E = strtold(&textBuff[24], (char**)NULL);
    textBuff[36] = rem;
    /*--------------------------------------- J_LOW: 5 */
    rem = textBuff[41];
    textBuff[41] = 0;
    PBIN->LOW.J = (float)strtod(&textBuff[36], (char**)NULL);
    textBuff[41] = rem;
    /* [41] = blank. */
    /*--------------------------------------- LABEL_LOW: 10 */
    n = 52;
    while ((textBuff[n-1] == ' ') && (n > 42))
        n--;
    rem = textBuff[n];
    textBuff[n] = 0;
    strcpy(PBIN->LOW.LABEL, &textBuff[42]);
    textBuff[n] = rem;
    /*--------------------------------------- E_UPP: 12 */
    rem = textBuff[64];
    textBuff[64] = 0;
    PBIN->UPP.E = strtold(&textBuff[52], (char**)NULL);
    textBuff[64] = rem;
    /*--------------------------------------- J_UPP: 5 */
    rem = textBuff[69];
    textBuff[69] = 0;
    PBIN->UPP.J = (float)strtod(&textBuff[64], (char**)NULL);
    textBuff[69] = rem;
    /* [69] = blank. */
    /*--------------------------------------- LABEL_UPP: 10 */
    n = 80;
    while ((textBuff[n-1] == ' ') && (n > 70))
        n--;
    rem = textBuff[n];
    textBuff[n] = 0;
    strcpy(PBIN->UPP.LABEL, &textBuff[70]);
    textBuff[n] = rem;
    /*--------------------------------------- REF: 4 */
    n = 102;
    while ((textBuff[n-1] == ' ') && (n > 98))
        n--;
    rem = textBuff[n];
    textBuff[n] = 0;
    strcpy(PBIN->REF, &textBuff[98]);
    textBuff[n] = rem;
    /*---------------------------------------*/
    return 0;
}

/*----------------------------------------------------------------------------*/
int pack_data(FILE* msg)
{
    #define kLineLen  (160)
    FILE              *textDatabaseFP,          /* Input file. */
                      *pbinDatabaseFP;          /* Output file. */
    char              lineBuff[kLineLen + 4];   /* +4 for LF, CR and Z. */
    pieterPacked      pbinBuff;
    unsigned long     linesParsed;
    unsigned long     checksum;  /* Used to test the original testfile. */
    int               e = 0;

    /* We first check whether the binary database is already there. If so,  
       we just assume it is all right then and we return +1 to the caller. */
    pbinDatabaseFP = fopen(kST_bin_base, "rb");
    if (pbinDatabaseFP != NULL)
        {
        linesParsed = 0L; /* Count the number of pbinBuff structs. */
        while (1 == fread(&pbinBuff, sizeof(pbinBuff), 1, pbinDatabaseFP))
            linesParsed++;
        fclose(pbinDatabaseFP);
        if (linesParsed == kKuruczLinesTotal)
            {
            if (msg != NULL)
                fprintf(msg, "Binary database '%s' seems present.\n", kST_bin_base);
            return 1;
            }
        if (msg != NULL)
            fprintf(msg, "Overwriting wrong binary database...\n");
        }
    textDatabaseFP = fopen(kST_txt_base, "rt");         /* OPEN FILE TEXTFILE. */
    if (textDatabaseFP != NULL)
        {
        if (msg != NULL)
            fprintf(msg, "Reading textfile '%s'.\n", kST_txt_base);
        pbinDatabaseFP = fopen(kST_bin_base, "wb");     /* CREATE NEW FILE. */
        if (pbinDatabaseFP != NULL)
            {
            if (msg != NULL)
                fprintf(msg, "Writing binary file '%s' (%d bytes per record).\n",
                             kST_bin_base, (int)sizeof(pbinBuff));
            linesParsed = 0UL;
            checksum = 0UL;
            while (fgets(lineBuff, kLineLen + 2, textDatabaseFP) != NULL) /* +1 +1 */
                {
                e = processLine(lineBuff, &pbinBuff, &checksum);   /* READ FROM FILE. */
                if (e)
                    {
                    if (msg != NULL)
                        fprintf(msg, "Error: processLine()=%d, at line #%ld!\n", e, 1L+linesParsed);
                    e = -4; goto abortParse;
                    }
                if ((pbinBuff.ELEM < 100) || (pbinBuff.ELEM > 9201))    /* Just warn: */
                    if (msg) fprintf(msg, "Error: element out of range, line #%ld!\n", 1L+linesParsed);
                if ((pbinBuff.WLVAC < 2.50) || (pbinBuff.WLVAC > 754956.3) ||
                    (pbinBuff.LOGGF > 3.75) || (pbinBuff.LOGGF < -18.0))
                    {
                    if (msg != NULL)
                        fprintf(msg, "#%6ld: ELEM=%6.2f; WLVAC=%11.4f; LOGGF=%7.3f; REF=%s;\n\
         E_LOW=%12.3Lf J_LOW=%5.2f %s;\n\
         E_UPP=%12.3Lf J_UPP=%5.2f %s.\n", 1L+linesParsed,
                        (double)pbinBuff.ELEM / 100.0, pbinBuff.WLVAC, pbinBuff.LOGGF, pbinBuff.REF,
                        pbinBuff.LOW.E, pbinBuff.LOW.J, pbinBuff.LOW.LABEL,
                        pbinBuff.UPP.E, pbinBuff.UPP.J, pbinBuff.UPP.LABEL);
                    }
                fwrite(&pbinBuff, sizeof(pbinBuff), 1, pbinDatabaseFP); /* WRITE TO FILE. */
                linesParsed++;                                          /* PLATFORM DEPENDANT! */
                }
            if (linesParsed != kKuruczLinesTotal)
                {
                if (msg != NULL)
                    fprintf(msg, "Wrong number of lines in database: linesParsed=%ld.\n", linesParsed);
                e = -3;
                }
            else if (checksum != kKuruczChecksum)
                {
                if (msg != NULL)
                    fprintf(msg, "Wrong checksum in original textfile: %ld.\n", checksum);
                e = -33;
                }
            else if (msg != NULL)
                fprintf(msg, "\
Parsed right amount of lines (%ld),\n\
binary file is expected to be %ld bytes in size,\n\
checksum tested (%lld).\n", linesParsed, linesParsed * sizeof(pbinBuff), (long long)checksum);
abortParse: fclose(pbinDatabaseFP);
            }
        else
            {
            if (msg != NULL)
                fprintf(msg, "Cannot create file '%s'.\n", kST_bin_base);
            e = -2;
            }
        fclose(textDatabaseFP);
        }
    else
        {
        if (msg != NULL)
            fprintf(msg, "Cannot open file '%s'.\n", kST_txt_base);
        e = -1;
        }
    return e;
}

/*----------------------------------------------------------------------------*/
void printLineEssential(pieterPacked* line, FILE* to)
{
    if (to != NULL)                            /* Called in composition-funcs.*/
        fprintf(to, "ELEM=%6.2f; WLVAC=%11.4f; LOGGF=%7.3f; REF=%4s; \
E_LOW=%12.3Lf J_LOW=%5.2f %10s; E_UPP=%12.3Lf J_UPP=%5.2f %10s.\n",
            (double)line->ELEM / 100.0, line->WLVAC, line->LOGGF, line->REF,
            line->LOW.E, line->LOW.J, line->LOW.LABEL,
            line->UPP.E, line->UPP.J, line->UPP.LABEL);
}

/*-------------------------------------------------------------------------------------*/
void moreLineData(pieterPacked* essential, derivedData* derived)
{
    double sig, sigquad, nair;      /* No energy levels in eV, we don't need those.    */
    /*---------------------------------------------------------------------------------*/
    /* Dispersion of air / "air wavelength": B. Edlen, Metrologia V. 2, p. 71 (1966):  */
    /*---------------------------------------------------------------------------------*/
    sig = 1.0e3 / essential->WLVAC;                         /* from nm to 1/micrometer */
    sigquad = sig * sig * 1.0e8;
    nair = 1.0000834213 + 2406030.0 / (1.3e10 - sigquad) + 15997.0 / (3.89e9 - sigquad);
    if(essential->WLVAC < 200.0)                            /* nm */
        {
        derived->wlVacCorrected = essential->WLVAC;         /* WLVAC was only for vac. */
        derived->wlAirCorrected = essential->WLVAC / nair;  /* when below 200 nm!      */
        }
    else
        {
        derived->wlVacCorrected = essential->WLVAC * nair;  /* Not absolutely correct, */
        derived->wlAirCorrected = essential->WLVAC;         /* but accuracy should be  */
        }                                                   /* sufficient.             */
    /*---------------------------------------------------------------------------------*/
    /* Adapted from source file "read_pac.c" from Peter L. Smith, Claas Heise,         */
    /* Jim R. Esmond, Robert L. Kurucz (Hannover University / Harvard).                */
    /* Reference for all constants & conversions: Spectrophysics, Anne P. Thorne,      */
    /* 2nd ed., Chapman and Hall, London - New York (1988). Conversion constants for   */
    /* corrections from Astrophysical Quantities, C.W.Allen, 3rd Ed., 1976.            */
    /* A-value: A_{21} = (6.670e13 * g1/g2 / lambda/lambda) * (pow(10.0,log(gf)) / g1) */
    /*                   g1: (2*j_low + 1), g2: (2*j_upp + 1), log_gf = log (g1*f).    */
    /*---------------------------------------------------------------------------------*/
    derived->waveNumbers = essential->UPP.E - essential->LOW.E;         /* cm^-1       */
    /* Note these energy-differences are sometimes NEGATIVE!!(??)                      */
    /* (Except from sign, corresponding wavelengths seem however to be ok.)            */
    if (derived->waveNumbers < 0.0)
        {
        derived->waveNumbers = -derived->waveNumbers;
        /* printf("(\"moreLineData()\" negated frequency)\n");   */
        }
    /*       photonEnergy not used directly (yet/any longer):
    derived->photonEnergy = 1.2398544e-4 * derived->waveNumbers; // cm^-1 -> eV
    derived->frequency = 2.417965e14 * derived->photonEnergy;    // s^-1
    */
    /*  Calculate frequency directly from waveNumbers to reduce rounding errors: */
    derived->frequency = 29979245442.96 * derived->waveNumbers;

    /* Thus (Pieter, 2008):
    
       A-value = (6.670e13 * g1/g2 / lambda/lambda) * (pow(10.0,log(gf)) / g1)
            g1 = 2*j_low + 1
            g2 = 2*j_upp + 1
        log_gf = log(g1*f)
        lambda = corrected wavelength in vacuum?!
        
        
        A-value = (6.670e13 * (2*j_low + 1) / (2*j_upp + 1) / lambda / lambda)
                * (pow(10.0,log(gf)) / (2*j_low + 1))
        
        Hmmmmmm...... many doubts whether the following was correct:
    
    */
    derived->Avalue = 6.670e13 / (2.0 * essential->UPP.J + 1.0) / 
                      derived->wlVacCorrected / derived->wlVacCorrected *
                                                pow(10.0, essential->LOGGF);
}

/*-------------------------------------------------------------------------------------*/
void printLineDerived(derivedData* line, FILE* to)
{
    if (to != NULL)
        fprintf(to, "           vacCorr=%11.4f; airCorr=%11.4f; frq=%.10Le; A=%.4e.\n",
                line->wlVacCorrected, line->wlAirCorrected, line->frequency, line->Avalue);
}