utils.cpp

Go to the documentation of this file.

 
/* ==============================================================================================================================
 
   This file is part of CoastalME, the Coastal Modelling Environment.
 
   CoastalME is free software; you can redistribute it and/or modify it under the terms of the GNU General Public  License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
 
   This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 
==============================================================================================================================*/
#include <assert.h>
 
#ifdef _WIN32
#include <windows.h>       // Needed for CalcProcessStats()
#include <psapi.h>
#include <io.h>            // For isatty()
#elif defined __GNUG__
#include <sys/resource.h>  // Needed for CalcProcessStats()
#include <unistd.h>        // For isatty()
#include <sys/types.h>
#include <sys/wait.h>
#endif
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
#include <stdio.h>
 
#include <cstdlib>
 
#include <cctype>
using std::tolower;
 
#include <cmath>
using std::floor;
 
#include <ctime>
using std::clock;
using std::clock_t;
using std::difftime;
using std::localtime;
using std::time;
 
#include <ios>
using std::fixed;
 
#include <iostream>
using std::cerr;
using std::cout;
using std::endl;
using std::ios;
 
#include <iomanip>
using std::put_time;
using std::resetiosflags;
using std::setprecision;
using std::setw;
 
#include <string>
using std::to_string;
 
#include <sstream>
using std::stringstream;
 
#include <algorithm>
using std::transform;
 
#include <gdal.h>
 
#include "cme.h"
#include "simulation.h"
#include "coast.h"
#include "2di_point.h"
 
//===============================================================================================================================
//===============================================================================================================================
int CSimulation::nHandleCommandLineParams(int nArg, char const* pcArgv[])
{
   if ((!isatty(fileno(stdout))) || (!isatty(fileno(stderr))))
      // Running with stdout or stderr not a tty, so either redirected or running as a background job. Ignore all command line parameters
      return RTN_OK;
 
   // Process the parameters following the name of the executable
   for (int i = 1; i < nArg; i++)
   {
      string strArg = pcArgv[i];
      strArg = strTrim(&strArg);
 
#ifdef _WIN32
      // Swap any forward slashes to backslashes
      strArg = pstrChangeToBackslash(&strArg);
#endif
 
      if (strArg.find("--gdal") != string::npos)
      {
         // User wants to know what GDAL raster drivers are available
         cout << GDAL_DRIVERS << endl
              << endl;
 
         for (int j = 0; j < GDALGetDriverCount(); j++)
         {
            GDALDriverH hDriver = GDALGetDriver(j);
 
            string strTmp(GDALGetDriverShortName(hDriver));
            strTmp.append("          ");
            strTmp.append(GDALGetDriverLongName(hDriver));
 
            cout << strTmp << endl;
         }
 
         return (RTN_HELP_ONLY);
      }
 
      else
      {
         if (strArg.find("--about") != string::npos)
         {
            // User wants information about CoastalME
            cout << ABOUT << endl;
            cout << THANKS << endl;
 
            return (RTN_HELP_ONLY);
         }
 
         else
         {
            if (strArg.find("--home") != string::npos)
            {
               // Read in user defined runtime directory
               // string strTmp;
 
               // Find the position of '='
               size_t const pos = strArg.find('=');
 
               // Was '=' found?
               if (pos != string::npos)
               {
                  // Yes, so get the substring after '=' and assign it to the global variable
                  m_strCMEIni = strArg.substr(pos + 1);
               }
 
               else
               {
                  // No
                  cout << "No '=' found in the input string" << endl;
               }
 
               return (RTN_OK);
            }
 
            // TODO 049 Handle other command line parameters e.g. path to .ini file, path to datafile
            else
            {
               // Display usage information
               cout << USAGE << endl;
               cout << USAGE1 << endl;
               cout << USAGE2 << endl;
               cout << USAGE3 << endl;
               cout << USAGE4 << endl;
               cout << USAGE5 << endl;
 
               return (RTN_HELP_ONLY);
            }
         }
      }
   }
 
   return RTN_OK;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceStart(void)
{
   cout << endl
        << PROGRAM_NAME << " for " << PLATFORM << " " << strGetBuild() << endl;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::StartClock(void)
{
   // First start the 'CPU time' clock ticking
   if (static_cast<clock_t>(-1) == clock())
   {
      // There's a problem with the clock, but continue anyway
      LogStream << NOTE << "CPU time not available" << endl;
      m_dCPUClock = -1;
   }
 
   else
   {
      // All OK, so get the time in m_dClkLast (this is needed to check for clock rollover on long runs)
      m_dClkLast = static_cast<double>(clock());
      m_dClkLast -= CLOCK_T_MIN; // necessary if clock_t is signed to make m_dClkLast unsigned
   }
 
   // And now get the actual time we started
   m_tSysStartTime = time(nullptr);
}
 
//===============================================================================================================================
//===============================================================================================================================
bool CSimulation::bFindExeDir(char const* pcArg)
{
   string strTmp;
   char szBuf[BUF_SIZE] = "";
 
#ifdef _WIN32
 
   if (0 != GetModuleFileName(NULL, szBuf, BUF_SIZE))
      strTmp = szBuf;
 
   else
      // It failed, so try another approach
      strTmp = pcArg;
 
#else
   // char* pResult = getcwd(szBuf, BUF_SIZE);          // Used to use this, but what if cwd is not the same as the CoastalME dir?
 
   if (-1 != readlink("/proc/self/exe", szBuf, BUF_SIZE))
      strTmp = szBuf;
 
   else
      // It failed, so try another approach
      strTmp = pcArg;
 
#endif
 
   // Neither approach has worked, so give up
   if (strTmp.empty())
      return false;
 
   // It's OK, so trim off the executable's name
   int const nPos = static_cast<int>(strTmp.find_last_of(PATH_SEPARATOR));
   m_strCMEDir = strTmp.substr(0, nPos + 1); // Note that this must be terminated with a backslash
 
   return true;
}
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceLicence(void)
{
   cout << COPYRIGHT << endl
        << endl;
   cout << LINE << endl;
   cout << DISCLAIMER1 << endl;
   cout << DISCLAIMER2 << endl;
   cout << DISCLAIMER3 << endl;
   cout << DISCLAIMER4 << endl;
   cout << DISCLAIMER5 << endl;
   cout << DISCLAIMER6 << endl;
   cout << LINE << endl
        << endl;
 
   cout << START_NOTICE << strGetComputerName() << " at " << put_time(localtime(&m_tSysStartTime), "%T on %A %d %B %Y") << endl;
   cout << INITIALIZING_NOTICE << endl;
}
 
//===============================================================================================================================
//===============================================================================================================================
double CSimulation::dGetTimeMultiplier(string const* strIn)
{
   // First decide what the time units are
   int const nTimeUnits = nDoTimeUnits(strIn);
 
   // Then return the correct multiplier, since m_dTimeStep is in hours
   switch (nTimeUnits)
   {
   case TIME_UNKNOWN:
      return TIME_UNKNOWN;
      break;
 
   case TIME_HOURS:
      return 1; // Multiplier for hours
      break;
 
   case TIME_DAYS:
      return 24; // Multiplier for days -> hours
      break;
 
   case TIME_MONTHS:
      return 24 * 30.416667; // Multiplier for months -> hours (assume 30 + 5/12 day months, no leap years)
      break;
 
   case TIME_YEARS:
      return 24 * 365.25; // Multiplier for years -> hours
      break;
   }
 
   return 0;
}
 
//===============================================================================================================================
//===============================================================================================================================
int CSimulation::nDoSimulationTimeMultiplier(string const* strIn)
{
   // First decide what the time units are
   int const nTimeUnits = nDoTimeUnits(strIn);
 
   // Next set up the correct multiplier, since m_dTimeStep is in hours
   switch (nTimeUnits)
   {
   case TIME_UNKNOWN:
      return RTN_ERR_TIMEUNITS;
      break;
 
   case TIME_HOURS:
      m_dDurationUnitsMult = 1; // Multiplier for hours
      m_strDurationUnits = "hours";
      break;
 
   case TIME_DAYS:
      m_dDurationUnitsMult = 24; // Multiplier for days -> hours
      m_strDurationUnits = "days";
      break;
 
   case TIME_MONTHS:
      m_dDurationUnitsMult = 24 * 30.416667; // Multiplier for months -> hours (assume 30 + 5/12 day months, no leap years)
      m_strDurationUnits = "months";
      break;
 
   case TIME_YEARS:
      m_dDurationUnitsMult = 24 * 365.25; // Multiplier for years -> hours
      m_strDurationUnits = "years";
      break;
   }
 
   return RTN_OK;
}
 
//===============================================================================================================================
//===============================================================================================================================
int CSimulation::nDoTimeUnits(string const* strIn)
{
   if (strIn->find("hour") != string::npos)
      return TIME_HOURS;
 
   else if (strIn->find("day") != string::npos)
      return TIME_DAYS;
 
   else if (strIn->find("month") != string::npos)
      return TIME_MONTHS;
 
   else if (strIn->find("year") != string::npos)
      return TIME_YEARS;
 
   else
      return TIME_UNKNOWN;
}
 
//===============================================================================================================================
//===============================================================================================================================
bool CSimulation::bOpenLogFile(void)
{
   if (m_nLogFileDetail == 0)
   {
      LogStream.open("/dev/null", ios::out | ios::trunc);
      cout << "Warning: log file is not writting" << endl;
   }
 
   else
      LogStream.open(m_strLogFile.c_str(), ios::out | ios::trunc);
 
   if (!LogStream)
   {
      // Error, cannot open log file
      cerr << ERR << "cannot open " << m_strLogFile << " for output" << endl;
      return false;
   }
 
   return true;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadBasementDEM(void) const
{
   // Tell the user what is happening
#ifdef _WIN32
   cout << READING_BASEMENT << pstrChangeToForwardSlash(&m_strInitialBasementDEMFile) << endl;
#else
   cout << READING_BASEMENT << m_strInitialBasementDEMFile << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceAllocateMemory(void)
{
   cout << ALLOCATE_MEMORY << endl;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceAddLayers(void)
{
   // Tell the user what is happening
   cout << ADD_LAYERS << endl;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadRasterFiles(void)
{
   cout << READING_RASTER_FILES << endl;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadVectorFiles(void)
{
   cout << READING_VECTOR_FILES << endl;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadLGIS(void) const
{
   // Tell the user what is happening
   if (! m_strInitialLandformFile.empty())
#ifdef _WIN32
      cout << READING_LANDFORM_FILE << pstrChangeToForwardSlash(&m_strInitialLandformFile) << endl;
 
#else
      cout << READING_LANDFORM_FILE << m_strInitialLandformFile << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadICGIS(void) const
{
   // Tell the user what is happening
   if (! m_strInterventionClassFile.empty())
#ifdef _WIN32
      cout << READING_INTERVENTION_CLASS_FILE << pstrChangeToForwardSlash(&m_strInterventionClassFile) << endl;
 
#else
      cout << READING_INTERVENTION_CLASS_FILE << m_strInterventionClassFile << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadIHGIS(void) const
{
   // Tell the user what is happening
   if (! m_strInterventionHeightFile.empty())
#ifdef _WIN32
      cout << READING_INTERVENTION_HEIGHT_FILE << pstrChangeToForwardSlash(&m_strInterventionHeightFile) << endl;
 
#else
      cout << READING_INTERVENTION_HEIGHT_FILE << m_strInterventionHeightFile << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadDeepWaterWaveValuesGIS(void) const
{
   // Tell the user what is happening
   if (! m_strDeepWaterWavesInputFile.empty())
#ifdef _WIN32
      cout << READING_DEEP_WATER_WAVE_FILE << pstrChangeToForwardSlash(&m_strDeepWaterWavesInputFile) << endl;
 
#else
      cout << READING_DEEP_WATER_WAVE_FILE << m_strDeepWaterWavesInputFile << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadSedimentEventInputValuesGIS(void) const
{
   // Tell the user what is happening
   if (! m_strSedimentInputEventFile.empty())
#ifdef _WIN32
      cout << READING_SED_INPUT_EVENT_FILE << pstrChangeToForwardSlash(&m_strSedimentInputEventFile) << endl;
 
#else
      cout << READING_SED_INPUT_EVENT_FILE << m_strSedimentInputEventFile << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadFloodLocationGIS(void) const
{
   // Tell the user what is happening
   if (! m_strFloodLocationShapefile.empty())
#ifdef _WIN32
      cout << READING_FLOOD_LOCATION << pstrChangeToForwardSlash(&m_strFloodLocationShapefile) << endl;
 
#else
      cout << READING_FLOOD_LOCATION << m_strFloodLocationShapefile << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadInitialSuspSedGIS(void) const
{
   // Tell the user what is happening
#ifdef _WIN32
   cout << READING_SUSPENDED_SEDIMENT_FILE << pstrChangeToForwardSlash(&m_strInitialSuspSedimentFile) << endl;
#else
   cout << READING_SUSPENDED_SEDIMENT_FILE << m_strInitialSuspSedimentFile << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadInitialFineUnconsSedGIS(int const nLayer) const
{
   // Tell the user what is happening
#ifdef _WIN32
   cout << READING_UNCONS_FINE_SEDIMENT_FILE << nLayer + 1 << "): " << pstrChangeToForwardSlash(&m_VstrInitialFineUnconsSedimentFile[nLayer]) << endl;
#else
   cout << READING_UNCONS_FINE_SEDIMENT_FILE << nLayer + 1 << "): " << m_VstrInitialFineUnconsSedimentFile[nLayer] << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadInitialSandUnconsSedGIS(int const nLayer) const
{
   // Tell the user what is happening
#ifdef _WIN32
   cout << READING_UNCONS_SAND_SEDIMENT_FILE << nLayer + 1 << "): " << pstrChangeToForwardSlash(&m_VstrInitialSandUnconsSedimentFile[nLayer]) << endl;
#else
   cout << READING_UNCONS_SAND_SEDIMENT_FILE << nLayer + 1 << "): " << m_VstrInitialSandUnconsSedimentFile[nLayer] << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadInitialCoarseUnconsSedGIS(int const nLayer) const
{
   // Tell the user what is happening
#ifdef _WIN32
   cout << READING_UNCONS_COARSE_SEDIMENT_FILE << nLayer + 1 << "): " << pstrChangeToForwardSlash(&m_VstrInitialCoarseUnconsSedimentFile[nLayer]) << endl;
#else
   cout << READING_UNCONS_COARSE_SEDIMENT_FILE << nLayer + 1 << "): " << m_VstrInitialCoarseUnconsSedimentFile[nLayer] << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadInitialFineConsSedGIS(int const nLayer) const
{
   // Tell the user what is happening
#ifdef _WIN32
   cout << READING_CONS_FINE_SEDIMENT_FILE << nLayer + 1 << "): " << pstrChangeToForwardSlash(&m_VstrInitialFineConsSedimentFile[nLayer]) << endl;
#else
   cout << READING_CONS_FINE_SEDIMENT_FILE << nLayer + 1 << "): " << m_VstrInitialFineConsSedimentFile[nLayer] << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadInitialSandConsSedGIS(int const nLayer) const
{
   // Tell the user what is happening
#ifdef _WIN32
   cout << READING_CONS_SAND_SEDIMENT_FILE << nLayer + 1 << "): " << pstrChangeToForwardSlash(&m_VstrInitialSandConsSedimentFile[nLayer]) << endl;
#else
   cout << READING_CONS_SAND_SEDIMENT_FILE << nLayer + 1 << "): " << m_VstrInitialSandConsSedimentFile[nLayer] << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadInitialCoarseConsSedGIS(int const nLayer) const
{
   // Tell the user what is happening
#ifdef _WIN32
   cout << READING_CONS_COARSE_SEDIMENT_FILE << nLayer + 1 << "): " << pstrChangeToForwardSlash(&m_VstrInitialCoarseConsSedimentFile[nLayer]) << endl;
#else
   cout << READING_CONS_COARSE_SEDIMENT_FILE << nLayer + 1 << "): " << m_VstrInitialCoarseConsSedimentFile[nLayer] << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadTideData(void) const
{
#ifdef _WIN32
   cout << READING_TIDE_DATA_FILE << pstrChangeToForwardSlash(&m_strTideDataFile) << endl;
#else
   cout << READING_TIDE_DATA_FILE << m_strTideDataFile << endl;
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceReadSCAPEShapeFunctionFile(void)
{
   cout << READING_SCAPE_SHAPE_FUNCTION_FILE << endl;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceInitializing(void)
{
   // Tell the user what is happening
   cout << INITIALIZING << endl;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceIsRunning(void)
{
   cout << RUN_NOTICE << endl;
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strListRasterFiles(void) const
{
   string strTmp;
 
   if (m_bBasementElevSave)
   {
      strTmp.append(RASTER_BASEMENT_ELEVATION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bSedimentTopSurfSave)
   {
      strTmp.append(RASTER_SEDIMENT_TOP_CODE);
      strTmp.append(", ");
   }
 
   if (m_bTopSurfSave)
   {
      strTmp.append(RASTER_TOP_CODE);
      strTmp.append(", ");
   }
 
   if (m_bSeaDepthSave)
   {
      strTmp.append(RASTER_SEA_DEPTH_NAME);
      strTmp.append(", ");
   }
 
   if (m_bAvgSeaDepthSave)
   {
      strTmp.append(RASTER_AVG_SEA_DEPTH_CODE);
      strTmp.append(", ");
   }
 
   if (m_bSeaMaskSave)
   {
      strTmp.append(RASTER_INUNDATION_MASK_CODE);
      strTmp.append(", ");
   }
 
   if (m_bWaveHeightSave)
   {
      strTmp.append(RASTER_WAVE_HEIGHT_CODE);
      strTmp.append(", ");
   }
 
   if (m_bWaveAngleSave)
   {
      strTmp.append(RASTER_WAVE_ORIENTATION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bAvgWaveHeightSave)
   {
      strTmp.append(RASTER_AVG_WAVE_HEIGHT_CODE);
      strTmp.append(", ");
   }
 
   if (m_bBeachProtectionSave)
   {
      strTmp.append(RASTER_BEACH_PROTECTION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bPotentialPlatformErosionSave)
   {
      strTmp.append(RASTER_POTENTIAL_PLATFORM_EROSION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bPotentialPlatformErosionMaskSave)
   {
      strTmp.append(RASTER_POTENTIAL_PLATFORM_EROSION_MASK_CODE);
      strTmp.append(", ");
   }
 
   if (m_bBeachDepositionSave)
   {
      strTmp.append(RASTER_BEACH_DEPOSITION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bTotalBeachDepositionSave)
   {
      strTmp.append(RASTER_TOTAL_BEACH_DEPOSITION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bBeachMaskSave)
   {
      strTmp.append(RASTER_BEACH_MASK_CODE);
      strTmp.append(", ");
   }
 
   if (m_bActualPlatformErosionSave)
   {
      strTmp.append(RASTER_ACTUAL_PLATFORM_EROSION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bTotalPotentialPlatformErosionSave)
   {
      strTmp.append(RASTER_TOTAL_POTENTIAL_PLATFORM_EROSION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bTotalActualPlatformErosionSave)
   {
      strTmp.append(RASTER_TOTAL_ACTUAL_PLATFORM_EROSION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bPotentialBeachErosionSave)
   {
      strTmp.append(RASTER_POTENTIAL_BEACH_EROSION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bActualBeachErosionSave)
   {
      strTmp.append(RASTER_ACTUAL_BEACH_EROSION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bTotalPotentialBeachErosionSave)
   {
      strTmp.append(RASTER_TOTAL_POTENTIAL_BEACH_EROSION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bTotalActualBeachErosionSave)
   {
      strTmp.append(RASTER_TOTAL_ACTUAL_BEACH_EROSION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bLandformSave)
   {
      strTmp.append(RASTER_LANDFORM_CODE);
      strTmp.append(", ");
   }
 
   if (m_bLocalSlopeSave)
   {
      strTmp.append(RASTER_LOCAL_SLOPE_CODE);
      strTmp.append(", ");
   }
 
   if (m_bInterventionClassSave)
   {
      strTmp.append(RASTER_INTERVENTION_CLASS_CODE);
      strTmp.append(", ");
   }
 
   if (m_bInterventionHeightSave)
   {
      strTmp.append(RASTER_INTERVENTION_HEIGHT_CODE);
      strTmp.append(", ");
   }
 
   if (m_bHaveFineSediment && m_bSuspSedSave)
   {
      strTmp.append(RASTER_SUSP_SED_CODE);
      strTmp.append(", ");
   }
 
   if (m_bHaveFineSediment && m_bAvgSuspSedSave)
   {
      strTmp.append(RASTER_AVG_SUSP_SED_CODE);
      strTmp.append(", ");
   }
 
   if (m_bHaveFineSediment && m_bFineUnconsSedSave)
   {
      strTmp.append(RASTER_FINE_UNCONS_CODE);
      strTmp.append(", ");
   }
 
   if (m_bHaveSandSediment && m_bSandUnconsSedSave)
   {
      strTmp.append(RASTER_SAND_UNCONS_CODE);
      strTmp.append(", ");
   }
 
   if (m_bHaveCoarseSediment && m_bCoarseUnconsSedSave)
   {
      strTmp.append(RASTER_COARSE_UNCONS_CODE);
      strTmp.append(", ");
   }
 
   if (m_bHaveFineSediment && m_bFineConsSedSave)
   {
      strTmp.append(RASTER_FINE_CONS_CODE);
      strTmp.append(", ");
   }
 
   if (m_bHaveSandSediment && m_bSandConsSedSave)
   {
      strTmp.append(RASTER_SAND_CONS_CODE);
      strTmp.append(", ");
   }
 
   if (m_bHaveCoarseSediment && m_bCoarseConsSedSave)
   {
      strTmp.append(RASTER_COARSE_CONS_CODE);
      strTmp.append(", ");
   }
 
   if (m_bRasterCoastlineSave)
   {
      strTmp.append(RASTER_COAST_CODE);
      strTmp.append(", ");
   }
 
   if (m_bRasterNormalProfileSave)
   {
      strTmp.append(RASTER_COAST_NORMAL_CODE);
      strTmp.append(", ");
   }
 
   if (m_bActiveZoneSave)
   {
      strTmp.append(RASTER_ACTIVE_ZONE_CODE);
      strTmp.append(", ");
   }
 
   if (m_bRasterPolygonSave)
   {
      strTmp.append(RASTER_POLYGON_CODE);
      strTmp.append(", ");
   }
 
   if (m_bPotentialPlatformErosionMaskSave)
   {
      strTmp.append(RASTER_POTENTIAL_PLATFORM_EROSION_MASK_CODE);
      strTmp.append(", ");
   }
 
   if (m_bSedimentInputEventSave)
   {
      strTmp.append(RASTER_SEDIMENT_INPUT_EVENT_CODE);
      strTmp.append(", ");
   }
 
   // Remove the trailing comma and space
   if (strTmp.size() > 2)
      strTmp.resize(strTmp.size() - 2);
 
   return strTmp;
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strListVectorFiles(void) const
{
   string strTmp;
 
   if (m_bCoastSave)
   {
      strTmp.append(VECTOR_COAST_CODE);
      strTmp.append(", ");
   }
 
   if (m_bNormalsSave)
   {
      strTmp.append(VECTOR_NORMALS_CODE);
      strTmp.append(", ");
   }
 
   if (m_bInvalidNormalsSave)
   {
      strTmp.append(VECTOR_INVALID_NORMALS_CODE);
      strTmp.append(", ");
   }
 
   if (m_bWaveAngleAndHeightSave)
   {
      strTmp.append(VECTOR_WAVE_ANGLE_AND_HEIGHT_CODE);
      strTmp.append(", ");
   }
 
   if (m_bAvgWaveAngleAndHeightSave)
   {
      strTmp.append(VECTOR_AVG_WAVE_ANGLE_AND_HEIGHT_CODE);
      strTmp.append(", ");
   }
 
   if (m_bCoastCurvatureSave)
   {
      strTmp.append(VECTOR_COAST_CURVATURE_CODE);
      strTmp.append(", ");
   }
 
   if (m_bWaveEnergySinceCollapseSave)
   {
      strTmp.append(VECTOR_WAVE_ENERGY_SINCE_COLLAPSE_CODE);
      strTmp.append(", ");
   }
 
   if (m_bMeanWaveEnergySave)
   {
      strTmp.append(VECTOR_MEAN_WAVE_ENERGY_CODE);
      strTmp.append(", ");
   }
 
   if (m_bBreakingWaveHeightSave)
   {
      strTmp.append(VECTOR_BREAKING_WAVE_HEIGHT_CODE);
      strTmp.append(", ");
   }
 
   if (m_bPolygonNodeSave)
   {
      strTmp.append(VECTOR_POLYGON_NODE_CODE);
      strTmp.append(", ");
   }
 
   if (m_bPolygonBoundarySave)
   {
      strTmp.append(VECTOR_POLYGON_BOUNDARY_CODE);
      strTmp.append(", ");
   }
 
   if (m_bCliffNotchSave)
   {
      strTmp.append(VECTOR_CLIFF_NOTCH_SIZE_CODE);
      strTmp.append(", ");
   }
 
   if (m_bShadowBoundarySave)
   {
      strTmp.append(VECTOR_SHADOW_BOUNDARY_CODE);
      strTmp.append(", ");
   }
 
   if (m_bShadowDowndriftBoundarySave)
   {
      strTmp.append(VECTOR_DOWNDRIFT_BOUNDARY_CODE);
      strTmp.append(", ");
   }
 
   if (m_bDeepWaterWaveAngleAndHeightSave)
   {
      strTmp.append(VECTOR_DEEP_WATER_WAVE_ANGLE_AND_HEIGHT_CODE);
      strTmp.append(", ");
   }
 
   // Remove the trailing comma and space
   if (strTmp.size() > 2)
      strTmp.resize(strTmp.size() - 2);
 
   return strTmp;
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strListTSFiles(void) const
{
   string strTmp;
 
   if (m_bSeaAreaTSSave)
   {
      strTmp.append(TIME_SERIES_SEA_AREA_CODE);
      strTmp.append(", ");
   }
 
   if (m_bStillWaterLevelTSSave)
   {
      strTmp.append(TIME_SERIES_STILL_WATER_LEVEL_CODE);
      strTmp.append(", ");
   }
 
   if (m_bActualPlatformErosionTSSave)
   {
      strTmp.append(TIME_SERIES_PLATFORM_EROSION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bCliffCollapseErosionTSSave)
   {
      strTmp.append(TIME_SERIES_CLIFF_COLLAPSE_EROSION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bCliffCollapseDepositionTSSave)
   {
      strTmp.append(TIME_SERIES_CLIFF_COLLAPSE_DEPOSITION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bCliffCollapseNetTSSave)
   {
      strTmp.append(TIME_SERIES_CLIFF_COLLAPSE_NET_CODE);
      strTmp.append(", ");
   }
 
   if (m_bBeachErosionTSSave)
   {
      strTmp.append(TIME_SERIES_BEACH_EROSION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bBeachDepositionTSSave)
   {
      strTmp.append(TIME_SERIES_BEACH_DEPOSITION_CODE);
      strTmp.append(", ");
   }
 
   if (m_bBeachSedimentChangeNetTSSave)
   {
      strTmp.append(TIME_SERIES_BEACH_CHANGE_NET_CODE);
      strTmp.append(", ");
   }
 
   if (m_bSuspSedTSSave)
   {
      strTmp.append(TIME_SERIES_SUSPENDED_SEDIMENT_CODE);
      strTmp.append(", ");
   }
 
   if (m_bFloodSetupSurgeTSSave)
   {
      strTmp.append(TIME_SERIES_FLOOD_SETUP_SURGE_CODE);
      strTmp.append(", ");
   }
 
   if (m_bFloodSetupSurgeRunupTSSave)
   {
      strTmp.append(TIME_SERIES_FLOOD_SETUP_SURGE_RUNUP_CODE);
      strTmp.append(", ");
   }
 
   // Remove the trailing comma and space
   if (strTmp.size() > 2)
      strTmp.resize(strTmp.size() - 2);
 
   return strTmp;
}
 
//===============================================================================================================================
//===============================================================================================================================
bool CSimulation::bSetUpTSFiles(void)
{
   string strTSFile;
 
   if (m_bSeaAreaTSSave)
   {
      // Start with wetted area
      strTSFile = m_strOutPath;
      strTSFile.append(TIME_SERIES_SEA_AREA_NAME);
      strTSFile.append(CSVEXT);
 
      // Open wetted time-series CSV file
      SeaAreaTSStream.open(strTSFile.c_str(), ios::out | ios::trunc);
 
      if (!SeaAreaTSStream)
      {
         // Error, cannot open wetted area  time-series file
         cerr << ERR << "cannot open " << strTSFile << " for output" << endl;
         return false;
      }
   }
 
   if (m_bStillWaterLevelTSSave)
   {
      // Now still water level
      strTSFile = m_strOutPath;
      strTSFile.append(TIME_SERIES_STILL_WATER_LEVEL_NAME);
      strTSFile.append(CSVEXT);
 
      // Open still water level time-series CSV file
      StillWaterLevelTSStream.open(strTSFile.c_str(), ios::out | ios::trunc);
 
      if (!StillWaterLevelTSStream)
      {
         // Error, cannot open still water level time-series file
         cerr << ERR << "cannot open " << strTSFile << " for output" << endl;
         return false;
      }
   }
 
   if (m_bActualPlatformErosionTSSave)
   {
      // Erosion (fine, sand, coarse)
      strTSFile = m_strOutPath;
      strTSFile.append(TIME_SERIES_PLATFORM_EROSION_NAME);
      strTSFile.append(CSVEXT);
 
      // Open erosion time-series CSV file
      PlatformErosionTSStream.open(strTSFile.c_str(), ios::out | ios::trunc);
 
      if (! PlatformErosionTSStream)
      {
         // Error, cannot open erosion time-series file
         cerr << ERR << "cannot open " << strTSFile << " for output" << endl;
         return false;
      }
   }
 
   if (m_bCliffCollapseErosionTSSave)
   {
      // Erosion due to cliff collapse
      strTSFile = m_strOutPath;
      strTSFile.append(TIME_SERIES_CLIFF_COLLAPSE_EROSION_NAME);
      strTSFile.append(CSVEXT);
 
      // Open cliff collapse erosion time-series CSV file
      CliffCollapseErosionTSStream.open(strTSFile.c_str(), ios::out | ios::trunc);
 
      if (!CliffCollapseErosionTSStream)
      {
         // Error, cannot open cliff collapse erosion time-series file
         cerr << ERR << "cannot open " << strTSFile << " for output" << endl;
         return false;
      }
   }
 
   if (m_bCliffCollapseDepositionTSSave)
   {
      // Deposition due to cliff collapse
      strTSFile = m_strOutPath;
      strTSFile.append(TIME_SERIES_CLIFF_COLLAPSE_DEPOSITION_NAME);
      strTSFile.append(CSVEXT);
 
      // Open cliff collapse deposition time-series CSV file
      CliffCollapseDepositionTSStream.open(strTSFile.c_str(), ios::out | ios::trunc);
 
      if (!CliffCollapseDepositionTSStream)
      {
         // Error, cannot open cliff collapse deposition time-series file
         cerr << ERR << "cannot open " << strTSFile << " for output" << endl;
         return false;
      }
   }
 
   if (m_bCliffCollapseNetTSSave)
   {
      // Net change in unconsolidated sediment due to cliff collapse
      strTSFile = m_strOutPath;
      strTSFile.append(TIME_SERIES_CLIFF_COLLAPSE_NET_NAME);
      strTSFile.append(CSVEXT);
 
      // Open net cliff collapse time-series CSV file
      CliffCollapseNetChangeTSStream.open(strTSFile.c_str(), ios::out | ios::trunc);
 
      if (!CliffCollapseNetChangeTSStream)
      {
         // Error, cannot open net cliff collapse time-series file
         cerr << ERR << "cannot open " << strTSFile << " for output" << endl;
         return false;
      }
   }
 
   if (m_bBeachErosionTSSave)
   {
      // Beach erosion
      strTSFile = m_strOutPath;
      strTSFile.append(TIME_SERIES_BEACH_EROSION_NAME);
      strTSFile.append(CSVEXT);
 
      // Open beach erosion time-series CSV file
      BeachErosionTSStream.open(strTSFile.c_str(), ios::out | ios::trunc);
 
      if (!BeachErosionTSStream)
      {
         // Error, cannot open beach erosion time-series file
         cerr << ERR << "cannot open " << strTSFile << " for output" << endl;
         return false;
      }
   }
 
   if (m_bBeachDepositionTSSave)
   {
      // Beach deposition
      strTSFile = m_strOutPath;
      strTSFile.append(TIME_SERIES_BEACH_DEPOSITION_NAME);
      strTSFile.append(CSVEXT);
 
      // Open beach deposition time-series CSV file
      BeachDepositionTSStream.open(strTSFile.c_str(), ios::out | ios::trunc);
 
      if (!BeachDepositionTSStream)
      {
         // Error, cannot open beach deposition time-series file
         cerr << ERR << "cannot open " << strTSFile << " for output" << endl;
         return false;
      }
   }
 
   if (m_bBeachSedimentChangeNetTSSave)
   {
      // Beach sediment change
      strTSFile = m_strOutPath;
      strTSFile.append(TIME_SERIES_BEACH_CHANGE_NET_NAME);
      strTSFile.append(CSVEXT);
 
      // Open net beach sediment change time-series CSV file
      BeachSedimentNetChangeTSStream.open(strTSFile.c_str(), ios::out | ios::trunc);
 
      if (!BeachSedimentNetChangeTSStream)
      {
         // Error, cannot open beach sediment change time-series file
         cerr << ERR << "cannot open " << strTSFile << " for output" << endl;
         return false;
      }
   }
 
   if (m_bSuspSedTSSave)
   {
      // Sediment load
      strTSFile = m_strOutPath;
      strTSFile.append(TIME_SERIES_SUSPENDED_SEDIMENT_NAME);
      strTSFile.append(CSVEXT);
 
      // Open sediment load time-series CSV file
      FineSedSuspensionTSStream.open(strTSFile.c_str(), ios::out | ios::trunc);
 
      if (!FineSedSuspensionTSStream)
      {
         // Error, cannot open sediment load time-series file
         cerr << ERR << "cannot open " << strTSFile << " for output" << endl;
         return false;
      }
   }
 
   if (m_bFloodSetupSurgeTSSave)
   {
      // Sediment load
      strTSFile = m_strOutPath;
      strTSFile.append(TIME_SERIES_FLOOD_SETUP_SURGE_CODE);
      strTSFile.append(CSVEXT);
 
      // Open sediment load time-series CSV file
      FloodSetupSurgeTSStream.open(strTSFile.c_str(), ios::out | ios::trunc);
 
      if (!FloodSetupSurgeTSStream)
      {
         // Error, cannot open sediment load time-series file
         cerr << ERR << "cannot open " << strTSFile << " for output" << endl;
         return false;
      }
   }
 
   if (m_bFloodSetupSurgeRunupTSSave)
   {
      // Sediment load
      strTSFile = m_strOutPath;
      strTSFile.append(TIME_SERIES_FLOOD_SETUP_SURGE_RUNUP_CODE);
      strTSFile.append(CSVEXT);
 
      // Open sediment load time-series CSV file
      FloodSetupSurgeRunupTSStream.open(strTSFile.c_str(), ios::out | ios::trunc);
 
      if (!FloodSetupSurgeRunupTSStream)
      {
         // Error, cannot open sediment load time-series file
         cerr << ERR << "cannot open " << strTSFile << " for output" << endl;
         return false;
      }
   }
 
   return true;
}
 
//===============================================================================================================================
//===============================================================================================================================
bool CSimulation::bTimeToQuit(void)
{
   // Add timestep to the total time simulated so far
   m_dSimElapsed += m_dTimeStep;
 
   if (m_dSimElapsed >= m_dSimDuration)
   {
      // It is time to quit
      m_dSimElapsed = m_dSimDuration;
      AnnounceProgress();
      return true;
   }
 
   // Not quitting, so increment the timestep count, and recalc total timesteps
   m_ulIter++;
   m_ulTotTimestep = static_cast<unsigned long>(dRound(m_dSimDuration / m_dTimeStep));
 
   // Check to see if we have done CLOCK_CHECK_ITERATION timesteps: if so, it is time to reset the CPU time running total in case the clock() function later rolls over
   if (0 == m_ulIter % CLOCK_CHECK_ITERATION)
      DoCPUClockReset();
 
   // Not yet time to quit
   return false;
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strGetComputerName(void)
{
   string strComputerName;
 
#ifdef _WIN32
   // Being compiled to run under Windows, either by MS VC++, Borland C++, or Cygwin
   strComputerName = getenv("COMPUTERNAME");
#else
   // Being compiled for another platform; assume for Linux-Unix
   char szHostName[BUF_SIZE] = "";
   gethostname(szHostName, BUF_SIZE);
 
   strComputerName = szHostName;
 
   if (strComputerName.empty())
      strComputerName = "Unknown Computer";
 
#endif
 
   return strComputerName;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::DoCPUClockReset(void)
{
   if (static_cast<clock_t>(-1) == clock())
   {
      // Error
      LogStream << "CPU time not available" << endl;
      m_dCPUClock = -1;
      return;
   }
 
   // OK, so carry on
   double dClkThis = static_cast<double>(clock());
   dClkThis -= CLOCK_T_MIN; // necessary when clock_t is signed, to make dClkThis unsigned
 
   if (dClkThis < m_dClkLast)
   {
      // Clock has 'rolled over'
      m_dCPUClock += (CLOCK_T_RANGE + 1 - m_dClkLast); // this elapsed before rollover
      m_dCPUClock += dClkThis;                         // this elapsed after rollover
 
#ifdef CLOCKCHECK
      // For debug purposes
      LogStream << "Rolled over: dClkThis=" << dClkThis << " m_dClkLast=" << m_dClkLast << endl
                << "\t"
                << " before rollover=" << (CLOCK_T_RANGE + 1 - m_dClkLast) << endl
                << "\t"
                << " after rollover=" << dClkThis << endl
                << "\t"
                << " ADDED=" << (CLOCK_T_RANGE + 1 - m_dClkLast + dClkThis) << endl;
#endif
   }
 
   else
   {
      // No rollover
      m_dCPUClock += (dClkThis - m_dClkLast);
 
#ifdef CLOCKCHECK
      // For debug purposes
      LogStream << "No rollover: dClkThis=" << dClkThis << " m_dClkLast=" << m_dClkLast << " ADDED=" << dClkThis - m_dClkLast << endl;
#endif
   }
 
   // Reset for next time
   m_dClkLast = dClkThis;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceSimEnd(void)
{
   cout << endl
        << FINAL_OUTPUT << endl;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::CalcTime(double const dRunLength)
{
   // Reset CPU count for last time
   DoCPUClockReset();
 
   if (! bFPIsEqual(m_dCPUClock, -1.0, TOLERANCE))
   {
      // Calculate CPU time in secs
      double const dDuration = m_dCPUClock / CLOCKS_PER_SEC;
 
      // And write CPU time out to OutStream and LogStream
      OutStream << "CPU time elapsed: " << strDispTime(dDuration, false, true);
      LogStream << "CPU time elapsed: " << strDispTime(dDuration, false, true);
 
      // Calculate CPU time per timestep
      double const dPerTimestep = dDuration / static_cast<double>(m_ulTotTimestep);
 
      // And write CPU time per timestep to OutStream and LogStream
      OutStream << fixed << setprecision(4) << " (" << dPerTimestep << " per timestep)" << endl;
      LogStream << fixed << setprecision(4) << " (" << dPerTimestep << " per timestep)" << endl;
 
      // Calculate ratio of CPU time to time simulated
      OutStream << resetiosflags(ios::floatfield);
      OutStream << fixed << setprecision(0) << "In terms of CPU time, this is ";
      LogStream << resetiosflags(ios::floatfield);
      LogStream << fixed << setprecision(0) << "In terms of CPU time, this is ";
 
      if (dDuration > dRunLength)
      {
         OutStream << dDuration / dRunLength << " x slower than reality" << endl;
         LogStream << dDuration / dRunLength << " x slower than reality" << endl;
      }
 
      else
      {
         OutStream << dRunLength / dDuration << " x faster than reality" << endl;
         LogStream << dRunLength / dDuration << " x faster than reality" << endl;
      }
   }
 
   // Calculate run time
   double const dDuration = difftime(m_tSysEndTime, m_tSysStartTime);
 
   // And write run time out to OutStream and LogStream
   OutStream << "Run time elapsed: " << strDispTime(dDuration, false, false);
   LogStream << "Run time elapsed: " << strDispTime(dDuration, false, false);
 
   // Calculate run time per timestep
   double const dPerTimestep = dDuration / static_cast<double>(m_ulTotTimestep);
 
   // And write run time per timestep to OutStream and LogStream
   OutStream << resetiosflags(ios::floatfield);
   OutStream << " (" << fixed << setprecision(4) << dPerTimestep << " per timestep)" << endl;
   LogStream << resetiosflags(ios::floatfield);
   LogStream << " (" << fixed << setprecision(4) << dPerTimestep << " per timestep)" << endl;
 
   // Calculate ratio of run time to time simulated
   OutStream << fixed << setprecision(0) << "In terms of run time, this is ";
   LogStream << fixed << setprecision(0) << "In terms of run time, this is ";
 
   if (dDuration > dRunLength)
   {
      OutStream << dDuration / dRunLength << " x slower than reality" << endl;
      LogStream << dDuration / dRunLength << " x slower than reality" << endl;
   }
 
   else
   {
      OutStream << dRunLength / dDuration << " x faster than reality" << endl;
      LogStream << dRunLength / dDuration << " x faster than reality" << endl;
   }
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strDispSimTime(const double dTimeIn)
{
   // Make sure no negative times
   double dTmpTime = tMax(dTimeIn, 0.0);
 
   string strTime;
 
   // Constants
   double const dHoursInYear = 24 * 365; // it was 365.25
   double const dHoursInDay = 24;
 
   // Display years
   if (dTmpTime >= dHoursInYear)
   {
      double const dYears = floor(dTmpTime / dHoursInYear);
      dTmpTime -= (dYears * dHoursInYear);
 
      strTime = to_string(static_cast<int>(dYears));
      strTime.append("y ");
   }
 
   else
      strTime = "0y ";
 
   // Display Julian days
   if (dTmpTime >= dHoursInDay)
   {
      double const dJDays = floor(dTmpTime / dHoursInDay);
      dTmpTime -= (dJDays * dHoursInDay);
 
      stringstream ststrTmp;
      ststrTmp << FillToWidth('0', 3) << static_cast<int>(dJDays);
      strTime.append(ststrTmp.str());
      strTime.append("d ");
   }
 
   else
      strTime.append("000d ");
 
   // Display hours
   stringstream ststrTmp;
   ststrTmp << FillToWidth('0', 2) << static_cast<int>(dTmpTime);
   strTime.append(ststrTmp.str());
   strTime.append("h");
 
   return strTime;
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strDispTime(const double dTimeIn, const bool bRound, const bool bFrac)
{
   // Make sure no negative times
   double dTime = tMax(dTimeIn, 0.0);
 
   string strTime;
 
   if (bRound)
      dTime = dRound(dTime);
 
   unsigned long ulTimeIn = static_cast<unsigned long>(floor(dTime));
   dTime -= static_cast<double>(ulTimeIn);
 
   // Hours
   if (ulTimeIn >= 3600)
   {
      // Display some hours
      unsigned long const ulHours = ulTimeIn / 3600ul;
      ulTimeIn -= (ulHours * 3600ul);
 
      strTime = to_string(ulHours);
      strTime.append(":");
   }
 
   else
      strTime = "0:";
 
   // Minutes
   if (ulTimeIn >= 60)
   {
      // display some minutes
      unsigned long const ulMins = ulTimeIn / 60ul;
      ulTimeIn -= (ulMins * 60ul);
 
      stringstream ststrTmp;
      ststrTmp << FillToWidth('0', 2) << ulMins;
      strTime.append(ststrTmp.str());
      strTime.append(":");
   }
 
   else
      strTime.append("00:");
 
   // Seconds
   stringstream ststrTmp;
   ststrTmp << FillToWidth('0', 2) << ulTimeIn;
   strTime.append(ststrTmp.str());
 
   if (bFrac)
   {
      // Fractions of a second
      strTime.append(".");
      ststrTmp.clear();
      ststrTmp.str(string());
      ststrTmp << FillToWidth('0', 2) << static_cast<unsigned long>(dTime * 100);
      strTime.append(ststrTmp.str());
   }
 
   return strTime;
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strGetBuild(void)
{
   string strBuild("(");
   strBuild.append(__TIME__);
   strBuild.append(" ");
   strBuild.append(__DATE__);
#ifdef _DEBUG
   strBuild.append(" DEBUG");
#endif
   strBuild.append(" build)");
 
   return strBuild;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AnnounceProgress(void)
{
   if (isatty(fileno(stdout)))
   {
      // Stdout is connected to a tty, so not running as a background job
      static double sdElapsed = 0;
      static double sdToGo = 0;
      time_t const tNow = time(nullptr);
 
      // Calculate time elapsed and remaining
      sdElapsed = difftime(tNow, m_tSysStartTime);
      sdToGo = (sdElapsed * m_dSimDuration / m_dSimElapsed) - sdElapsed;
 
      // Tell the user about progress (note need to make several separate calls to cout here, or MS VC++ compiler appears to get confused)
      cout << SIMULATING << strDispSimTime(m_dSimElapsed);
      cout << fixed << setprecision(3) << setw(9) << 100 * m_dSimElapsed / m_dSimDuration;
      cout << "%   (elapsed " << strDispTime(sdElapsed, false, false) << " remaining ";
 
      cout << strDispTime(sdToGo, false, false) << ")  ";
 
      // Add a 'marker' for GIS saves etc.
      if (m_bSaveGISThisIter)
         cout << setw(9) << "GIS" + to_string(m_nGISSave);
 
      else if (m_bSedimentInputThisIter)
         cout << setw(9) << "SED INPUT";
 
      else
         cout << setw(9) << SPACE;
 
      cout.flush();
   }
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::CalcProcessStats(void)
{
   string const NA = "Not available";
 
   OutStream << endl;
   OutStream << "Process statistics" << endl;
   OutStream << "------------------" << endl;
 
#ifdef _WIN32
   // First, find out which version of Windows we are running under
   OSVERSIONINFOEX osvi;
   BOOL bOsVersionInfoEx;
 
   ZeroMemory(&osvi, sizeof(OSVERSIONINFOEX)); // fill this much memory with zeros
   osvi.dwOSVersionInfoSize = sizeof(OSVERSIONINFOEX);
 
   if (!(bOsVersionInfoEx = GetVersionEx((OSVERSIONINFO*)&osvi)))
   {
      // OSVERSIONINFOEX didn't work so try OSVERSIONINFO instead
      osvi.dwOSVersionInfoSize = sizeof(OSVERSIONINFO);
 
      if (!GetVersionEx((OSVERSIONINFO*)&osvi))
      {
         // That didn't work either, too risky to proceed so give up
         OutStream << NA << endl;
         return;
      }
   }
 
   // OK, we have Windows version so display it
   OutStream << "Running under                                \t: ";
 
   switch (osvi.dwPlatformId)
   {
   case VER_PLATFORM_WIN32_NT:
      if (osvi.dwMajorVersion <= 4)
         OutStream << "Windows NT ";
 
      else if (5 == osvi.dwMajorVersion && 0 == osvi.dwMinorVersion)
         OutStream << "Windows 2000 ";
 
      else if (5 == osvi.dwMajorVersion && 1 == osvi.dwMinorVersion)
         OutStream << "Windows XP ";
 
      else if (6 == osvi.dwMajorVersion && 0 == osvi.dwMinorVersion)
         OutStream << "Windows Vista ";
 
      else if (6 == osvi.dwMajorVersion && 1 == osvi.dwMinorVersion)
         OutStream << "Windows 7 ";
 
      else if (6 == osvi.dwMajorVersion && 2 == osvi.dwMinorVersion)
         OutStream << "Windows 8 ";
 
      else if (6 == osvi.dwMajorVersion && 3 == osvi.dwMinorVersion)
         OutStream << "Windows 8.1 ";
 
      else if (10 == osvi.dwMajorVersion && 0 == osvi.dwMinorVersion)
         OutStream << "Windows 10 ";
 
      else if (11 == osvi.dwMajorVersion && 0 == osvi.dwMinorVersion)
         OutStream << "Windows 11 ";
 
      else
         OutStream << "unknown Windows version ";
 
      // Display version, service pack (if any), and build number
      if (osvi.dwMajorVersion <= 4)
         OutStream << "version " << osvi.dwMajorVersion << "." << osvi.dwMinorVersion << " " << osvi.szCSDVersion << " (Build " << (osvi.dwBuildNumber & 0xFFFF) << ")" << endl;
 
      else
         OutStream << osvi.szCSDVersion << " (Build " << (osvi.dwBuildNumber & 0xFFFF) << ")" << endl;
 
      break;
 
   case VER_PLATFORM_WIN32_WINDOWS:
      if (4 == osvi.dwMajorVersion && 0 == osvi.dwMinorVersion)
      {
         OutStream << "Windows 95";
 
         if ('C' == osvi.szCSDVersion[1] || 'B' == osvi.szCSDVersion[1])
            OutStream << " OSR2";
 
         OutStream << endl;
      }
 
      else if (4 == osvi.dwMajorVersion && 10 == osvi.dwMinorVersion)
      {
         OutStream << "Windows 98";
 
         if ('A' == osvi.szCSDVersion[1])
            OutStream << "SE";
 
         OutStream << endl;
      }
 
      else if (4 == osvi.dwMajorVersion && 90 == osvi.dwMinorVersion)
         OutStream << "Windows Me" << endl;
 
      else
         OutStream << "unknown 16-bit Windows version " << endl;
 
      break;
 
   case VER_PLATFORM_WIN32s:
      OutStream << "Win32s" << endl;
      break;
   }
 
   // Now get process timimgs: this only works under 32-bit windows
   if (VER_PLATFORM_WIN32_NT == osvi.dwPlatformId)
   {
      FILETIME ftCreate, ftExit, ftKernel, ftUser;
 
      if (GetProcessTimes(GetCurrentProcess(), &ftCreate, &ftExit, &ftKernel, &ftUser))
      {
         ULARGE_INTEGER ul;
         ul.LowPart = ftUser.dwLowDateTime;
         ul.HighPart = ftUser.dwHighDateTime;
         OutStream << "Time spent executing user code               \t: " << strDispTime(static_cast<double>(ul.QuadPart) * 1e-7, false) << endl;
         ul.LowPart = ftKernel.dwLowDateTime;
         ul.HighPart = ftKernel.dwHighDateTime;
         OutStream << "Time spent executing kernel code             \t: " << strDispTime(static_cast<double>(ul.QuadPart) * 1e-7, false) << endl;
      }
   }
 
   else
      OutStream << "Process timings                              \t: " << NA << endl;
 
   // Finally get more process statistics: this needs psapi.dll, so only proceed if it is present on this system
   HINSTANCE hDLL = LoadLibrary("psapi.dll");
 
   if (hDLL != NULL)
   {
      // The dll has been found
      typedef BOOL(__stdcall * DLLPROC)(HANDLE, PPROCESS_MEMORY_COUNTERS, DWORD);
      DLLPROC ProcAdd;
 
      // Try to get the address of the function we will call
      ProcAdd = (DLLPROC)GetProcAddress(hDLL, "GetProcessMemoryInfo");
 
      if (ProcAdd)
      {
         // Address was found
         PROCESS_MEMORY_COUNTERS pmc;
 
         // Now call the function
         if ((ProcAdd)(GetCurrentProcess(), &pmc, sizeof(pmc)))
         {
            OutStream << "Peak working set size                        \t: " << pmc.PeakWorkingSetSize / (1024.0 * 1024.0) << " Mb" << endl;
            OutStream << "Current working set size                     \t: " << pmc.WorkingSetSize / (1024.0 * 1024.0) << " Mb" << endl;
            OutStream << "Peak paged pool usage                        \t: " << pmc.QuotaPeakPagedPoolUsage / (1024.0 * 1024.0) << " Mb" << endl;
            OutStream << "Current paged pool usage                     \t: " << pmc.QuotaPagedPoolUsage / (1024.0 * 1024.0) << " Mb" << endl;
            OutStream << "Peak non-paged pool usage                    \t: " << pmc.QuotaPeakNonPagedPoolUsage / (1024.0 * 1024.0) << " Mb" << endl;
            OutStream << "Current non-paged pool usage                 \t: " << pmc.QuotaNonPagedPoolUsage / (1024.0 * 1024.0) << " Mb" << endl;
            OutStream << "Peak pagefile usage                          \t: " << pmc.PeakPagefileUsage / (1024.0 * 1024.0) << " Mb" << endl;
            OutStream << "Current pagefile usage                       \t: " << pmc.PagefileUsage / (1024.0 * 1024.0) << " Mb" << endl;
            OutStream << "No. of page faults                           \t: " << pmc.PageFaultCount << endl;
         }
      }
 
      // Free the memory used by the dll
      FreeLibrary(hDLL);
   }
 
#elif defined __GNUG__
   rusage ru;
 
   if (getrusage(RUSAGE_SELF, &ru) >= 0)
   {
      OutStream << "Time spent executing user code               \t: " << strDispTime(static_cast<double>(ru.ru_utime.tv_sec), false, true) << endl;
      // OutStream << "ru_utime.tv_usec                             \t: " << ru.ru_utime.tv_usec << endl;
      OutStream << "Time spent executing kernel code             \t: " << strDispTime(static_cast<double>(ru.ru_stime.tv_sec), false, true) << endl;
      // OutStream << "ru_stime.tv_usec                             \t: " << ru.ru_stime.tv_usec << endl;
      // OutStream << "Maximum resident set size                    \t: " << ru.ru_maxrss/1024.0 << " Mb" << endl;
      // OutStream << "ixrss (???)                                  \t: " << ru.ru_ixrss << endl;
      // OutStream << "Sum of rm_asrss (???)                        \t: " << ru.ru_idrss << endl;
      // OutStream << "isrss (???)                                  \t: " << ru.ru_isrss << endl;
      OutStream << "No. of page faults not requiring physical I/O\t: " << ru.ru_minflt << endl;
      OutStream << "No. of page faults requiring physical I/O    \t: " << ru.ru_majflt << endl;
      // OutStream << "No. of times swapped out of main memory      \t: " << ru.ru_nswap << endl;
      // OutStream << "No. of times performed input (read request)  \t: " << ru.ru_inblock << endl;
      // OutStream << "No. of times performed output (write request)\t: " << ru.ru_oublock << endl;
      // OutStream << "No. of signals received                      \t: " << ru.ru_nsignals << endl;
      OutStream << "No. of voluntary context switches            \t: " << ru.ru_nvcsw << endl;
      OutStream << "No. of involuntary context switches          \t: " << ru.ru_nivcsw << endl;
   }
 
   else
      OutStream << NA << endl;
 
#else
   OutStream << NA << endl;
#endif
 
   OutStream << endl;
 
#ifdef _OPENMP
#pragma omp parallel
   {
      if (0 == omp_get_thread_num())
      {
         OutStream << "Number of OpenMP threads                     \t: " << omp_get_num_threads() << endl;
         OutStream << "Number of OpenMP processors                  \t: " << omp_get_num_procs() << endl;
 
         LogStream << "Number of OpenMP threads                     \t: " << omp_get_num_threads() << endl;
         LogStream << "Number of OpenMP processors                  \t: " << omp_get_num_procs() << endl;
      }
   }
#endif
 
   time_t const tRunTime = m_tSysEndTime - m_tSysStartTime;
   struct tm* ptmRunTime = gmtime(&tRunTime);
 
   OutStream << "Time required for simulation                 \t: " << put_time(ptmRunTime, "%T") << endl;
   LogStream << "Time required for simulation                 \t: " << put_time(ptmRunTime, "%T") << endl;
 
   double const dSpeedUp = m_dSimDuration * 3600 / static_cast<double>(tRunTime);
   OutStream << setprecision(0);
   OutStream << "Time simulated / time required for simulation\t: " << dSpeedUp << " x faster than reality" << endl;
 
   LogStream << setprecision(0);
   LogStream << "Time simulated / time required for simulation\t: " << dSpeedUp << " x faster than reality" << endl;
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strGetErrorText(int const nErr)
{
   string strErr;
 
   switch (nErr)
   {
   case RTN_USER_ABORT:
      strErr = "run ended by user";
      break;
 
   case RTN_ERR_BADPARAM:
      strErr = "error in command-line parameter";
      break;
 
   case RTN_ERR_INI:
      strErr = "error reading initialization file";
      break;
 
   case RTN_ERR_CMEDIR:
      strErr = "error in directory name";
      break;
 
   case RTN_ERR_RUNDATA:
      strErr = "error reading run details file";
      break;
 
   case RTN_ERR_SCAPE_SHAPE_FUNCTION_FILE:
      strErr = "error reading SCAPE shape function file";
      break;
 
   case RTN_ERR_TIDEDATAFILE:
      strErr = "error reading tide data file";
      break;
 
   case RTN_ERR_LOGFILE:
      strErr = "error creating log file";
      break;
 
   case RTN_ERR_OUTFILE:
      strErr = "error creating text output file";
      break;
 
   case RTN_ERR_TSFILE:
      strErr = "error creating time series file";
      break;
 
   case RTN_ERR_DEMFILE:
      strErr = "error reading initial DEM file";
      break;
 
   case RTN_ERR_RASTER_FILE_READ:
      strErr = "error reading raster GIS file";
      break;
 
   case RTN_ERR_VECTOR_FILE_READ:
      strErr = "error reading vector GIS file";
      break;
 
   case RTN_ERR_MEMALLOC:
      strErr = "error allocating memory";
      break;
 
   case RTN_ERR_RASTER_GIS_OUT_FORMAT:
      strErr = "problem with raster GIS output format";
      break;
 
   case RTN_ERR_VECTOR_GIS_OUT_FORMAT:
      strErr = "problem with vector GIS output format";
      break;
 
   case RTN_ERR_TEXT_FILE_WRITE:
      strErr = "error writing text output file";
      break;
 
   case RTN_ERR_RASTER_FILE_WRITE:
      strErr = "error writing raster GIS output file";
      break;
 
   case RTN_ERR_VECTOR_FILE_WRITE:
      strErr = "error writing vector GIS output file";
      break;
 
   case RTN_ERR_TIMESERIES_FILE_WRITE:
      strErr = "error writing time series output file";
      break;
 
   case RTN_ERR_LINETOGRID:
      strErr = "error putting linear feature onto raster grid";
      break;
 
   case RTN_ERR_NOSEACELLS:
      strErr = "no sea cells found";
      break;
 
   case RTN_ERR_GRIDTOLINE:
      strErr = "error when searching grid for linear feature";
      break;
 
   case RTN_ERR_NOCOAST:
      strErr = "no coastlines found. Is the SWL correct?";
      break;
 
   case RTN_ERR_PROFILEWRITE:
      strErr = "error writing coastline-normal profiles";
      break;
 
   case RTN_ERR_TIMEUNITS:
      strErr = "error in time units";
      break;
 
   case RTN_ERR_NO_SOLUTION_FOR_ENDPOINT:
      strErr = "no solution when finding end point for coastline-normal line";
      break;
 
   // case RTN_ERR_PROFILE_ENDPOINT_AT_GRID_EDGE:
   // strErr = "end point for coastline-normal line is at the grid edge";
   // break;
   case RTN_ERR_PROFILE_ENDPOINT_IS_INLAND:
      strErr = "end point for coastline-normal line is not in the contiguous sea";
      break;
 
   case RTN_ERR_CLIFFNOTCH:
      strErr = "cliff notch is above sediment top elevation";
      break;
 
   case RTN_ERR_CLIFFDEPOSIT:
      strErr = "unable to deposit sediment from cliff collapse";
      break;
 
   case RTN_ERR_PROFILESPACING:
      strErr = "coastline-normal profiles are too closely spaced";
      break;
 
   case RTN_ERR_NO_PROFILES_1:
      strErr = "no coastline-normal profiles created, check the SWL";
      break;
 
   case RTN_ERR_NO_PROFILES_2:
      strErr = "no coastline-normal profiles created during rasterization";
      break;
 
   case RTN_ERR_EDGE_OF_GRID:
      strErr = "hit grid edge when eroding beach";
      break;
 
   case RTN_ERR_NO_SEAWARD_END_OF_PROFILE_BEACH_EROSION:
      strErr = "could not locate seaward end of profile when creating Dean profile for beach erosion";
      break;
 
   case RTN_ERR_NO_SEAWARD_END_OF_PROFILE_UPCOAST_BEACH_DEPOSITION:
      strErr = "could not locate seaward end of profile when creating Dean profile for beach deposition (up-coast)";
      break;
 
   case RTN_ERR_NO_SEAWARD_END_OF_PROFILE_DOWNCOAST_BEACH_DEPOSITION:
      strErr = "could not locate seaward end of profile when creating Dean profile for beach deposition (down-coast)";
      break;
 
   case RTN_ERR_LANDFORM_TO_GRID:
      strErr = "updating grid with landforms";
      break;
 
   case RTN_ERR_NO_TOP_LAYER:
      strErr = "no top layer of sediment";
      break;
 
   case RTN_ERR_NO_ADJACENT_POLYGON:
      strErr = "problem with polygon-to-polygon sediment routing sequence";
      break;
 
   case RTN_ERR_BAD_MULTILINE:
      strErr = "inconsistent multiline";
      break;
 
   case RTN_ERR_CANNOT_INSERT_POINT:
      strErr = "cannot insert point into multiline";
      break;
 
   case RTN_ERR_CANNOT_ASSIGN_COASTAL_LANDFORM:
      strErr = "cannot assign coastal landform";
      break;
 
   case RTN_ERR_SHADOW_ZONE_FLOOD_FILL_NOGRID:
      strErr = "start point for cell-by-cell fill of wave shadow zone is outside grid";
      break;
 
   case RTN_ERR_SHADOW_ZONE_FLOOD_START_POINT:
      strErr = "could not find start point for cell-by-cell fill of wave shadow zone";
      break;
 
   case RTN_ERR_CSHORE_EMPTY_PROFILE:
      strErr = "empty profile during during CShore wave propagation";
      break;
 
   case RTN_ERR_CSHORE_FILE_INPUT:
      strErr = "creating file for CShore input";
      break;
 
   case RTN_ERR_READING_CSHORE_FILE_OUTPUT:
      strErr = "reading CShore output file";
      break;
 
   case RTN_ERR_WAVE_INTERPOLATION_LOOKUP:
      strErr = "during wave interpolation lookup";
      break;
 
   case RTN_ERR_GRIDCREATE:
      strErr = "while running GDALGridCreate()";
      break;
 
   case RTN_ERR_COAST_CANT_FIND_EDGE_CELL:
      strErr = "cannot find edge cell while constructing grid-edge profile";
      break;
 
   case RTN_ERR_CSHORE_ERROR:
      strErr = "CShore did not finish correctly";
      break;
 
   case RTN_ERR_NO_CELL_UNDER_COASTLINE:
      strErr = "Could not find cell under coastline";
      break;
 
   case RTN_ERR_OPEN_DEEP_WATER_WAVE_DATA:
      strErr = "opening deep sea wave time series file";
      break;
 
   case RTN_ERR_READING_DEEP_WATER_WAVE_DATA:
      strErr = "reading deep sea wave time series file";
      break;
 
   case RTN_ERR_BOUNDING_BOX:
      strErr = "finding edges of the bounding box";
      break;
 
   case RTN_ERR_READING_SEDIMENT_INPUT_EVENT:
      strErr = "reading sediment input event time series file";
      break;
 
   case RTN_ERR_SEDIMENT_INPUT_EVENT:
      strErr = "simulating sediment input event";
      break;
 
   case RTN_ERR_SEDIMENT_INPUT_EVENT_LOCATION:
      strErr = "location of sediment input event is outside grod";
      break;
 
   case RTN_ERR_WAVESTATION_LOCATION:
      strErr = "location of wavestation is outside grid";
      break;
 
   case RTN_ERR_CLIFF_NOT_IN_POLYGON:
      strErr = "cliff not in polygon";
      break;
 
   case RTN_ERR_CELL_MARKED_PROFILE_COAST_BUT_NOT_PROFILE:
      strErr = "Cell marked as profile coast but not as profile";
      break;
 
   case RTN_ERR_TRACING_FLOOD:
      strErr = "error tracing flood line on grid";
      break;
 
   case RTN_ERR_NO_START_FINISH_POINTS_TRACING_COAST:
      strErr = "error tracing coastline on grid, no coast start-finish points found";
      break;
 
   case RTN_ERR_NO_VALID_COAST:
      strErr = "error tracing coastline on grid, no valid coast found";
      break;
 
   case RTN_ERR_REPEATING_WHEN_TRACING_COAST:
      strErr = "error tracing coastline on grid, coast search just repeats";
      break;
 
   case RTN_ERR_ZERO_LENGTH_COAST:
      strErr = "error tracing coastline on grid, zero-length coast found";
      break;
 
   case RTN_ERR_COAST_TOO_SMALL:
      strErr = "error tracing coastline on grid, coast below minimum permitted length";
      break;
 
   case RTN_ERR_IGNORING_COAST:
      strErr = "error tracing coastline on grid, coast ignored";
      break;
 
   case RTN_ERR_TOO_LONG_TRACING_COAST:
      strErr = "error tracing coastline on grid, too many times round tracing loop";
      break;
 
   case RTN_ERR_CELL_NOT_FOUND_IN_HIT_PROFILE_DIFFERENT_COASTS:
      strErr = "intersection cell not found in hit profile";
      break;
 
   case RTN_ERR_UNKNOWN:
      strErr = "unknown error";
      break;
 
   default:
      // should never get here
      strErr = "totally unknown error";
   }
 
   return strErr;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::DoSimulationEnd(int const nRtn)
{
   // If we don't know the time that the run ended (e.g. because it did not finish correctly), then get it now
   if (m_tSysEndTime == 0)
      m_tSysEndTime = time(nullptr);
 
   switch (nRtn)
   {
   case (RTN_OK):
      // normal ending
      cout << RUN_END_NOTICE << put_time(localtime(&m_tSysEndTime), "%T %A %d %B %Y") << endl;
      break;
 
   case (RTN_HELP_ONLY):
   case (RTN_CHECK_ONLY):
      return;
 
   default:
      // Aborting because of some error
      cerr << RUN_END_NOTICE << "iteration " << m_ulIter << ERROR_NOTICE << nRtn << " (" << strGetErrorText(nRtn) << ") on " << put_time(localtime(&m_tSysEndTime), "%T %A %d %B %Y") << endl;
 
      if (m_ulIter > 1)
      {
         // If the run has actually started, then output all GIS files: this is very helpful in tracking down problems
         m_bSaveGISThisIter = true;
         m_nGISSave = 998; // Will get incremented to 999 when we write the files
         bSaveAllRasterGISFiles();
         bSaveAllVectorGISFiles();
      }
 
      // Write the error message to the logfile and to stdout
      if (LogStream && LogStream.is_open())
      {
         LogStream << ERR << strGetErrorText(nRtn) << " (error code " << nRtn << ") on " << put_time(localtime(&m_tSysEndTime), "%T %A %d %B %Y") << endl;
         LogStream.flush();
      }
 
      if (OutStream && OutStream.is_open())
      {
         OutStream << ERR << strGetErrorText(nRtn) << " (error code " << nRtn << ") on " << put_time(localtime(&m_tSysEndTime), "%T %A %d %B %Y") << endl;
         OutStream.flush();
      }
   }
 
#ifdef __GNUG__
 
   if (isatty(fileno(stdout)))
   {
      // Stdout is connected to a tty, so not running as a background job
      // cout << endl
      // << PRESS_KEY;
      // cout.flush();
      // getchar();
   }
   else
   {
      // Stdout is not connected to a tty, so must be running in the background; if we have something entered for the email address, then send an email
      if (! m_strMailAddress.empty())
      {
         cout << SEND_EMAIL << m_strMailAddress << endl;
 
         string strCmd("echo \"");
 
         stringstream ststrTmp;
         ststrTmp << put_time(localtime(&m_tSysEndTime), "%T on %A %d %B %Y") << endl;
 
         // Send an email using Linux/Unix mail command
         if (RTN_OK == nRtn)
         {
            // Finished normally
            strCmd.append("Simulation ");
            strCmd.append(m_strRunName);
            strCmd.append(", running on ");
            strCmd.append(strGetComputerName());
            strCmd.append(", completed normally at ");
            strCmd.append(ststrTmp.str());
            strCmd.append("\" | mail -s \"");
            strCmd.append(PROGRAM_NAME);
            strCmd.append(": normal completion\" ");
            strCmd.append(m_strMailAddress);
         }
 
         else
         {
            // Error, so give some information to help debugging
            strCmd.append("Simulation ");
            strCmd.append(m_strRunName);
            strCmd.append(", running on ");
            strCmd.append(strGetComputerName());
            strCmd.append(", aborted with error code ");
            strCmd.append(to_string(nRtn));
            strCmd.append(": ");
            strCmd.append(strGetErrorText(nRtn));
            strCmd.append(" at timestep ");
            strCmd.append(to_string(m_ulIter));
            strCmd.append(" (");
            strCmd.append(strDispSimTime(m_dSimElapsed));
            strCmd.append(").\n\nThis message sent at ");
            strCmd.append(ststrTmp.str());
            strCmd.append("\" | mail -s \"");
            strCmd.append(PROGRAM_NAME);
            strCmd.append(": ERROR\" ");
            strCmd.append(m_strMailAddress);
         }
 
         int const nRet = system(strCmd.c_str());
 
         if (WEXITSTATUS(nRet) != 0)
            cerr << ERR << EMAIL_ERROR << endl;
      }
   }
 
#endif
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::pstrChangeToBackslash(string const* strIn)
{
   string strOut(*strIn);
   strOut.replace(strOut.begin(), strOut.end(), '/', '\\');
   return strOut;
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::pstrChangeToForwardSlash(string const* strIn)
{
   string strOut(*strIn);
   strOut.replace(strOut.begin(), strOut.end(), '\\', '/');
   return strOut;
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strTrimLeft(string const* strIn)
{
   // Trim leading spaces
   size_t const nStartpos = strIn->find_first_not_of(" \t");
 
   if (nStartpos == string::npos)
      return *strIn;
 
   else
      return strIn->substr(nStartpos);
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strTrimRight(string const* strIn)
{
   string strTmp(*strIn);
 
   // Remove any stray carriage returns (can happen if file was edited in Windows)
   strTmp.erase(remove(strTmp.begin(), strTmp.end(), '\r'), strTmp.end());
 
   // Trim trailing spaces
   size_t const nEndpos = strTmp.find_last_not_of(" \t");
 
   if (nEndpos == string::npos)
      return strTmp;
 
   else
      return strTmp.substr(0, nEndpos + 1);
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strTrim(string const* strIn)
{
   string strTmp = *strIn;
 
   // Remove any stray carriage returns (can happen if file was edited in Windows)
   strTmp.erase(remove(strTmp.begin(), strTmp.end(), '\r'), strTmp.end());
 
   // Trim trailing spaces
   size_t nPos = strTmp.find_last_not_of(" \t");
 
   if (nPos != string::npos)
      strTmp.resize(nPos + 1);
 
   // Trim leading spaces
   nPos = strTmp.find_first_not_of(" \t");
 
   if (nPos != string::npos)
      strTmp = strTmp.substr(nPos);
 
   return strTmp;
}
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strToLower(string const* strIn)
{
   string strOut = *strIn;
   transform(strIn->begin(), strIn->end(), strOut.begin(), tolower);
   return strOut;
}
 
//===============================================================================================================================
// Returns the upper case version of an string, leaving the original unchanged
//===============================================================================================================================
// string CSimulation::strToUpper(string const* strIn)
// {
// string strOut = *strIn;
// transform(strIn->begin(), strIn->end(), strOut.begin(), toupper);
// return strOut;
// }
 
//===============================================================================================================================
//===============================================================================================================================
string CSimulation::strRemoveSubstr(string* pStrIn, string const* pStrSub)
{
   size_t const nPos = pStrIn->find(*pStrSub);
 
   if (nPos != string::npos)
   {
      // OK, found the substring
      pStrIn->replace(nPos, pStrSub->size(), "");
      return strTrim(pStrIn);
   }
 
   else
   {
      // If not found, return the string unchanged
      return *pStrIn;
   }
}
 
//===============================================================================================================================
//===============================================================================================================================
vector<string>* CSimulation::VstrSplit(string const* s, char const delim, vector<string>* elems)
{
   stringstream ss(*s);
   string item;
 
   while (getline(ss, item, delim))
   {
      if (!item.empty())
         elems->push_back(item);
   }
 
   return elems;
}
 
//===============================================================================================================================
//===============================================================================================================================
vector<string> CSimulation::VstrSplit(string const* s, char const delim)
{
   vector<string> elems;
   VstrSplit(s, delim, &elems);
   return elems;
}
 
// //===============================================================================================================================
// //! Calculates the vector cross product of three points
// //===============================================================================================================================
// double CSimulation::dCrossProduct(double const dX1, double const dY1, double const dX2, double const dY2, double const dX3, double const dY3)
// {
//    // Based on code at http://debian.fmi.uni-sofia.bg/~sergei/cgsr/docs/clockwise.htm
// return (dX2 - dX1) * (dY3 - dY2) - ((dY2 - dY1) * (dX3 - dX2));
// }
 
// //===============================================================================================================================
// //! Calculates the mean of a pointer to a vector of doubles
// //===============================================================================================================================
// double CSimulation::dGetMean(vector<double> const* pV)
// {
// double dSum = accumulate(pV->begin(), pV->end(), 0.0);
// double dMean = dSum / static_cast<double>(pV->size());
// return dMean;
// }
 
// //===============================================================================================================================
// //! Calculates the standard deviation of a pointer to a vector of doubles. From http://stackoverflow.com/questions/7616511/calculate-mean-and-standard-deviation-from-a-vector-of-samples-in-c-using-boos
// //===============================================================================================================================
// double CSimulation::dGetStdDev(vector<double> const* pV)
// {
// double dSum = accumulate(pV->begin(), pV->end(), 0.0);
// double dMean = dSum / static_cast<double>(pV->size());
//
// double dSqSum = inner_product(pV->begin(), pV->end(), pV->begin(), 0.0);
// double dStdDev = sqrt(dSqSum / static_cast<double>(pV->size()) - dMean * dMean);
//
// return dStdDev;
// }
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AppendEnsureNoGap(vector<CGeom2DIPoint>* pVPtiPoints, CGeom2DIPoint const* pPti)
{
   int const nX = pPti->nGetX();
   int const nY = pPti->nGetY();
   int const nXLast = pVPtiPoints->back().nGetX();
   int const nYLast = pVPtiPoints->back().nGetY();
   int const nXDiff = nX - nXLast;
   int const nYDiff = nY - nYLast;
   int const nXDiffA = tAbs(nXDiff);
   int const nYDiffA = tAbs(nYDiff);
   int const nDiff = tMax(nXDiffA, nYDiffA);
 
   if (nDiff > 1)
   {
      // We have a gap
      double
          dXInc = 0,
          dYInc = 0;
 
      if (nXDiffA > 1)
         dXInc = static_cast<double>(nXDiff) / nDiff;
 
      if (nYDiffA > 1)
         dYInc = static_cast<double>(nYDiff) / nDiff;
 
      for (int n = 1; n < nDiff; n++)
      {
         CGeom2DIPoint const Pti(nXLast + nRound(n * dXInc), nYLast + nRound(n * dYInc));
         pVPtiPoints->push_back(Pti);
      }
   }
 
   pVPtiPoints->push_back(CGeom2DIPoint(nX, nY));
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::CalcDeanProfile(vector<double>* pdVDeanProfile, double const dInc, double const dDeanTopElev, double const dA, bool const bDeposition, int const nSeawardOffset, double const dStartCellElev)
{
   double dDistFromProfileStart = 0;
 
   if (bDeposition)
   {
      // This Dean profile is for deposition i.e. seaward displacement of the profile
      pdVDeanProfile->at(0) = dStartCellElev; // Is talus-top elev for cliffs, coast elevation for coasts
 
      for (int n = 1; n < static_cast<int>(pdVDeanProfile->size()); n++)
      {
         if (n <= nSeawardOffset)
            // As we extend the profile seaward, the elevation of any points coastward of the new coast point of the Dean profile are set to the elevation of the original coast or the talus top (is this realistic for talus?)
            pdVDeanProfile->at(n) = dStartCellElev;
 
         else
         {
            double const dDistBelowTop = dA * pow(dDistFromProfileStart, DEAN_POWER);
            pdVDeanProfile->at(n) = dDeanTopElev - dDistBelowTop;
 
            dDistFromProfileStart += dInc;
         }
      }
   }
 
   else
   {
      // This Dean profile is for erosion i.e. landward displacement of the profile
      for (int n = 0; n < static_cast<int>(pdVDeanProfile->size()); n++)
      {
         double const dDistBelowTop = dA * pow(dDistFromProfileStart, DEAN_POWER);
         pdVDeanProfile->at(n) = dDeanTopElev - dDistBelowTop;
 
         dDistFromProfileStart += dInc;
      }
   }
}
 
//===============================================================================================================================
//===============================================================================================================================
double CSimulation::dSubtractProfiles(vector<double> const* pdVFirstProfile, vector<double> const* pdVSecondProfile, vector<bool> const* pbVIsValid)
{
   double dTotElevDiff = 0;
 
   // Note that this assumes that all three vectors are of equal length, should really check this
   for (int n = 0; n < static_cast<int>(pdVFirstProfile->size()); n++)
   {
      if (pbVIsValid->at(n))
      {
         double const dProfileDiff = pdVFirstProfile->at(n) - pdVSecondProfile->at(n);
 
         dTotElevDiff += dProfileDiff;
      }
   }
 
   //    // DEBUG CODE -----------------------------------------------------
   // LogStream << endl;
   // LogStream << "First profile = ";
   // for (int n = 0; n < static_cast<int>(pdVFirstProfile->size()); n++)
   // {
   // LogStream << pdVFirstProfile->at(n) << " ";
   // }
   // LogStream << endl;
   // LogStream << "Second profile = ";
   // for (int n = 0; n < static_cast<int>(pdVFirstProfile->size()); n++)
   // {
   // LogStream << pdVSecondProfile->at(n) << " ";
   // }
   // LogStream << endl;
   // LogStream << "Difference = ";
   // for (int n = 0; n < static_cast<int>(pdVFirstProfile->size()); n++)
   // {
   // LogStream << pdVFirstProfile->at(n) - pdVSecondProfile->at(n) << " ";
   // }
   // LogStream << endl;
   //    // DEBUG CODE -----------------------------------------------------
 
   return dTotElevDiff;
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::CalcDepthOfClosure(void)
{
   double
       dDeepWaterWaveHeight,
       dDeepWaterPeriod;
 
   if (m_bSingleDeepWaterWaveValues)
   {
      dDeepWaterWaveHeight = m_dAllCellsDeepWaterWaveHeight;
      dDeepWaterPeriod = m_dAllCellsDeepWaterWavePeriod;
   }
 
   else
   {
      dDeepWaterWaveHeight = m_dMaxUserInputWaveHeight;
      dDeepWaterPeriod = m_dMaxUserInputWavePeriod;
   }
 
   // TODO 051 Calculate depth of closure using 'average of the maximum values observed during a typical year'
   // dL = 2.28 * Hsx − (68.5 * Hsx^2 / (g * Tsx^2))
   // where:
   // Hsx is the nearshore storm wave height that is exceeded only 12 hours each year
   // Tsx is the associated wave period
   // from Hallermeier, R.J. (1978). Uses for a calculated limit depth to beach erosion. Proc. 16th Coastal Engineering Conf., ASCE, New York. Pp 1493 - 1512
   //
   // For the time being, and since we assume wave height and period constant just use the actual wave height and period to calculate the depth of closure
   // m_dDepthOfClosure = (2.28 * dDeepWaterWaveHeight) - (68.5 * dDeepWaterWaveHeight * dDeepWaterWaveHeight / (m_dG * dDeepWaterPeriod * dDeepWaterPeriod));
 
   // An alternative (which produces smaller depth of closure estimates) is Birkemeier (1985) TODO 007 Full reference needed
   // dL = 1.75 * Hsx - (57.9 * Hsx^2/ (g * Tsx^2))
   m_dDepthOfClosure = (1.75 * dDeepWaterWaveHeight) - (57.9 * dDeepWaterWaveHeight * dDeepWaterWaveHeight / (m_dG * dDeepWaterPeriod * dDeepWaterPeriod));
}
 
// //===============================================================================================================================
// //! Tests a reference to a string to see if it is numeric (modified from https://tfetimes.com/c-determine-if-a-string-is-numeric/)
// //===============================================================================================================================
// bool CSimulation::bIsNumeric(string const*strIn)
// {
// return all_of(strIn->begin(), strIn->end(), isdigit);
// }
 
//===============================================================================================================================
//===============================================================================================================================
bool CSimulation::bParseDate(string const* strDate, int& nDay, int& nMonth, int& nYear)
{
   vector<string> VstrTmp = VstrSplit(strDate, SLASH);
 
   if (VstrTmp.size() < 3)
   {
      cerr << "date string must include day, month, and year '" << strDate << "'" << endl;
      return false;
   }
 
   // Sort out day
   if (! bIsStringValidInt(VstrTmp[0]))
   {
      cerr << "invalid integer for day in date '" << strDate << "'" << endl;
      return false;
   }
 
   nDay = stoi(VstrTmp[0]);
 
   if ((nDay < 1) || (nDay > 31))
   {
      cerr << "day must be between 1 and 31 in date '" << strDate << "'" << endl;
      return false;
   }
 
   // Sort out month
   if (! bIsStringValidInt(VstrTmp[1]))
   {
      cerr << "invalid integer for month in date '" << strDate << "'" << endl;
      return false;
   }
 
   nMonth = stoi(VstrTmp[1]);
 
   if ((nMonth < 1) || (nMonth > 12))
   {
      cerr << "month must be between 1 and 12 in date '" << strDate << "'" << endl;
      return false;
   }
 
   // Sort out year
   if (! bIsStringValidInt(VstrTmp[2]))
   {
      cerr << "invalid integer for year in date '" << strDate << "'" << endl;
      return false;
   }
 
   nYear = stoi(VstrTmp[2]);
 
   if (nYear < 0)
   {
      cerr << "year must be > 0 in date '" << strDate << "'" << endl;
      return false;
   }
 
   return true;
}
 
//===============================================================================================================================
//===============================================================================================================================
bool CSimulation::bParseTime(string const* strTime, int& nHour, int& nMin, int& nSec)
{
   vector<string> VstrTmp = VstrSplit(strTime, DASH);
 
   if (VstrTmp.size() < 3)
   {
      cerr << "time string must include hours, minutes, and seconds '" << strTime << "'" << endl;
      return false;
   }
 
   // Sort out hour
   if (! bIsStringValidInt(VstrTmp[0]))
   {
      cerr << "invalid integer for hours in time '" << strTime << "'" << endl;
      return false;
   }
 
   nHour = stoi(VstrTmp[0]);
 
   if ((nHour < 0) || (nHour > 23))
   {
      cerr << "hour must be between 0 and 23 in time '" << strTime << "'" << endl;
      return false;
   }
 
   // Sort out minutes
   if (! bIsStringValidInt(VstrTmp[1]))
   {
      cerr << "invalid integer for minutes in time '" << strTime << "'" << endl;
      return false;
   }
 
   nMin = stoi(VstrTmp[1]);
 
   if ((nMin < 0) || (nMin > 59))
   {
      cerr << "minutes must be betwen 0 and 59 in time '" << strTime << "'" << endl;
      return false;
   }
 
   // Sort out seconds
   if (! bIsStringValidInt(VstrTmp[2]))
   {
      cerr << "invalid integer for seconds in time '" << strTime << "'" << endl;
      return false;
   }
 
   nSec = stoi(VstrTmp[2]);
 
   if ((nSec < 0) || (nSec > 59))
   {
      cerr << "seconds must be between 0 and 59 in time '" << strTime << "'" << endl;
      return false;
   }
 
   return true;
}
 
//===============================================================================================================================
//===============================================================================================================================
unsigned long CSimulation::ulConvertToTimestep(string const* pstrIn) const
{
   unsigned long ulTimeStep = 0;
 
   // Convert to lower case, remove leading and trailing whitespace
   string strDate = strToLower(pstrIn);
   strDate = strTrim(&strDate);
 
   if (strDate.find("hour") != string::npos)
   {
      // OK, this is a number of hours (a relative time, from the start of simulation)
      vector<string> VstrTmp = VstrSplit(&strDate, SPACE);
 
      if ((VstrTmp.size() < 2) || (! bIsStringValidInt(VstrTmp[0])))
      {
         cerr << "Error in number of hours '" + strDate + "' for sediment input event" << endl;
         return SEDIMENT_INPUT_EVENT_ERROR;
      }
 
      double const dHours = stod(strTrim(&VstrTmp[0]));
 
      if (dHours > m_dSimDuration)
      {
         cerr << "Sediment input event '" + strDate + "' occurs after end of simulation" << endl;
         return SEDIMENT_INPUT_EVENT_ERROR;
      }
 
      ulTimeStep = static_cast<unsigned long>(dRound(dHours / m_dTimeStep));
   }
 
   else if (strDate.find("day") != string::npos)
   {
      // OK, this is a number of days (a relative time, from the start of simulation)
      vector<string> VstrTmp = VstrSplit(&strDate, SPACE);
 
      if ((VstrTmp.size() < 2) || (! bIsStringValidInt(VstrTmp[0])))
      {
         cerr << "Error in number of days '" + strDate + "' for sediment input event" << endl;
         return SEDIMENT_INPUT_EVENT_ERROR;
      }
 
      double const dHours = stod(strTrim(&VstrTmp[0])) * 24;
 
      if (dHours > m_dSimDuration)
      {
         cerr << "Sediment input event '" + strDate + "' occurs after end of simulation" << endl;
         return SEDIMENT_INPUT_EVENT_ERROR;
      }
 
      ulTimeStep = static_cast<unsigned long>(dRound(dHours / m_dTimeStep));
   }
 
   else
   {
      // This is an absolute time/date in the format hh-mm-ss dd/mm/yyyy
      vector<string> VstrTmp = VstrSplit(&strDate, SPACE);
 
      if (VstrTmp.size() < 2)
      {
         cerr << "Error in time/date '" + strDate + "' of sediment input event" << endl;
         return SEDIMENT_INPUT_EVENT_ERROR;
      }
 
      int nHour = 0;
      int nMin = 0;
      int nSec = 0;
 
      // OK, first sort out the time
      if (! bParseTime(&VstrTmp[0], nHour, nMin, nSec))
      {
         cerr << "Error in time '" + VstrTmp[0] + "' of sediment input event" << endl;
         return SEDIMENT_INPUT_EVENT_ERROR;
      }
 
      int nDay = 0;
      int nMonth = 0;
      int nYear = 0;
 
      // Now sort out the time
      if (! bParseDate(&VstrTmp[1], nDay, nMonth, nYear))
      {
         cerr << "Error in date '" + VstrTmp[1] + "' of sediment input event" << endl;
         return SEDIMENT_INPUT_EVENT_ERROR;
      }
 
      // This is modified from https://stackoverflow.com/questions/14218894/number-of-days-between-two-dates-c
      struct tm tmSimStart = {};
      tmSimStart.tm_sec = m_nSimStartSec;
      tmSimStart.tm_min = m_nSimStartMin;
      tmSimStart.tm_hour = m_nSimStartHour;
      tmSimStart.tm_mday = m_nSimStartDay;
      tmSimStart.tm_mon = m_nSimStartMonth - 1;
      tmSimStart.tm_year = m_nSimStartYear - 1900;
 
      struct tm tmSimEvent = {};
      tmSimEvent.tm_sec = nSec;
      tmSimEvent.tm_min = nMin;
      tmSimEvent.tm_hour = nHour;
      tmSimEvent.tm_mday = nDay;
      tmSimEvent.tm_mon = nMonth - 1;
      tmSimEvent.tm_year = nYear - 1900;
 
      time_t const tStart = mktime(&tmSimStart);
      time_t const tEvent = mktime(&tmSimEvent);
 
      if (tStart == (time_t)(-1))
      {
         cerr << "Error in simulation start time/date" << endl;
         return SEDIMENT_INPUT_EVENT_ERROR;
      }
 
      if (tEvent == (time_t)(-1))
      {
         cerr << "Error in time/date '" + strDate + "' of sediment input event" << endl;
         return SEDIMENT_INPUT_EVENT_ERROR;
      }
 
      double const dHours = difftime(tEvent, tStart) / (60 * 60);
 
      if (dHours < 0)
      {
         cerr << "Sediment input event '" + strDate + "' occurs before start of simulation" << endl;
         return SEDIMENT_INPUT_EVENT_ERROR;
      }
 
      if (dHours > m_dSimDuration)
      {
         cerr << "Sediment input event '" + strDate + "' occurs after end of simulation" << endl;
         return SEDIMENT_INPUT_EVENT_ERROR;
      }
 
      ulTimeStep = static_cast<unsigned long>(dHours / m_dTimeStep);
   }
 
   return ulTimeStep;
}
 
//===============================================================================================================================
//===============================================================================================================================
bool CSimulation::bIsInterventionCell(int const nX, int const nY) const
{
   if (m_pRasterGrid->m_Cell[nX][nY].pGetLandform()->nGetLFCategory() == LF_CAT_INTERVENTION)
      return true;
 
   return false;
}
 
//===============================================================================================================================
//===============================================================================================================================
int CSimulation::nGetCoastPolygonSize(void) const
{
   return static_cast<int>(m_pVCoastPolygon.size());
}
 
//===============================================================================================================================
//===============================================================================================================================
CGeomCoastPolygon* CSimulation::pGetPolygon(int const nPoly) const
{
   // TODO 055 No check to see if nPoly < m_pVCoastPolygon.size()
   return m_pVCoastPolygon[nPoly];
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::AppendPolygon(CGeomCoastPolygon* pPolygon)
{
   m_pVCoastPolygon.push_back(pPolygon);
}
 
//===============================================================================================================================
//===============================================================================================================================
void CSimulation::DoEndOfRunDeletes(void)
{
   // Clear all vector coastlines, profiles, and polygons
   for (int i = 0; i < static_cast<int>(m_pVCoastPolygon.size()); i++)
      delete m_pVCoastPolygon[i];
 
   m_pVCoastPolygon.clear();
   m_VCoast.clear();
 
   // m_VFloodWaveSetup.clear();
   m_VFloodWaveSetupSurge.clear();
   m_VFloodWaveSetupSurgeRunup.clear();
}