#include <assert.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>

#include <deque>
#include <fstream>
#include <sstream>
#include <string>

using namespace std;

struct Record {
    double d_time;
    double d_weight;
    double d_level;
    double d_weightedLevel;
    double d_threeYearWeightedLevel;
};

int main(int argc, char **argv)
{
    assert(2 == argc);

    const string inFile       = argv[1];
    const string cmdFn        = "tmp.sea_level.gnuplot";

    const string rawFn        = "tmp.sea_level_raw.plot";
    const string threeYrFn    = "tmp.sea_level_3yr.plot";
    const string tenYrDeltaFn = "tmp.sea_level_3yrDelta.plot";
    const string initialFn    = "tmp.sea_level_1900.plot";

    fstream in(inFile,             ios_base::in);
    fstream cmdOut(cmdFn,          ios_base::out);

    fstream rawOut(rawFn,          ios_base::out);
    fstream threeYrOut(threeYrFn,  ios_base::out);
    fstream riseOut(tenYrDeltaFn,  ios_base::out);
    fstream initialRise(initialFn, ios_base::out);

    assert(cmdOut.     good());
    assert(in.         good());
    assert(rawOut.     good());
    assert(threeYrOut. good());
    assert(riseOut.    good());
    assert(initialRise.good());

    const string chartName    = "Sea Level";

    // create the gnu plot command
    {
        cmdOut << "plot ";

        cmdOut << "'";
        cmdOut << rawFn;
        cmdOut << "' with lines lw 2 lc rgb '#00ff00' title '";
        cmdOut << chartName;
        cmdOut << ' ';
        cmdOut << "Raw Relative to Jan 4th 1993 (mm)',";

        cmdOut << "'";
        cmdOut << threeYrFn;
        cmdOut << "' with lines lw 2 lc rgb '#0000ff' title '";
        cmdOut << chartName;
        cmdOut << ' ';
        cmdOut << "(Three Year Smooth) (mm)',";

        cmdOut << "'";
        cmdOut << tenYrDeltaFn;
        cmdOut << "' with lines lw 2 lc rgb '#ff0000' title '";
        cmdOut << chartName;
        cmdOut << ' ';
        cmdOut << "Rise Rate (mm / 10 years)',";

        cmdOut << "'";
        cmdOut << initialFn;
        cmdOut << "' with lines lw 2 lc rgb '#a52a2a' title '";
        cmdOut << chartName;
        cmdOut << ' ';
        cmdOut << "Rise Rate 1900 - 1993 (mm / 10 years)'\n";
    }

    string line;
    deque<Record> records;

    const double nullDouble = -1e20, threshold = -0.9e20;
    const double rise1900_1993 = 115;    // Read manually from a chart

    double startTime = -1, endTime = -1;
    double startLevel = nullDouble, level1900 = nullDouble;
    while (getline(in, line)) {
        if ("HDR" == line.substr(0, 3)) {
            // skip lines that begin with "HDR"

            continue;
        }

        Record newRec = { 0, 0, nullDouble, 0, nullDouble };

        int altimeterType, iDum;
        double rDum;

        istringstream ss(line);
	ss >> altimeterType;
        ss >> iDum;
        ss >> newRec.d_time;
        ss >> iDum;
        ss >> newRec.d_weight;
        ss >> rDum;
        ss >> rDum;
        ss >> rDum;
        ss >> rDum;
        ss >> rDum;
        ss >> rDum;
        ss >> rDum;
	assert(!ss.eof());
        ss >> newRec.d_level;
	assert(ss.eof());
        assert(0 == altimeterType || 999 == altimeterType);
        assert(threshold < newRec.d_level);
        if (0 != altimeterType) {
            // something's weird about the lines that start with '999'

            continue;
        }

        if (records.empty()) {
            startTime  = newRec.d_time;
	    assert(abs(startTime - 1993) < 0.1);
            startLevel = (int) newRec.d_level;
            level1900 = startLevel - rise1900_1993;
        }
        endTime = newRec.d_time;

        newRec.d_weightedLevel = newRec.d_weight * newRec.d_level;
        records.push_back(newRec);
    }
    assert(1990 < startTime);
    assert(1990 < endTime);
    assert(startTime < endTime);
    assert(threshold < startLevel);
    assert(threshold < level1900);
    assert(1000 < records.size());

    const ptrdiff_t numRecords = records.size();
    for (int uu = 0; uu < numRecords; ++uu) {
        Record&       current = records[uu];
        const double  time    = current.d_time;

        rawOut << time << '\t' << (current.d_level - startLevel) << '\n';

        if (1.5 <= time - startTime && 1.5 <= endTime - time) {
            double totalWeight = 0, totalWeightedLevel = 0,
                                                         totalWeightedTime = 0;

            for (int vv = 0; vv < numRecords; ++vv) {
                if (abs(time - records[vv].d_time) <= 1.5) {
                    const Record& vRec = records[vv];
                    totalWeight        += vRec.d_weight;
                    totalWeightedLevel += vRec.d_weightedLevel;
		    totalWeightedTime  += vRec.d_weight * vRec.d_time;
                }
            }

            if (0 < totalWeight) {
                current.d_threeYearWeightedLevel =
                                              totalWeightedLevel / totalWeight;

                const double level = totalWeightedLevel / totalWeight;
                const double weightedTime = totalWeightedTime / totalWeight;

                threeYrOut << weightedTime << '\t' << (level - startLevel) <<
                                                                          '\n';
            }
        }
    }

    for (unsigned uu = 0; uu < numRecords; ++uu) {
        Record&       current = records[uu];
        const double  time    = current.d_time;

        if (threshold < current.d_threeYearWeightedLevel) {
            for (int vv = uu; vv < numRecords; ++vv) {
                Record& vRec = records[vv];

                if (vRec.d_time - time >= 3.0) {
                    if (threshold < vRec.d_threeYearWeightedLevel) {
                        const double deltaTime = vRec.d_time - time;
                        const double diff3 = vRec.d_threeYearWeightedLevel -
                                              current.d_threeYearWeightedLevel;

                        riseOut << ((vRec.d_time + time) / 2) << '\t' <<
                                              (10 * diff3 / deltaTime) << '\n';
                    }
                    break;
                }
            }
        }
    }

    const double initialRiseRate = 10 * (startLevel - level1900) /
                                                          (startTime - 1900.0);
    initialRise << startTime << '\t' << initialRiseRate << '\n';
    initialRise << endTime   << '\t' << initialRiseRate << '\n';

    assert(in.         eof());
    assert(cmdOut.     good());

    assert(rawOut.     good());
    assert(threeYrOut. good());
    assert(riseOut.    good());
    assert(initialRise.good());

    in.         close();
    cmdOut.     close();

    rawOut.     close();
    threeYrOut. close();
    riseOut.    close();
    initialRise.close();

    system(("gnuplot -p " + cmdFn).c_str());
}