/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#include "sac.h"

#ifdef GCTL_SEISMIC_MATHGL

gctl::SAC_PLOT::SAC_PLOT()
{
    sig_ptr = nullptr; hd_ptr = nullptr;
}

gctl::SAC_PLOT::~SAC_PLOT(){}

void gctl::SAC_PLOT::link_source(SIG_UNIT *sig, SAC_HD *hd, std::string filename)
{
    sig_ptr = sig; hd_ptr = hd; file = filename;
    return;
}

int gctl::SAC_PLOT::Draw(mglGraph *gr)
{
    gctl::array<float> x;
    x.sequence((float) 0.0, hd_ptr->delta, sig_ptr->val.size());

    mglData x_d, y_d;
    x_d.Set(x.get(), x.size());
    y_d.Set(sig_ptr->val.get(), sig_ptr->val.size());

    float mean = sig_ptr->val.mean();
    std::string s_nk = hd_ptr->knetwk; s_nk.erase(s_nk.find_last_not_of(" ") + 1);
    std::string s_st = hd_ptr->kstnm; s_st.erase(s_st.find_last_not_of(" ") + 1);
    std::string s_cp = hd_ptr->kcmpnm; s_cp.erase(s_cp.find_last_not_of(" ") + 1);
    std::string leng = s_nk + "." + s_st + "." + s_cp;

    gr->SetSize(800, 400);
    gr->SetRanges(x.front(), x.back(), hd_ptr->xminimum - 0.5*(mean - hd_ptr->xminimum), hd_ptr->xmaximum + 0.5*(hd_ptr->xmaximum - mean));
    gr->AddLegend(leng.c_str(), "h");
    gr->Title(sig_ptr->t0.time_str().c_str());
    gr->Box();
    gr->Axis("x");
    gr->Axis("y");
    gr->Grid("xy", "h:");
    gr->Legend(2);
    gr->Label('x', "Time (s)", 0.0);
    gr->Plot(x_d, y_d, "B");

    if (file != "null")
    {
        std::string tp = file.substr(file.find_last_of('.') + 1);
        if (tp == "eps") gr->WriteEPS(file.c_str());
        if (tp == "png") gr->WritePNG(file.c_str());
        if (tp == "jpg") gr->WriteJPEG(file.c_str());
    }
    return 0;
}

#endif // GCTL_SEISMIC_MATHGL

gctl::SAC::SAC(){}

gctl::SAC::~SAC(){}

gctl::SAC::SAC(std::string filename, bool headonly) {read(filename, headonly);}

void gctl::SAC::read(std::string filename, bool headonly)
{
    std::string suf_str = filename.substr(filename.find_last_of('.') + 1);
    if (suf_str != "SAC" && suf_str != "SAC_ASC")
    {
        throw std::runtime_error("Invalid file suffix.");
    }

    std::ifstream fs;
    if (suf_str == "SAC_ASC")
    {
        gctl::open_infile(fs, filename);

        fs >> head.delta     >> head.depmin    >> head.depmax   >> head.unused1  >> head.odelta;    //0 - 4
        fs >> head.b         >> head.e         >> head.o        >> head.a        >> head.internal;  //5 - 9
        fs >> head.t0        >> head.t1        >> head.t2       >> head.t3       >> head.t4;        //10 - 14
        fs >> head.t5        >> head.t6        >> head.t7       >> head.t8       >> head.t9;        //15 - 19
        fs >> head.f         >> head.resp0     >> head.resp1    >> head.resp2    >> head.resp3;     //20 - 24
        fs >> head.resp4     >> head.resp5     >> head.resp6    >> head.resp7    >> head.resp8;     //25 - 29
        fs >> head.resp9     >> head.stla      >> head.stlo     >> head.stel     >> head.stdp;      //30 - 34
        fs >> head.evla      >> head.evlo      >> head.evel     >> head.evdp     >> head.mag;       //35 - 39
        fs >> head.user0     >> head.user1     >> head.user2    >> head.user3    >> head.user4;     //40 - 44
        fs >> head.user5     >> head.user6     >> head.user7    >> head.user8    >> head.user9;     //45 - 49
        fs >> head.dist      >> head.az        >> head.baz      >> head.gcarc    >> head.sb; 
        fs >> head.sdelta    >> head.depmen    >> head.cmpaz    >> head.cmpinc   >> head.xminimum;   
        fs >> head.xmaximum  >> head.yminimum  >> head.ymaximum >> head.adjtm    >> head.unused2;   
        fs >> head.unused3   >> head.unused4   >> head.unused5  >> head.unused6  >> head.unused7;  
        fs >> head.nzyear    >> head.nzjday    >> head.nzhour   >> head.nzmin    >> head.nzsec;     
        fs >> head.nzmsec    >> head.nvhdr     >> head.norid    >> head.nevid    >> head.npts;      
        fs >> head.nsnpts    >> head.nwfid     >> head.nxsize   >> head.nysize   >> head.unused8;  
        fs >> head.iftype    >> head.idep      >> head.iztype   >> head.unused9  >> head.iinst;     
        fs >> head.istreg    >> head.ievreg    >> head.ievtyp   >> head.iqual    >> head.isynth;    
        fs >> head.imagtyp   >> head.imagsrc   >> head.ibody    >> head.unused10 >> head.unused11;  
        fs >> head.unused12  >> head.unused13  >> head.unused14 >> head.unused15 >> head.unused16;  
        fs >> head.leven     >> head.lpspol    >> head.lovrok   >> head.lcalda   >> head.unused17;  
        fs >> head.kstnm  >> head.kevnm;           
        fs >> head.khole  >> head.ko     >> head.ka;               
        fs >> head.kt0    >> head.kt1    >> head.kt2;              
        fs >> head.kt3    >> head.kt4    >> head.kt5;              
        fs >> head.kt6    >> head.kt7    >> head.kt8;              
        fs >> head.kt9    >> head.kf     >> head.kuser0;           
        fs >> head.kuser1 >> head.kuser2 >> head.kcmpnm;
        fs >> head.knetwk >> head.kdatrd >> head.kinst;

        if (!headonly)
        {
            signal.val.resize(head.npts);
            for (int i = 0; i < head.npts; i++)
            {
                fs >> signal.val[i];
            }
        }
        fs.close();
    }
    else
    {
        gctl::open_infile(fs, filename, "", std::ios::binary);

        fs.read((char*)&head, sizeof(SAC_HD));
        head.kstnm[7] = '\0';  head.kevnm[15] = '\0';
        head.khole[7] = '\0';  head.ko[7]     = '\0'; head.ka[7]  = '\0';               
        head.kt0[7]   = '\0';  head.kt1[7]    = '\0'; head.kt2[7] = '\0';              
        head.kt3[7]   = '\0';  head.kt4[7]    = '\0'; head.kt5[7] = '\0';              
        head.kt6[7]   = '\0';  head.kt7[7]    = '\0'; head.kt8[7] = '\0';              
        head.kt9[7]   = '\0';  head.kf[7]     = '\0'; head.kuser0[7] = '\0';           
        head.kuser1[7] = '\0'; head.kuser2[7] = '\0'; head.kcmpnm[7] = '\0';
        head.knetwk[7] = '\0'; head.kdatrd[7] = '\0'; head.kinst[7]  = '\0';

        if (!headonly)
        {
            signal.val.resize(head.npts);
            array<float> tmp_s(head.npts);

            fs.read((char*)tmp_s.get(), head.npts*sizeof(float));
            for (size_t i = 0; i < head.npts; i++)
            {
                signal.val[i] = (double) tmp_s[i];
            }
        }
        fs.close();
    }

    signal.delta = head.delta;
    signal.t0.set_time(head.nzyear, head.nzjday, head.nzhour, head.nzmin, head.nzsec, head.nzmsec);
    signal.te = signal.t0;
    signal.te.add_duration(head.e);
    sync();
    return;
}

void gctl::SAC::save(std::string filename) const
{
    std::string suf_str = filename.substr(filename.find_last_of('.') + 1);
    if (suf_str != "SAC" && suf_str != "SAC_ASC")
    {
        throw std::runtime_error("[GCTL_SEISMIC] Invalid file suffix.");
    }

    if (head.npts != signal.val.size())
    {
        throw std::runtime_error("[GCTL_SEISMIC] Record length and the head section don't match.");
    }

    std::ofstream fs;
    if (suf_str == "SAC_ASC")
    {
        gctl::open_outfile(fs, filename);

        fs << std::setprecision(7) << std::setw(15) << head.delta     << std::setw(15) << head.depmin    << std::setw(15) << head.depmax   << std::setw(15) << head.unused1  << std::setw(15) << head.odelta << "\n";    //0 - 4
        fs << std::setprecision(7) << std::setw(15) << head.b         << std::setw(15) << head.e         << std::setw(15) << head.o        << std::setw(15) << head.a        << std::setw(15) << head.internal << "\n";  //5 - 9
        fs << std::setprecision(7) << std::setw(15) << head.t0        << std::setw(15) << head.t1        << std::setw(15) << head.t2       << std::setw(15) << head.t3       << std::setw(15) << head.t4 << "\n";        //10 - 14
        fs << std::setprecision(7) << std::setw(15) << head.t5        << std::setw(15) << head.t6        << std::setw(15) << head.t7       << std::setw(15) << head.t8       << std::setw(15) << head.t9 << "\n";        //15 - 19
        fs << std::setprecision(7) << std::setw(15) << head.f         << std::setw(15) << head.resp0     << std::setw(15) << head.resp1    << std::setw(15) << head.resp2    << std::setw(15) << head.resp3 << "\n";     //20 - 24
        fs << std::setprecision(7) << std::setw(15) << head.resp4     << std::setw(15) << head.resp5     << std::setw(15) << head.resp6    << std::setw(15) << head.resp7    << std::setw(15) << head.resp8 << "\n";     //25 - 29
        fs << std::setprecision(7) << std::setw(15) << head.resp9     << std::setw(15) << head.stla      << std::setw(15) << head.stlo     << std::setw(15) << head.stel     << std::setw(15) << head.stdp << "\n";      //30 - 34
        fs << std::setprecision(7) << std::setw(15) << head.evla      << std::setw(15) << head.evlo      << std::setw(15) << head.evel     << std::setw(15) << head.evdp     << std::setw(15) << head.mag << "\n";       //35 - 39
        fs << std::setprecision(7) << std::setw(15) << head.user0     << std::setw(15) << head.user1     << std::setw(15) << head.user2    << std::setw(15) << head.user3    << std::setw(15) << head.user4 << "\n";     //40 - 44
        fs << std::setprecision(7) << std::setw(15) << head.user5     << std::setw(15) << head.user6     << std::setw(15) << head.user7    << std::setw(15) << head.user8    << std::setw(15) << head.user9 << "\n";     //45 - 49
        fs << std::setprecision(7) << std::setw(15) << head.dist      << std::setw(15) << head.az        << std::setw(15) << head.baz      << std::setw(15) << head.gcarc    << std::setw(15) << head.sb << "\n"; 
        fs << std::setprecision(7) << std::setw(15) << head.sdelta    << std::setw(15) << head.depmen    << std::setw(15) << head.cmpaz    << std::setw(15) << head.cmpinc   << std::setw(15) << head.xminimum << "\n";   
        fs << std::setprecision(7) << std::setw(15) << head.xmaximum  << std::setw(15) << head.yminimum  << std::setw(15) << head.ymaximum << std::setw(15) << head.adjtm    << std::setw(15) << head.unused2 << "\n";   
        fs << std::setprecision(7) << std::setw(15) << head.unused3   << std::setw(15) << head.unused4   << std::setw(15) << head.unused5  << std::setw(15) << head.unused6  << std::setw(15) << head.unused7 << "\n";  
        fs << std::setprecision(7) << std::setw(10) << head.nzyear    << std::setw(10) << head.nzjday    << std::setw(10) << head.nzhour   << std::setw(10) << head.nzmin    << std::setw(10) << head.nzsec << "\n";     
        fs << std::setprecision(7) << std::setw(10) << head.nzmsec    << std::setw(10) << head.nvhdr     << std::setw(10) << head.norid    << std::setw(10) << head.nevid    << std::setw(10) << head.npts << "\n";      
        fs << std::setprecision(7) << std::setw(10) << head.nsnpts    << std::setw(10) << head.nwfid     << std::setw(10) << head.nxsize   << std::setw(10) << head.nysize   << std::setw(10) << head.unused8 << "\n";  
        fs << std::setprecision(7) << std::setw(10) << head.iftype    << std::setw(10) << head.idep      << std::setw(10) << head.iztype   << std::setw(10) << head.unused9  << std::setw(10) << head.iinst << "\n";     
        fs << std::setprecision(7) << std::setw(10) << head.istreg    << std::setw(10) << head.ievreg    << std::setw(10) << head.ievtyp   << std::setw(10) << head.iqual    << std::setw(10) << head.isynth << "\n";    
        fs << std::setprecision(7) << std::setw(10) << head.imagtyp   << std::setw(10) << head.imagsrc   << std::setw(10) << head.ibody    << std::setw(10) << head.unused10 << std::setw(10) << head.unused11 << "\n";  
        fs << std::setprecision(7) << std::setw(10) << head.unused12  << std::setw(10) << head.unused13  << std::setw(10) << head.unused14 << std::setw(10) << head.unused15 << std::setw(10) << head.unused16 << "\n";  
        fs << std::setprecision(7) << std::setw(10) << head.leven     << std::setw(10) << head.lpspol    << std::setw(10) << head.lovrok   << std::setw(10) << head.lcalda   << std::setw(10) << head.unused17 << "\n";  
        fs << head.kstnm  << " " << head.kevnm << " \n";           
        fs << head.khole  << " " << head.ko     << " " << head.ka << " \n";               
        fs << head.kt0    << " " << head.kt1    << " " << head.kt2 << " \n";              
        fs << head.kt3    << " " << head.kt4    << " " << head.kt5 << " \n";              
        fs << head.kt6    << " " << head.kt7    << " " << head.kt8 << " \n";              
        fs << head.kt9    << " " << head.kf     << " " << head.kuser0 << " \n";           
        fs << head.kuser1 << " " << head.kuser2 << " " << head.kcmpnm << " \n";
        fs << head.knetwk << " " << head.kdatrd << " " << head.kinst << " \n";

        int i = 0;
        while (i < head.npts)
        {
            fs << std::setprecision(7) << std::setw(15) << signal.val[i];
            if ((i+1)%5 == 0) fs << "\n";
            i++;
        }
        fs.close();
    }
    else
    {
        gctl::open_outfile(fs, filename, "", std::ios::binary);

        fs.write((char*)&head, sizeof(SAC_HD));

        array<float> tmp_s(head.npts);
        for (size_t i = 0; i < head.npts; i++)
        {
            tmp_s[i] = (float) signal.val[i];
        }

        fs.write((char*)tmp_s.get(), head.npts*sizeof(float));
        fs.close();
    }
    return;
}

void gctl::SAC::info() const
{
    std::clog << "SAC Info\n-------------\n";
    std::clog << "Station: " << head.kstnm << "\nNetwork: " << head.knetwk << "\nCompnent: " << head.kcmpnm << "\n";
    std::clog << "Start-Time: " << signal.t0.time_str() << "\nEnd-Time: " << signal.te.time_str() << "\n";
    std::clog << "Delta: " << head.delta << " s\nDuration: " << (double) (head.npts-1)*head.delta << " s\nSingal-Length: " << head.npts << "\n";
    std::clog << "-------------\n";
    return;
}

void gctl::SAC::sync()
{
    head.nzyear = signal.t0.year;
    head.nzjday = signal.t0.julday;
    head.nzhour = signal.t0.hour;
    head.nzmin = signal.t0.mint;
    head.nzsec = signal.t0.sec;
    head.nzmsec = signal.t0.msec;
    head.npts = signal.val.size();
    head.delta = signal.delta;
    head.e = (head.npts - 1)*head.delta;

    double xmin = signal.val[0];
    double xmax = signal.val[0];
    for (size_t i = 1; i < signal.val.size(); i++)
    {
        xmin = std::min(xmin, signal.val[i]);
        xmax = std::max(xmax, signal.val[i]);
    }
    
    head.xmaximum = xmax;
    head.xminimum = xmin;
    return;
}

void gctl::SAC::cut_section(SAC &ret, double st, double et) const
{
    signal.cut_section(ret.signal, st, et);
    ret.head = head;
    ret.sync();
    return;
}

void gctl::SAC::cut_section(SAC &ret, const UTC_TIME &st, double len) const
{
    double stt = st.diff_sec(signal.t0);
    double ett = stt + len;
    return cut_section(ret, stt, ett);
}

void gctl::SAC::stack_sac(const gctl::array<SAC*> &list)
{
    int ls = list.size();
    if (list.empty())
    {
        throw std::runtime_error("[GCTL_SEISMIC] Invalid stacking list for gctl::SAC::stack_sac(...)");
    }

    array<SIG_UNIT*> sigs(list.size());
    for (size_t i = 0; i < ls; i++)
    {
        sigs[i] = &list[i]->signal;
    }

    signal.stack_sig_units(sigs);

    head = list[0]->head;
    strcpy(head.kstnm, "STACK  ");
    sync();
    return;
}

void gctl::SAC::down_sampling(int factor)
{
    signal.down_sampling(factor);
    sync();
    return;
}

void gctl::SAC::remove_mean()
{
    signal.remove_mean();
    sync();
    return;
}

void gctl::SAC::remove_outliers(double factor)
{
    signal.remove_outliers(factor);
    sync();
    return;
}

void gctl::SAC::normalize(double norm)
{
    signal.normalize(norm);
    sync();
    return;
}

#ifdef GCTL_SEISMIC_MATHGL

int gctl::SAC::plot(std::string filename)
{
    plt.link_source(&signal, &head, filename);
    mglFLTK gr(&plt, "SAC plot");
    return gr.Run();
}

#endif // GCTL_SEISMIC_MATHGL