Command.cpp

/**
 * @file Command.cpp
 * @brief Main interaction class, contains all comands logics that are activated py the interpreter.
 * @copyright Copyright (C) 2025 ForeFire, Fire Team, SPE, CNRS/Universita di Corsica.
 * @license This program is free software; See LICENSE file for details. (See LICENSE file).
 * @author Jean‑Baptiste Filippi — 2025
 */

#include "Command.h"
#include "colormap.h"
#include <sstream>
#include <dirent.h>
#include <cmath> 
#include <fstream>
#define STB_IMAGE_WRITE_IMPLEMENTATION
#include "stb_image_write.h"

using namespace std;

namespace libforefire
{

    size_t Command::currentLevel = 0;

    bool Command::init = true;
    bool Command::currentFrontCompleted = false;

    double Command::startTime = 0;
    double Command::endTime = 0;

    bool Command::firstCommand = true;
    size_t Command::refTabs = 0;

    FFPoint *Command::lastReadLoc = 0;
    FireNode *Command::previousNode = 0;
    FireNode *Command::leftLinkNode = 0;
    FireNode *Command::rightLinkNode = 0;

    double Command::bmapOutputUpdate = 0;
    int Command::numBmapOutputs = 0;
    double Command::refTime = 0;
    int Command::numAtmoIterations = 0;

    const string Command::stringError = "1234567890";
    const FFPoint Command::pointError = FFPoint(1234567890., 1234567890., 0);
    const FFVector Command::vectorError = FFVector(1234567890., 1234567890.);

    vector<string> Command::outputDirs;

    Command::Session Command::currentSession =
        {
            SimulationParameters::GetInstance(),
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            &cout,
            0,
            0,
    };

    const Command::commandMap Command::translator = Command::makeCmds();

    // Defaults constructor and destructor for the 'Command' abstract class
    Command::Command()
    {
    }

    Command::~Command()
    {
    }

    void Command::increaseLevel()
    {
        currentLevel++;
    }

    void Command::decreaseLevel()
    {
        currentLevel--;
    }

    int Command::createDomain(const string &arg, size_t &numTabs)
    {

        size_t n = argCount(arg);
        if (n == 2)
        {
           
       // if there are 2 ption and it is "pgdNcFile" and a timestamp then it extracts NS, SW and t then calls again createDomain(..) with the 3 options
            string pgdNcFile = getString("pgdNcFile", arg);
            string timeStampDomain = getString("ISOdate", arg);
            cout << "pgdNcFile: " << pgdNcFile << endl;
            cout << "timeStampDomain: " << timeStampDomain << endl;
            if (pgdNcFile != stringError)
            {
                
              
                    int year, yday;
                    double secs;
                    std::vector<double> xhatValues, yhatValues;
                    double deltaY = 0.0, deltaX = 0.0;
                    SimulationParameters *simParam = SimulationParameters::GetInstance();
                    if (simParam->ISODateDecomposition(timeStampDomain, secs, year, yday))
                    {
                        simParam->setInt("refYear", year);
                        simParam->setInt("refDay", yday);
                        //simParam->setInt("refTime", secs);
                        simParam->setParameter("ISOdate", timeStampDomain);
                        cout<<"refYear: " << year << " refDay: " << yday << " refTime: " << secs << endl;
                    }
                    else
                    {
                        cout << "Error: Invalid date format "<< timeStampDomain<<" Expected YYYY-MM-DDTHH:MM:SSZ" << endl;
                        throw BadOption();
                    }

                    try
                    {
                    NcFile dataFile(pgdNcFile.c_str(), NcFile::read);
                    if (!dataFile.isNull())
                    {
                        // Retrieve and process the "domain" variable attributes.
                        NcVar XHAT = dataFile.getVar("XHAT");
                        NcVar YHAT = dataFile.getVar("YHAT");
                        
                        
                        // Assuming XHAT and YHAT are NcVar objects and that you have a way
                        // to read their data into a C++ container.
                       
                        size_t nXElements = XHAT.getDim(0).getSize();
                        xhatValues.resize(nXElements);
                        size_t nYElements = YHAT.getDim(0).getSize();
                        yhatValues.resize(nYElements);
                        
                        YHAT.getVar(&yhatValues[0]);
                        XHAT.getVar(&xhatValues[0]);
                    }
                }
                catch (NcException &e)
                {
                    std::cerr << e.what() << std::endl;
                }

                if (xhatValues.size() >= 2 && yhatValues.size() >= 2) {
                    deltaX = xhatValues[1] - xhatValues[0];
                    deltaY = yhatValues[1] - yhatValues[0];
                } else {
                    std::cerr << "Error: Not enough data in XHAT or YHAT variable." << std::endl;
                    throw BadOption();
                } 
                 
                std::ostringstream domStream;
                if (!xhatValues.empty() && !yhatValues.empty()) {
                    // Compute the last coordinate values adding the delta.
                    double neX = xhatValues.back() + deltaX;
                    double neY = yhatValues.back() + deltaY;
                    
                    domStream << "FireDomain[sw=(" << xhatValues.front() << "," << yhatValues.front() 
                            << ",0);ne=(" << neX << "," << neY << ",0);t=" << secs << "]";
                    std::string domCommand = domStream.str();
                    
                    // Optionally, set or use 'dom' as needed.
                    cout << "Domain string created: " << domCommand << endl;
                    ExecuteCommand(domCommand);
                    return normal;
                } else {
                    cout << "Error: Unable to create domain string due to insufficient data." << endl;
                }
            

            }
        }
        if (n >= 3)
        {
            FFPoint SW = getPoint("sw", arg);
            FFPoint NE = getPoint("ne", arg);
            double t = getFloat("t", arg);

            setStartTime(t);
            setReferenceTime(t);

            /* creating the domain */
            if (currentSession.fd != 0)
            {

                if (currentSession.fd->getDomainID() == 1)
                { // ID 1 is for the Fortran MPI rank 1
                    
                    currentSession.params->setParameter("runmode", "masterMNH");

                    currentSession.fdp = new FireDomain(t, SW, NE);
                    currentSession.fdp->setTimeTable(currentSession.tt);
                    currentSession.ff = currentSession.fdp->getDomainFront();

                    ostringstream ffOutputsFNAME;
                    ffOutputsFNAME << currentSession.params->getParameter("caseDirectory") << '/'
                                   << currentSession.params->getParameter("fireOutputDirectory") << '/'
                                   << currentSession.params->getParameter("outputFiles")
                                   << "." << currentSession.fdp->getDomainID();
                    currentSession.params->setParameter("PffOutputsPattern", ffOutputsFNAME.str());
                    currentSession.outStrRepp = new StringRepresentation(currentSession.fdp);
                    currentSession.outStrRepp->setOutPattern(ffOutputsFNAME.str());
                    if (currentSession.params->getInt("outputsUpdate") != 0)
                    {
                        currentSession.tt->insert(new FFEvent(currentSession.outStrRepp));
                    }
                }
            }
            else
            {
                currentSession.fd = new FireDomain(t, SW, NE);

                /* setting up the pointer to the domain */

                /* creating the related timetable and simulator */
                currentSession.tt = new TimeTable(new FFEvent(getDomain()));
                currentSession.sim = new Simulator(currentSession.tt, currentSession.fd->outputs);
                getDomain()->setTimeTable(currentSession.tt);
                /* linking to the domain front */
                currentSession.ff = getDomain()->getDomainFront();
                /* managing the outputs */
                ostringstream ffOutputsPattern;
                ffOutputsPattern << currentSession.params->getParameter("caseDirectory") << '/'
                                 << currentSession.params->getParameter("fireOutputDirectory") << '/'
                                 << currentSession.params->getParameter("outputFiles")
                                 << "." << getDomain()->getDomainID();
                currentSession.params->setParameter("ffOutputsPattern", ffOutputsPattern.str());
                currentSession.outStrRep = new StringRepresentation(currentSession.fd);
                if (currentSession.params->getInt("outputsUpdate") != 0)
                {
                    currentSession.tt->insert(new FFEvent(currentSession.outStrRep));
                }
                /* increasing the level */
                increaseLevel();
            }
            /* local copy of the reference time */
            return normal;
        }
        else
        {
            throw MissingOption(3 - n);
        }
    }

    int Command::startFire(const string &arg, size_t &numTabs)
    {
        double t = getDomain()->getTime();
        SimulationParameters *simParam = SimulationParameters::GetInstance();
    
        // Process time and date in any case.
        double nt = getFloat("t", arg);
        if (nt != FLOATERROR)
            t = nt;
    
        string date = getString("date", arg);
        if (date != stringError)
        {
            int year, yday;
            double secs;
            simParam->ISODateDecomposition(date, secs, year, yday);
            t = simParam->SecsBetween(simParam->getDouble("refTime"),
                                       simParam->getInt("refYear"),
                                       simParam->getInt("refDay"),
                                       secs, year, yday);
            if (t<0){
                cout << "WARNING: Trying to set an ignition at date "<< date<<" before reference date at " << simParam->FormatISODate(simParam->getDouble("refTime"), simParam->getInt("refYear"), simParam->getInt("refDay")) << endl;
                t=0;     
            }
        }

        string type = "none";
        if (getString("polyencoded",arg) != stringError) type="polyencoded";
        if (getString("geojson",arg) != stringError)type="geojson";
        if (getString("kml",arg) != stringError)type="kml";
        if (getString("points",arg) != stringError)type="points";
        if (getString("loc",arg) != stringError)type="loc";
        if (getString("lonlat",arg) != stringError) type="lonlat";

        
        FireDomain *refDomain = getDomain();
        if (type == "geojson")
        {
            /*-----------------------------------------------------------------
             *  Decide whether the 'geojson' argument is an inline GeoJSON text
             *  or a filename.  If it's a filename, read the file so that the
             *  rest of the code always works on the string `geojsonText`.
             *-----------------------------------------------------------------*/
            std::string geojsonParam = getString("geojson", arg);
            std::string geojsonText;

            if (geojsonParam != stringError && !geojsonParam.empty())
            {
                /* Trim simple whitespace. */
                std::string trimmed = geojsonParam;
                trimmed.erase(0, trimmed.find_first_not_of(" \t\r\n"));
                trimmed.erase(trimmed.find_last_not_of(" \t\r\n") + 1);

                /* Triple‑quoted inline string?  e.g.  \"\"\"{ ... }\"\"\"   */
                bool tripleQuoted = trimmed.size() >= 6 &&
                                    trimmed.substr(0, 3) == "\"\"\"" &&
                                    trimmed.substr(trimmed.size() - 3) == "\"\"\"";

                if (tripleQuoted)
                {
                    geojsonText = trimmed.substr(3, trimmed.size() - 6); // strip outer quotes
                }
                else
                {
                    /* Try to open it as a file path. */
                    std::ifstream jf(trimmed.c_str());
                    if (jf)
                    {
                        geojsonText.assign((std::istreambuf_iterator<char>(jf)),
                                           std::istreambuf_iterator<char>());
                    }
                    else
                    {
                        /* Not a file (or failed to open) – treat as raw inline GeoJSON. */
                        geojsonText = trimmed;
                    }
                }
            }
            else
            {
                /* No explicit parameter value; fall back to the whole argument. */
                geojsonText = arg;
            }

            size_t vpos = geojsonText.find("\"valid_at\"");
            if (vpos != std::string::npos)
            {
                size_t q1 = geojsonText.find('"', vpos + 10);          // first quote after :
                size_t q2 = (q1 != std::string::npos) ? geojsonText.find('"', q1 + 1) : std::string::npos;
                if (q1 != std::string::npos && q2 != std::string::npos)
                {
                    std::string iso = geojsonText.substr(q1 + 1, q2 - q1 - 1);
                    int year, yday;
                    double secs;
                    if (simParam->ISODateDecomposition(iso, secs, year, yday))
                    {
                        t = simParam->SecsBetween(simParam->getDouble("refTime"),
                                                  simParam->getInt("refYear"),
                                                  simParam->getInt("refDay"),
                                                  secs, year, yday);
                          
                        if (t<0){
                                cout << "WARNING: Adding contour at t=0 because was trying to set an ignition at date "<< iso<<" before reference data date at " << simParam->FormatISODate(simParam->getDouble("refTime"), simParam->getInt("refYear"), simParam->getInt("refDay")) << endl;
                                t=0;     
                        }

                    }
                }
            }

            /* --------------------------------------------------------------
             *  Only parse inside the "coordinates" array, ignore everything
             *  else (timestamps, properties, …) to avoid spurious numbers.
             * -------------------------------------------------------------- */
            size_t coordKey = geojsonText.find("\"coordinates\"");
            if (coordKey == std::string::npos)
            {
                std::cout << "Error: \"coordinates\" key not found in GeoJSON." << std::endl;
                return normal;
            }
            size_t startBracket = geojsonText.find('[', coordKey);
            if (startBracket == std::string::npos)
            {
                std::cout << "Error: '[' after \"coordinates\" not found." << std::endl;
                return normal;
            }

            /* 2) Parse all coordinate triples belonging to rings.
             *    Depth == 2 (MultiPolygon) or 1 (Polygon) signals ring
             *    closure.                                               */
            std::vector<std::vector<FFPoint>> polygons;
            std::vector<FFPoint>               currentRing;
            std::string                        numbuf;
            std::vector<double>                triple;
            int                                depth = 0;

            double refLon     = getDomain()->getRefLongitude();
            double refLat     = getDomain()->getRefLatitude();
            double mPerDegLon = getDomain()->getMetersPerDegreesLon();
            double mPerDegLat = getDomain()->getMetersPerDegreeLat();

            auto flushNumber = [&](void)
            {
                if (!numbuf.empty())
                {
                    try
                    {
                        triple.push_back(std::stod(numbuf));
                    }
                    catch (...) { /* silently ignore bad numbers */ 
                    cout << "Error: Invalid number in GeoJSON coordinates: " << numbuf << std::endl;
                    }
                    numbuf.clear();

                    if (triple.size() == 3)
                    {
                        const double lon = triple[0];
                        const double lat = triple[1];
                        const double alt = triple[2];
                        const double x   = (lon - refLon) * mPerDegLon;
                        const double y   = (lat - refLat) * mPerDegLat;

                        /* Add the point to the current ring. */
                        currentRing.emplace_back(x, y, alt);
                        triple.clear();
                    }
                }
            };

            for (size_t idx = startBracket; idx < geojsonText.size(); ++idx)
            {
                char c = geojsonText[idx];

                if ((c >= '0' && c <= '9') || c == '-' || c == '+' || c == '.' ||
                    c == 'e' || c == 'E')
                {
                    numbuf.push_back(c);
                }
                else
                {
                    flushNumber();

                    if (c == '[')
                        ++depth;
                    else if (c == ']')
                    {
                        --depth;
                        /* Ring terminates when we just closed a bracket
                         * that returns us to depth 2 (MultiPolygon) or
                         * 1 (Polygon).                                   */
                        if ((depth == 2 || depth == 1) && !currentRing.empty())
                        {
                            /* Remove duplicate closing point, if present. */
                            if (currentRing.size() > 1 &&
                                currentRing.front().distance2D(currentRing.back()) < 1e-6)
                                currentRing.pop_back();

                            polygons.push_back(currentRing);
                            currentRing.clear();
                        }
                        /* Finished the outermost coordinates array?      */
                        if (depth == 0 && idx > startBracket)
                            break;
                    }
                }
            }
            flushNumber();
            if (!currentRing.empty())
                polygons.push_back(currentRing);

            if (polygons.empty())
            {
                std::cout << "Error: No coordinates found in GeoJSON string." << std::endl;
                return normal;
            }

     
            /* 3) Insert outer ring then inner rings as fire fronts.      */
            double fdepth = currentSession.params->getDouble("initialFrontDepth");
            double kappa  = 0.0;
            FFVector defaultVel(0, 0);

            FireFront  *contfront = currentSession.ff->getContFront();

            for (size_t p = 0; p < polygons.size(); ++p)
            {
                /* First polygon continues from contfront, subsequent ones
                 * are nested inside the current front.                   */
                if (p == 0)
                    currentSession.ff = refDomain->addFireFront(t, contfront);
                else
                    currentSession.ff = refDomain->addFireFront(t, currentSession.ff);

                FireNode *lastnode = nullptr;
                for (auto it = polygons[p].rbegin(); it != polygons[p].rend(); ++it)
                {
                    lastnode = refDomain->addFireNode(*it, defaultVel, t,
                                                      fdepth, kappa,
                                                      currentSession.ff, lastnode);
                }
                completeFront(currentSession.ff);
            }

            return normal;   /* GeoJSON fully handled – no fall-through. */
        }
        if (type == "polyencoded" ||
            type == "kml" || type == "points")
        {
            // Get the polygon points.
            std::vector<FFPoint> polygon = getPoly(type, arg);
            if (polygon.empty())
                return normal;  // Handle error or malformed input as needed.
    
            double fdepth = currentSession.params->getDouble("initialFrontDepth");
            double kappa = 0.0;
            FireNode *lastnode = nullptr;
    
            // For each point in the polygon, add a fire node.
            // Velocity here is set to a default; adjust as needed.
            FireFront *contfront = currentSession.ff->getContFront();
            currentSession.ff = refDomain->addFireFront(t, contfront);
            for ( FFPoint &p : polygon)
            {
                FFVector defaultVel(0, 0);
                lastnode = refDomain->addFireNode(p, defaultVel, t, fdepth, kappa, currentSession.ff, lastnode);
            }
            completeFront(currentSession.ff);
        }
        else if (type == "lonlat" || type == "loc")
        {
            // Process a single point.
            FFPoint pos = getPoint("lonlat", arg);
            if (pos == pointError)
                pos = getPoint("loc", arg);
            
    
            if (refDomain->striclyWithinDomain(pos) & !refDomain->isBurnt(pos,t))
            {
                // Optional: If a current fire front exists, finalize it.
                if (currentSession.ff != 0)
                {
                    completeFront(currentSession.ff);
                }
    
                double perimRes = refDomain->getPerimeterResolution() * 2;
                int fdom = getInt("domain", arg);
                if (fdom == INTERROR)
                    fdom = 0;
                double fdepth = currentSession.params->getDouble("initialFrontDepth");
                double kappa = 0.0;
                FireFront *contfront = currentSession.ff->getContFront();
                currentSession.ff = refDomain->addFireFront(t, contfront);
                // Build a triangle around the point.
                FFVector vel1(0, 1), vel2(1, -1), vel3(-1, -1);
                FFVector diffP1 = perimRes * vel1;
                FFVector diffP2 = perimRes * vel2;
                FFVector diffP3 = perimRes * vel3;

                FFPoint pos1 = pos + diffP1.toPoint();
                FFPoint pos2 = pos + diffP2.toPoint();
                FFPoint pos3 = pos + diffP3.toPoint();
    
                // Optionally adjust the velocities.
                vel1 *= 0.1;
                vel2 *= 0.1;
                vel3 *= 0.1;
    
                FireNode *lastnode = refDomain->addFireNode(pos1, vel1, t, fdepth, kappa, currentSession.ff, 0);
                lastnode = refDomain->addFireNode(pos2, vel2, t, fdepth, kappa, currentSession.ff, lastnode);
                refDomain->addFireNode(pos3, vel3, t, fdepth, kappa, currentSession.ff, lastnode);
                completeFront(currentSession.ff);
       
            }
        }
        else
        {
            // Fallback behavior in case type is not recognized.
            // This might log an error or simply do nothing.
        }
    
        return normal;
    }


    std::vector<FFPoint> Command::getPoly(const std::string &opt, const std::string &arg)
    {
        std::vector<FFPoint> points;

       if (opt == "kml")
        {
            // Example input:
            // "8.810264,41.968841,0.14 8.810117,41.968757,0.03 8.810017,41.968692,0.01"
            std::istringstream iss(arg);
            double refLon = getDomain()->getRefLongitude();
            double refLat = getDomain()->getRefLatitude();
            double mPerDegLon = getDomain()->getMetersPerDegreesLon();
            double mPerDegLat = getDomain()->getMetersPerDegreeLat();
            std::string token;
            while (iss >> token)
            {
                // Replace commas with spaces.
                std::replace(token.begin(), token.end(), ',', ' ');
                std::istringstream ptStream(token);
                double lon, lat, alt;
                if (!(ptStream >> lon >> lat >> alt))
                    continue;
                double x = (lon - refLon) * mPerDegLon;
                double y = (lat - refLat) * mPerDegLat;
                points.push_back(FFPoint(x, y, alt));
            }
        }
        else if (opt == "polyencoded")
        {
            // Polyline encoded string. Standard algorithm decodes lat and lon.
            // Altitude is not provided (set to 0).
            double refLon = getDomain()->getRefLongitude();
            double refLat = getDomain()->getRefLatitude();
            double mPerDegLon = getDomain()->getMetersPerDegreesLon();
            double mPerDegLat = getDomain()->getMetersPerDegreeLat();

            int index = 0, len = static_cast<int>(arg.size());
            int lat = 0, lng = 0;
            while (index < len)
            {
                int b, shift = 0, resultInt = 0;
                do
                {
                    b = arg[index++] - 63;
                    resultInt |= (b & 0x1f) << shift;
                    shift += 5;
                } while (b >= 0x20 && index < len);
                int dlat = (resultInt & 1) ? ~(resultInt >> 1) : (resultInt >> 1);
                lat += dlat;

                shift = 0;
                resultInt = 0;
                do
                {
                    b = arg[index++] - 63;
                    resultInt |= (b & 0x1f) << shift;
                    shift += 5;
                } while (b >= 0x20 && index < len);
                int dlng = (resultInt & 1) ? ~(resultInt >> 1) : (resultInt >> 1);
                lng += dlng;

                double latd = lat * 1e-5;
                double lngd = lng * 1e-5;
                double x = (lngd - refLon) * mPerDegLon;
                double y = (latd - refLat) * mPerDegLat;
                points.push_back(FFPoint(x, y, 0));
            }
        }
        else if (opt == "points")
        {
            // Example input: (x,y,z),(x,y,z),...
            // Coordinates are assumed to be already in Cartesian space.
            size_t pos = 0;
            while (true)
            {
                pos = arg.find('(', pos);
                if (pos == std::string::npos)
                    break;
                size_t end = arg.find(')', pos);
                if (end == std::string::npos)
                    break;
                std::string pointStr = arg.substr(pos + 1, end - pos - 1);
                std::replace(pointStr.begin(), pointStr.end(), ',', ' ');
                std::istringstream issPoint(pointStr);
                double x, y, z;
                if (!(issPoint >> x >> y >> z))
                {
                    pos = end + 1;
                    continue;
                }
                points.push_back(FFPoint(x, y, z));
                pos = end + 1;
            }
        }

        return points;
    }

    int Command::createFireFront(const string &arg, size_t &numTabs)
    {

        size_t n = argCount(arg);
        if (n >= 1)
        {
            double t = getFloat("t", arg);
            if (numTabs == currentLevel)
            {

                /* first completing the previous front if needed */
                if (currentSession.ff != 0)
                {
                    completeFront(currentSession.ff);
                    currentFrontCompleted = false;
                }
                /* creating the new front */
                FireFront *contfront = currentSession.ff->getContFront();
                currentSession.ff = getDomain()->addFireFront(t, contfront);
            }
            else if (numTabs > currentLevel)
            {

                /* first completing the previous front */
                completeFront(currentSession.ff);
                currentFrontCompleted = false;
                /* creation of an inner front to the current one */
                FireFront *contfront = currentSession.ff;
                currentSession.ff = getDomain()->addFireFront(t, contfront);
                currentLevel = numTabs;
            }
            else if (numTabs < currentLevel)
            {

                /* first completing the previous front */
                completeFront(currentSession.ff);
                currentFrontCompleted = false;
                /* creation of a new front at a higher level then the current one */
                FireFront *contfront = currentSession.ff;
                for (size_t k = 0; k < currentLevel - numTabs; k++)
                {
                    contfront = contfront->getContFront();
                }
                currentSession.ff = getDomain()->addFireFront(t, contfront);
                currentLevel = numTabs;
            }
        }
        else
        {
            throw MissingTime();
        }

        return normal;
    }

    int Command::addFireNode(const string &arg, size_t &numTabs)
    {

        if (lastReadLoc == 0)
            lastReadLoc = new FFPoint(-numeric_limits<double>::infinity(), -numeric_limits<double>::infinity(), 0);

        size_t n = argCount(arg);
        double perimRes = getDomain()->getPerimeterResolution();
        if (n >= 3)
        {
            FFPoint pos = getPoint("loc", arg);
            FFVector vel = getVector("vel", arg);
            double t = getFloat("t", arg);
            if(t == FLOATERROR)
                t = getDomain()->getTime();
            int fdom = getInt("domain", arg);
            if (fdom == INTERROR)
                fdom = 0;
            int id = getInt("id", arg);
            if (id == INTERROR)
                id = 0;
            double fdepth = getFloat("fdepth", arg);
            if (fdepth == FLOATERROR)
                fdepth = currentSession.params->getDouble("initialFrontDepth");
            double kappa = getFloat("kappa", arg);
            if (kappa == FLOATERROR)
                kappa = 0.;
            string state = getString("state", arg);
            if (state == stringError or state == "moving")
                state = "init";
            if (numTabs != currentLevel + 1)
            {
                cout << getDomain()->getDomainID() << ": WARNING : asked for a FireNode "
                     << " with wrong indentation, treating it the current fire front" << endl;
            }

            if (state == "link")
            {

                /* Creating a link node */
                previousNode = getDomain()->addFireNode(pos, vel, t, fdepth, kappa, currentSession.ff, previousNode, fdom, id, FireNode::link);
                if (getDomain()->commandOutputs)
                    cout << getDomain()->getDomainID() << ": INIT -> added " << previousNode->toString() << endl;
            }
            else if (state == "final")
            {

                /* Creating a link node */
                previousNode = getDomain()->addFireNode(pos, vel, t, fdepth, kappa, currentSession.ff, previousNode, fdom, id, FireNode::final);
                if (getDomain()->commandOutputs)
                    cout << getDomain()->getDomainID() << ": INIT -> added " << previousNode->toString() << endl;
            }
            else if (getDomain()->striclyWithinDomain(pos) and getDomain()->striclyWithinDomain(*lastReadLoc))
            {

                /* Both nodes are within the domain and represent real firenodes */
                /* checking the distance between the two nodes */
                double distanceBetweenNodes = lastReadLoc->distance2D(pos);
                int interNodes = (int)floor(distanceBetweenNodes / (2. * perimRes));
                FFPoint posinc = (1. / (interNodes + 1)) * (pos - previousNode->getLoc());
                FFVector velinc = (1. / (interNodes + 1)) * (vel - previousNode->getVel());
                double timeinc = (1. / (interNodes + 1)) * (t - previousNode->getTime());
                FFPoint ipos;
                FFVector ivel;
                double itime;
                for (int k = 0; k < interNodes; k++)
                {
                    ipos = previousNode->getLoc() + posinc;
                    ivel = previousNode->getVel() + velinc;
                    itime = previousNode->getTime() + timeinc;
                    previousNode = getDomain()->addFireNode(ipos, ivel, itime, fdepth, kappa, currentSession.ff, previousNode);
                    if (getDomain()->commandOutputs)
                        cout << getDomain()->getDomainID() << ": INIT -> added inter-node "
                             << previousNode->toString() << endl;
                }
                // creating the firenode
                previousNode = getDomain()->addFireNode(pos, vel, t, fdepth, kappa, currentSession.ff, previousNode, fdom, id);
                if (getDomain()->commandOutputs)
                    cout << getDomain()->getDomainID()
                         << ": INIT -> added " << previousNode->toString() << endl;
            }
            else if (!getDomain()->striclyWithinDomain(pos) and getDomain()->striclyWithinDomain(*lastReadLoc))
            {

                FFPoint linkPoint = getDomain()->findIntersectionWithFrontiers(
                    pos, *lastReadLoc);
                rightLinkNode = getDomain()->addLinkNode(linkPoint);
                /* checking the distance between the two nodes */
                double distanceBetweenNodes = linkPoint.distance2D(previousNode->getLoc());
                if (distanceBetweenNodes > 2. * perimRes)
                {
                    int interNodes = (int)floor(distanceBetweenNodes / (2. * perimRes));
                    FFPoint posinc = (1. / (interNodes + 1)) * (linkPoint - previousNode->getLoc());
                    FFPoint ipos;
                    for (int k = 0; k < interNodes; k++)
                    {
                        ipos = previousNode->getLoc() + posinc;
                        previousNode = getDomain()->addFireNode(ipos, vel, t, fdepth, kappa, currentSession.ff, previousNode);
                        if (getDomain()->commandOutputs)
                            cout << getDomain()->getDomainID() << ": INIT -> added inter-node "
                                 << previousNode->toString() << endl;
                    }
                }
                previousNode->insertAfter(rightLinkNode);
                if (leftLinkNode != 0)
                    getDomain()->relateLinkNodes(rightLinkNode, leftLinkNode);
            }
            else if (getDomain()->striclyWithinDomain(pos) and !getDomain()->striclyWithinDomain(*lastReadLoc))
            {

                if (lastReadLoc->getX() == -numeric_limits<double>::infinity())
                {

                    /* First node to be created */
                    previousNode = getDomain()->addFireNode(pos, vel, t, fdepth, kappa, currentSession.ff, previousNode, fdom, id);
                    if (getDomain()->commandOutputs)
                        cout << getDomain()->getDomainID() << ": INIT -> added " << previousNode->toString() << endl;
                }
                else
                {

                    /* Creating a link node */
                    FFPoint linkPoint =
                        getDomain()->findIntersectionWithFrontiers(
                            pos, *lastReadLoc);
                    leftLinkNode = getDomain()->addLinkNode(linkPoint);
                    if (rightLinkNode != 0)
                    {
                        currentSession.ff->addFireNode(leftLinkNode, rightLinkNode);
                    }
                    else
                    {
                        leftLinkNode->setFront(currentSession.ff);
                        currentSession.ff->addFireNode(leftLinkNode);
                    }
                    previousNode = leftLinkNode;
                    /* checking the distance between the two nodes */
                    double distanceBetweenNodes = leftLinkNode->getLoc().distance2D(pos);
                    if (distanceBetweenNodes > 2. * perimRes)
                    {
                        int interNodes = (int)floor(distanceBetweenNodes / (2. * perimRes));
                        FFPoint posinc = (1. / (interNodes + 1)) * (pos - leftLinkNode->getLoc());
                        FFPoint ipos;
                        for (int k = 0; k < interNodes; k++)
                        {
                            ipos = previousNode->getLoc() + posinc;
                            previousNode = getDomain()->addFireNode(ipos, vel, t, fdepth, kappa, currentSession.ff, previousNode);
                            if (getDomain()->commandOutputs)
                                cout << getDomain()->getDomainID() << ": INIT -> added inter-node "
                                     << previousNode->toString() << endl;
                        }
                    }
                    // creating the firenode
                    previousNode = getDomain()->addFireNode(pos, vel, t, fdepth, kappa, currentSession.ff, previousNode, fdom, id);
                    if (getDomain()->commandOutputs)
                        cout << getDomain()->getDomainID()
                             << ": INIT -> added " << previousNode->toString() << endl;
                }
            }

            lastReadLoc->setX(pos.getX());
            lastReadLoc->setY(pos.getY());

            return normal;
        }
        else
        {
            throw MissingOption(3 - n);
        }
    }

    void Command::completeFront(FireFront *ff)
    {
        if (ff->getHead() == 0)
            return;
        FireNode *firstNode = ff->getHead();
        FireNode *lastNode = firstNode->getPrev();
        if (firstNode->getState() == FireNode::init and lastNode->getState() == FireNode::init)
        {
            double perimRes = getDomain()->getPerimeterResolution();
            double distanceBetweenNodes = lastNode->distance2D(firstNode);
            double fdepth = 0.5 * (firstNode->getFrontDepth() + lastNode->getFrontDepth());
            double kappa = 0.5 * (firstNode->getCurvature() + lastNode->getCurvature());
            int numNewNodes = (int)floor(distanceBetweenNodes / (2. * perimRes));
            if (numNewNodes > 0)
            {
                FFPoint posinc = (1. / (numNewNodes + 1)) * (firstNode->getLoc() - lastNode->getLoc());
                FFVector velinc = (1. / (numNewNodes + 1)) * (firstNode->getVel() - lastNode->getVel());
                double timeinc = (1. / (numNewNodes + 1)) * (firstNode->getTime() - lastNode->getTime());
                FireNode *prevNode = lastNode;
                FFPoint pos;
                FFVector vel;
                double t;
                for (int k = 0; k < numNewNodes; k++)
                {
                    pos = lastNode->getLoc() + (k + 1) * posinc;
                    vel = lastNode->getVel() + (k + 1) * velinc;
                    t = lastNode->getTime() + (k + 1) * timeinc;
                    prevNode = getDomain()->addFireNode(pos, vel, t, fdepth, kappa, currentSession.ff, prevNode);
                    prevNode->computeNormal();
                }
            }
        }
        if (getDomain()->commandOutputs)
        {
            cout << getDomain()->getDomainID() << ": "
                 << "****************************************" << endl;
            cout << getDomain()->getDomainID() << ": "
                 << " BEFORE THE BEGINNING OF THE SIMULATION: " << endl
                 << ff->print(1);
            cout << getDomain()->getDomainID() << ": "
                 << "****************************************" << endl;
        }

        if (!currentSession.params->isValued("BMapFiles"))
        {
            /* testing the domain for burning matrix */
            double iniFrontDepth = currentSession.params->getDouble("initialFrontDepth");
            double iniBurningTime = currentSession.params->getDouble("initialBurningDuration");

            if (!currentSession.params->isValued("noInitialScan"))
            {
                getDomain()->frontInitialBurningScan(ff->getTime(), ff, iniFrontDepth, iniBurningTime);
            }
        }

        currentFrontCompleted = true;
        delete lastReadLoc;
        lastReadLoc = 0;
    }

    int Command::stepSimulation(const string &arg, size_t &numTabs)
    {
        if (getDomain() == 0)
            return normal;

        double dt = getFloat("dt", arg);
        endTime = startTime + dt;
        ostringstream etime;
        etime.precision(numeric_limits<double>::digits10);
        etime << "t=" << endTime;
        goTo(etime.str().c_str(), numTabs);

        return normal;
    }

    int Command::goTo(const string &arg, size_t &numTabs)
    {

        /* Advancing simulation to the prescribed time */
        endTime = getFloat("t", arg); 
        if (init)
        {
            /* finishing the initialization process */
     
            bmapOutputUpdate = currentSession.params->getDouble("bmapOutputUpdate");
            if (!currentFrontCompleted)
                completeFront(currentSession.ff);
            numAtmoIterations = currentSession.params->getInt("numAtmoIterations") - 1;
            init = false;
          
        }

        if (endTime > startTime)
        {
            getDomain()->setTime(startTime);
            if (getDomain()->commandOutputs)
            {
                cout.precision(numeric_limits<double>::digits10);
                cout << getDomain()->getDomainID() << ": "
                     << "***************************************************" << endl;
                cout << getDomain()->getDomainID() << ": "
                     << "   ADVANCING FOREFIRE SIMULATION FROM T="
                     << startTime << " to " << endTime << endl;
                cout << getDomain()->getDomainID() << ": "