/******************************************************************************
 *
 * Project:  PROJ
 * Purpose:  Test deformation model
 * Author:   Even Rouault <even dot rouault at spatialys dot com>
 *
 ******************************************************************************
 * Copyright (c) 2020, Even Rouault <even dot rouault at spatialys dot com>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/

#include "gtest_include.h"

// to be able to use internal::toString
#define FROM_PROJ_CPP
#include "proj/internal/internal.hpp"

#include "proj.h"
#include "proj_internal.h"

// Silence C4702 (unreachable code) due to some dummy implementation of the
// interfaces of defmodel.hpp
#ifdef _MSC_VER
#pragma warning(push)
#pragma warning(disable : 4702)
#endif

#define PROJ_COMPILATION
#define DEFORMATON_MODEL_NAMESPACE TestDeformationModel
#include "transformations/defmodel.hpp"

using namespace DEFORMATON_MODEL_NAMESPACE;

namespace {

constexpr double modelMinX = 158;
constexpr double modelMinY = -58;
constexpr double modelMaxX = 194;
constexpr double modelMaxY = -25;

static json getMinValidContent() {
    json j;
    j["file_type"] = "GeoTIFF";
    j["format_version"] = "1.0";
    j["source_crs"] = "EPSG:4959";
    j["target_crs"] = "EPSG:7907";
    j["definition_crs"] = "EPSG:4959";
    j["extent"]["type"] = "bbox";
    j["extent"]["parameters"] = {
        {"bbox", {modelMinX, modelMinY, modelMaxX, modelMaxY}}};
    j["time_extent"]["first"] = "1900-01-01T00:00:00Z";
    j["time_extent"]["last"] = "2050-01-01T00:00:00Z";
    j["components"] = json::array();

    return j;
}

// ---------------------------------------------------------------------------

constexpr int IDX_CONSTANT = 0;
constexpr int IDX_VELOCITY = 1;
constexpr int IDX_STEP = 2;
constexpr int IDX_REVERSE_STEP = 3;
constexpr int IDX_PIECEWISE = 4;
constexpr int IDX_EXPONENTIAL = 5;

static json getFullValidContent() {
    json j(getMinValidContent());
    j["name"] = "name";
    j["version"] = "version";
    j["publication_date"] = "2018-07-01T00:00:00Z";
    j["license"] = "license";
    j["description"] = "description";
    j["authority"]["name"] = "authority_name";
    j["authority"]["url"] = "authority_url";
    j["authority"]["address"] = "authority_address";
    j["authority"]["email"] = "authority_email";
    j["links"] = {{{"href", "href"},
                   {"rel", "rel"},
                   {"type", "type"},
                   {"title", "title"}}};
    j["reference_epoch"] = "2000-01-01T00:00:00Z";
    j["uncertainty_reference_epoch"] = "2018-12-15T00:00:00Z";
    j["horizontal_offset_method"] = "addition";
    j["horizontal_offset_unit"] = "metre";
    j["vertical_offset_unit"] = "metre";
    j["horizontal_uncertainty_type"] = "circular 95% confidence limit";
    j["horizontal_uncertainty_unit"] = "metre";
    j["vertical_uncertainty_type"] = "95% confidence limit";
    j["vertical_uncertainty_unit"] = "metre";
    j["components"] = {{
        {"description", "description"},
        {"displacement_type", "horizontal"},
        {"uncertainty_type", "none"},
        {"horizontal_uncertainty", 0.01},
        {"vertical_uncertainty", 0.02},
        {"extent",
         {{"type", "bbox"},
          {"parameters",
           {{"bbox", {modelMinX, modelMinY, modelMaxX, modelMaxY}}}}}},
        {"spatial_model",
         {
             {"type", "GeoTIFF"},
             {"interpolation_method", "bilinear"},
             {"filename", "nzgd2000-ndm-grid02.tif"},
             {"md5_checksum", "49fce8ab267be2c8d00d43683060a032"},
         }},
        {"time_function",
         {
             {"type", "constant"},
             {"parameters", json::object()},
         }},
    }};

    j["components"].push_back(j["components"][0]);
    j["components"][IDX_VELOCITY]["time_function"] = {
        {"type", "velocity"},
        {"parameters", {{"reference_epoch", "2000-01-01T00:00:00Z"}}},
    };

    j["components"].push_back(j["components"][0]);
    j["components"][IDX_STEP]["time_function"] = {
        {"type", "step"},
        {"parameters", {{"step_epoch", "2000-01-01T00:00:00Z"}}},
    };

    j["components"].push_back(j["components"][0]);
    j["components"][IDX_REVERSE_STEP]["time_function"] = {
        {"type", "reverse_step"},
        {"parameters", {{"step_epoch", "2000-01-01T00:00:00Z"}}},
    };

    j["components"].push_back(j["components"][0]);
    j["components"][IDX_PIECEWISE]["time_function"] = {
        {"type", "piecewise"},
        {"parameters",
         {{"before_first", "zero"},
          {"after_last", "constant"},
          {"model",
           {{{"epoch", "2016-01-01T00:00:00Z"}, {"scale_factor", 0.5}},
            {{"epoch", "2017-01-01T00:00:00Z"}, {"scale_factor", 1.0}},
            {{"epoch", "2017-01-01T00:00:00Z"}, {"scale_factor", 2.0}},
            {{"epoch", "2018-01-01T00:00:00Z"}, {"scale_factor", 1.0}}}}}}};

    j["components"].push_back(j["components"][0]);
    j["components"][IDX_EXPONENTIAL]["time_function"] = {
        {"type", "exponential"},
        {"parameters",
         {
             {"reference_epoch", "2000-01-01T00:00:00Z"},
             {"end_epoch", "2001-01-01T00:00:00Z"},
             {"relaxation_constant", 2.0},
             {"before_scale_factor", 0.0},
             {"initial_scale_factor", 1.0},
             {"final_scale_factor", 3.0},
         }},
    };

    return j;
}

// ---------------------------------------------------------------------------

TEST(defmodel, basic) {
    EXPECT_THROW(MasterFile::parse("foo"), ParsingException);
    EXPECT_THROW(MasterFile::parse("null"), ParsingException);
    EXPECT_THROW(MasterFile::parse("{}"), ParsingException);

    const auto jMinValid(getMinValidContent());
    {
        auto mf = MasterFile::parse(jMinValid.dump());
        EXPECT_EQ(mf->fileType(), "GeoTIFF");
        EXPECT_EQ(mf->formatVersion(), "1.0");
        EXPECT_EQ(mf->sourceCRS(), "EPSG:4959");
        EXPECT_EQ(mf->targetCRS(), "EPSG:7907");
        EXPECT_EQ(mf->definitionCRS(), "EPSG:4959");
        EXPECT_EQ(mf->extent().minx(), modelMinX);
        EXPECT_EQ(mf->extent().miny(), modelMinY);
        EXPECT_EQ(mf->extent().maxx(), modelMaxX);
        EXPECT_EQ(mf->extent().maxy(), modelMaxY);
        EXPECT_EQ(mf->timeExtent().first.toString(), "1900-01-01T00:00:00Z");
        EXPECT_EQ(mf->timeExtent().last.toString(), "2050-01-01T00:00:00Z");
    }

    // Check that removing one of each required key causes an exception
    for (const auto &kv : jMinValid.items()) {
        json jcopy(jMinValid);
        jcopy.erase(kv.key());
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jMinValid);
        jcopy["definition_crs"] = "EPSG:4326";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jMinValid);
        jcopy["file_type"] = 1;
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jMinValid);
        jcopy["extent"].erase("type");
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jMinValid);
        jcopy["extent"].erase("parameters");
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jMinValid);
        jcopy["extent"]["parameters"].clear();
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jMinValid);
        jcopy["extent"]["parameters"].erase("bbox");
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jMinValid);
        jcopy["extent"]["parameters"]["bbox"] = "foo";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jMinValid);
        jcopy["extent"]["parameters"]["bbox"] = {0, 1, 2};
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jMinValid);
        jcopy["extent"]["parameters"]["bbox"] = {0, 1, 2, "foo"};
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jMinValid);
        jcopy["time_extent"] = "foo";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jMinValid);
        jcopy["time_extent"].erase("first");
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jMinValid);
        jcopy["time_extent"].erase("last");
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }
}

// ---------------------------------------------------------------------------

TEST(defmodel, full) {
    const auto jFullValid(getFullValidContent());
    auto mf = MasterFile::parse(jFullValid.dump());
    EXPECT_EQ(mf->name(), "name");
    EXPECT_EQ(mf->version(), "version");
    EXPECT_EQ(mf->publicationDate(), "2018-07-01T00:00:00Z");
    EXPECT_EQ(mf->license(), "license");
    EXPECT_EQ(mf->description(), "description");
    EXPECT_EQ(mf->authority().name, "authority_name");
    EXPECT_EQ(mf->authority().url, "authority_url");
    EXPECT_EQ(mf->authority().address, "authority_address");
    EXPECT_EQ(mf->authority().email, "authority_email");
    EXPECT_EQ(mf->links().size(), 1U);
    EXPECT_EQ(mf->links()[0].href, "href");
    EXPECT_EQ(mf->links()[0].rel, "rel");
    EXPECT_EQ(mf->links()[0].type, "type");
    EXPECT_EQ(mf->links()[0].title, "title");
    EXPECT_EQ(mf->referenceEpoch(), "2000-01-01T00:00:00Z");
    EXPECT_EQ(mf->uncertaintyReferenceEpoch(), "2018-12-15T00:00:00Z");
    EXPECT_EQ(mf->horizontalOffsetUnit(), "metre");
    EXPECT_EQ(mf->verticalOffsetUnit(), "metre");
    EXPECT_EQ(mf->horizontalUncertaintyType(), "circular 95% confidence limit");
    EXPECT_EQ(mf->horizontalUncertaintyUnit(), "metre");
    EXPECT_EQ(mf->verticalUncertaintyType(), "95% confidence limit");
    EXPECT_EQ(mf->verticalUncertaintyUnit(), "metre");
    EXPECT_EQ(mf->horizontalOffsetMethod(), "addition");
    ASSERT_EQ(mf->components().size(), 6U);
    {
        const auto &comp = mf->components()[IDX_CONSTANT];
        EXPECT_EQ(comp.description(), "description");
        EXPECT_EQ(comp.displacementType(), "horizontal");
        EXPECT_EQ(comp.uncertaintyType(), "none");
        EXPECT_EQ(comp.horizontalUncertainty(), 0.01);
        EXPECT_EQ(comp.verticalUncertainty(), 0.02);
        EXPECT_EQ(comp.extent().minx(), modelMinX);
        EXPECT_EQ(comp.extent().miny(), modelMinY);
        EXPECT_EQ(comp.extent().maxx(), modelMaxX);
        EXPECT_EQ(comp.extent().maxy(), modelMaxY);
        EXPECT_EQ(comp.spatialModel().type, "GeoTIFF");
        EXPECT_EQ(comp.spatialModel().interpolationMethod, "bilinear");
        EXPECT_EQ(comp.spatialModel().filename, "nzgd2000-ndm-grid02.tif");
        EXPECT_EQ(comp.spatialModel().md5Checksum,
                  "49fce8ab267be2c8d00d43683060a032");
        ASSERT_NE(comp.timeFunction(), nullptr);
        ASSERT_EQ(comp.timeFunction()->type, "constant");
    }
    {
        const auto &comp = mf->components()[IDX_VELOCITY];
        ASSERT_NE(comp.timeFunction(), nullptr);
        ASSERT_EQ(comp.timeFunction()->type, "velocity");
        const auto velocity =
            static_cast<const Component::VelocityTimeFunction *>(
                comp.timeFunction());
        EXPECT_EQ(velocity->referenceEpoch.toString(), "2000-01-01T00:00:00Z");
    }
    {
        const auto &comp = mf->components()[IDX_STEP];
        ASSERT_NE(comp.timeFunction(), nullptr);
        ASSERT_EQ(comp.timeFunction()->type, "step");
        const auto step = static_cast<const Component::StepTimeFunction *>(
            comp.timeFunction());
        EXPECT_EQ(step->stepEpoch.toString(), "2000-01-01T00:00:00Z");
    }
    {
        const auto &comp = mf->components()[IDX_REVERSE_STEP];
        ASSERT_NE(comp.timeFunction(), nullptr);
        ASSERT_EQ(comp.timeFunction()->type, "reverse_step");
        const auto step =
            static_cast<const Component::ReverseStepTimeFunction *>(
                comp.timeFunction());
        EXPECT_EQ(step->stepEpoch.toString(), "2000-01-01T00:00:00Z");
    }
    {
        const auto &comp = mf->components()[IDX_PIECEWISE];
        ASSERT_NE(comp.timeFunction(), nullptr);
        ASSERT_EQ(comp.timeFunction()->type, "piecewise");
        const auto piecewise =
            static_cast<const Component::PiecewiseTimeFunction *>(
                comp.timeFunction());
        EXPECT_EQ(piecewise->beforeFirst, "zero");
        EXPECT_EQ(piecewise->afterLast, "constant");
        EXPECT_EQ(piecewise->model.size(), 4U);
        EXPECT_EQ(piecewise->model[0].epoch.toString(), "2016-01-01T00:00:00Z");
        EXPECT_EQ(piecewise->model[0].scaleFactor, 0.5);
    }
    {
        const auto &comp = mf->components()[IDX_EXPONENTIAL];
        ASSERT_NE(comp.timeFunction(), nullptr);
        ASSERT_EQ(comp.timeFunction()->type, "exponential");
        const auto exponential =
            static_cast<const Component::ExponentialTimeFunction *>(
                comp.timeFunction());
        EXPECT_EQ(exponential->referenceEpoch.toString(),
                  "2000-01-01T00:00:00Z");
        EXPECT_EQ(exponential->endEpoch.toString(), "2001-01-01T00:00:00Z");
        EXPECT_EQ(exponential->relaxationConstant, 2.0);
        EXPECT_EQ(exponential->beforeScaleFactor, 0.0);
        EXPECT_EQ(exponential->initialScaleFactor, 1.0);
        EXPECT_EQ(exponential->finalScaleFactor, 3.0);
    }
}

// ---------------------------------------------------------------------------

TEST(defmodel, error_cases) {
    const auto jFullValid(getFullValidContent());
    {
        json jcopy(jFullValid);
        jcopy["horizontal_offset_method"] = "unsupported";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jFullValid);
        jcopy["horizontal_offset_unit"] = "unsupported";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jFullValid);
        jcopy["vertical_offset_unit"] = "unsupported";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_CONSTANT]["spatial_model"]
             ["interpolation_method"] = "unsupported";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_CONSTANT]["displacement_type"] = "unsupported";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["model"] = "foo";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["before_first"] = "illegal";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["after_last"] = "illegal";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    {
        json jcopy(jFullValid);
        jcopy["components"][0]["time_function"]["type"] = "unknown";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    // Unsupported combination
    {
        json jcopy(jFullValid);
        jcopy["horizontal_offset_method"] = "geocentric";
        EXPECT_NO_THROW(MasterFile::parse(jcopy.dump()));
    }
    {
        json jcopy(jFullValid);
        jcopy["horizontal_offset_unit"] = "degree";
        EXPECT_NO_THROW(MasterFile::parse(jcopy.dump()));
    }
    {
        json jcopy(jFullValid);
        jcopy["horizontal_offset_method"] = "geocentric";
        jcopy["horizontal_offset_unit"] = "degree";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }

    // Unsupported combination
    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_VELOCITY]["spatial_model"]
             ["interpolation_method"] = "geocentric_bilinear";
        EXPECT_NO_THROW(MasterFile::parse(jcopy.dump()));
    }
    {
        json jcopy(jFullValid);
        jcopy["horizontal_offset_unit"] = "degree";
        EXPECT_NO_THROW(MasterFile::parse(jcopy.dump()));
    }
    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_VELOCITY]["spatial_model"]
             ["interpolation_method"] = "geocentric_bilinear";
        jcopy["horizontal_offset_unit"] = "degree";
        EXPECT_THROW(MasterFile::parse(jcopy.dump()), ParsingException);
    }
}

// ---------------------------------------------------------------------------

TEST(defmodel, ISO8601ToDecimalYear) {
    EXPECT_EQ(ISO8601ToDecimalYear("2000-01-01T00:00:00Z"), 2000.0);
    EXPECT_EQ(ISO8601ToDecimalYear("2000-02-29T12:00:00Z"),
              2000.0 + ((31 + 28) * 86400. + 12 * 3600) / (366 * 86400));
    EXPECT_EQ(ISO8601ToDecimalYear("2000-12-31T23:59:59Z"),
              2000.0 + (366 * 86400 - 1.) / (366 * 86400));
    EXPECT_EQ(ISO8601ToDecimalYear("2001-01-01T00:00:00Z"), 2001.0);
    EXPECT_EQ(ISO8601ToDecimalYear("2001-12-31T23:59:59Z"),
              2001.0 + (365 * 86400 - 1.) / (365 * 86400));
    EXPECT_THROW(ISO8601ToDecimalYear(""), ParsingException);
    EXPECT_THROW(ISO8601ToDecimalYear("0000-01-01T00:00:00Z"),
                 ParsingException);
    EXPECT_THROW(ISO8601ToDecimalYear("2001-02-29T00:00:00Z"),
                 ParsingException);
    EXPECT_THROW(ISO8601ToDecimalYear("2000-13-01T00:00:00Z"),
                 ParsingException);
    EXPECT_THROW(ISO8601ToDecimalYear("2000-01-32T00:00:00Z"),
                 ParsingException);
    EXPECT_THROW(ISO8601ToDecimalYear("2000-01-01T24:00:00Z"),
                 ParsingException);
    EXPECT_THROW(ISO8601ToDecimalYear("2000-01-01T00:60:00Z"),
                 ParsingException);
    EXPECT_THROW(ISO8601ToDecimalYear("2000-01-01T00:00:61Z"),
                 ParsingException);
}

// ---------------------------------------------------------------------------

TEST(defmodel, evaluate_constant) {
    const auto jFullValid(getFullValidContent());
    const auto mf = MasterFile::parse(jFullValid.dump());
    const auto &comp = mf->components()[IDX_CONSTANT];
    EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.0), 1.0);
    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.0), 1.0);
    EXPECT_EQ(comp.timeFunction()->evaluateAt(2001.0), 1.0);
}

// ---------------------------------------------------------------------------

TEST(defmodel, evaluate_velocity) {
    const auto jFullValid(getFullValidContent());
    const auto mf = MasterFile::parse(jFullValid.dump());
    const auto &comp = mf->components()[IDX_VELOCITY];
    EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.0), -1.0);
    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.0), 0.0);
    EXPECT_EQ(comp.timeFunction()->evaluateAt(2001.0), 1.0);
}

// ---------------------------------------------------------------------------

TEST(defmodel, evaluate_step) {
    const auto jFullValid(getFullValidContent());
    const auto mf = MasterFile::parse(jFullValid.dump());
    const auto &comp = mf->components()[IDX_STEP];
    EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.99), 0.0);
    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.00), 1.0);
    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.01), 1.0);
}

// ---------------------------------------------------------------------------

TEST(defmodel, evaluate_reverse_step) {
    const auto jFullValid(getFullValidContent());
    const auto mf = MasterFile::parse(jFullValid.dump());
    const auto &comp = mf->components()[IDX_REVERSE_STEP];
    EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.99), -1.0);
    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.00), 0.0);
    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.01), 0.0);
}

// ---------------------------------------------------------------------------

TEST(defmodel, evaluate_piecewise) {
    const auto jFullValid(getFullValidContent());
    {
        const auto mf = MasterFile::parse(jFullValid.dump());
        const auto &comp = mf->components()[IDX_PIECEWISE];
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.99), 0.0);
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2016.00), 0.5);
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2016.5), 0.75);
        EXPECT_NEAR(comp.timeFunction()->evaluateAt(2017 - 1e-9), 1.0, 1e-9);
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2017.0), 2.0);
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2017.5), 1.5);
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2018.0), 1.0);
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2019.0), 1.0);
    }

    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["before_first"] = "zero";
        const auto mf = MasterFile::parse(jcopy.dump());
        const auto &comp = mf->components()[IDX_PIECEWISE];
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.0);
    }

    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["before_first"] = "constant";
        const auto mf = MasterFile::parse(jcopy.dump());
        const auto &comp = mf->components()[IDX_PIECEWISE];
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.5);
    }

    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["before_first"] = "linear";
        const auto mf = MasterFile::parse(jcopy.dump());
        const auto &comp = mf->components()[IDX_PIECEWISE];
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.25);
    }

    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["after_last"] = "zero";
        const auto mf = MasterFile::parse(jcopy.dump());
        const auto &comp = mf->components()[IDX_PIECEWISE];
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2018.5), 0.0);
    }

    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["after_last"] = "constant";
        const auto mf = MasterFile::parse(jcopy.dump());
        const auto &comp = mf->components()[IDX_PIECEWISE];
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2018.5), 1.0);
    }

    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["after_last"] = "linear";
        const auto mf = MasterFile::parse(jcopy.dump());
        const auto &comp = mf->components()[IDX_PIECEWISE];
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2018.5), 0.5);
    }

    // No epoch
    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["model"]
                 .clear();
        const auto mf = MasterFile::parse(jcopy.dump());
        const auto &comp = mf->components()[IDX_PIECEWISE];
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.0);
    }

    // Just one epoch
    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["model"] = {
                 {{"epoch", "2016-01-01T00:00:00Z"}, {"scale_factor", 0.5}}};
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["before_first"] = "linear";
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["after_last"] = "linear";
        const auto mf = MasterFile::parse(jcopy.dump());
        const auto &comp = mf->components()[IDX_PIECEWISE];
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.5);
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2016.5), 0.5);
    }

    // Two identical epochs
    {
        json jcopy(jFullValid);
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["model"] = {
                 {{"epoch", "2016-01-01T00:00:00Z"}, {"scale_factor", 0.5}},
                 {{"epoch", "2016-01-01T00:00:00Z"}, {"scale_factor", 1.0}}};
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["before_first"] = "linear";
        jcopy["components"][IDX_PIECEWISE]["time_function"]["parameters"]
             ["after_last"] = "linear";
        const auto mf = MasterFile::parse(jcopy.dump());
        const auto &comp = mf->components()[IDX_PIECEWISE];
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2015.5), 0.5);
        EXPECT_EQ(comp.timeFunction()->evaluateAt(2016.5), 1.0);
    }
}

// ---------------------------------------------------------------------------

TEST(defmodel, evaluate_exponential) {
    const auto jFullValid(getFullValidContent());
    const auto mf = MasterFile::parse(jFullValid.dump());
    const auto &comp = mf->components()[IDX_EXPONENTIAL];
    EXPECT_EQ(comp.timeFunction()->evaluateAt(1999.99), 0.0);
    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.00), 1.0);
    EXPECT_EQ(comp.timeFunction()->evaluateAt(2000.50),
              1.0 + (3.0 - 1.0) * (1.0 - std::exp(-(2000.50 - 2000.00) / 2.0)));
    EXPECT_EQ(comp.timeFunction()->evaluateAt(2001.00),
              1.0 + (3.0 - 1.0) * (1.0 - std::exp(-(2001.00 - 2000.00) / 2.0)));
    EXPECT_EQ(comp.timeFunction()->evaluateAt(2002.00),
              1.0 + (3.0 - 1.0) * (1.0 - std::exp(-(2001.00 - 2000.00) / 2.0)));
}

// ---------------------------------------------------------------------------

inline double RadToDeg(double d) { return d / DEG_TO_RAD_CONSTANT; }

// ---------------------------------------------------------------------------

TEST(defmodel, evaluator_horizontal_unit_degree) {
    json j(getMinValidContent());
    j["horizontal_offset_method"] = "addition";
    j["horizontal_offset_unit"] = "degree";
    constexpr double tFactor = 0.5;
    constexpr double gridMinX = 160;
    constexpr double gridMinY = -50;
    constexpr double gridMaxX = 190;
    constexpr double gridMaxY = -30;
    j["components"] = {
        {{"displacement_type", "horizontal"},
         {"uncertainty_type", "none"},
         {"extent",
          {{"type", "bbox"},
           {"parameters",
            {{"bbox", {gridMinX, gridMinY, gridMaxX, gridMaxY}}}}}},
         {"spatial_model",
          {
              {"type", "GeoTIFF"},
              {"interpolation_method", "bilinear"},
              {"filename", "bla.tif"},
          }},
         {"time_function",
          {{"type", "piecewise"},
           {"parameters",
            {{"before_first", "zero"},
             {"after_last", "zero"},
             {"model",
              {{{"epoch", "2010-01-01T00:00:00Z"}, {"scale_factor", tFactor}},
               {{"epoch", "2020-01-01T00:00:00Z"},
                {"scale_factor", tFactor}}}}}}}}}};

    constexpr int iQueriedX = 1;
    constexpr int iQueriedY = 3;
    constexpr double longOffsetQueriedX = 0.01;
    constexpr double longOffsetQueriedXp1 = 0.02;
    constexpr double latOffsetQueriedY = 0.03;
    constexpr double latOffsetQueriedYp1 = 0.04;
    constexpr double zOffsetQueriedXY = 10.;
    constexpr double zOffsetQueriedXp1Y = 11.;
    constexpr double zOffsetQueriedXYp1 = 11.;
    constexpr double zOffsetQueriedXp1Yp1 = 12.;
    constexpr double gridResX = 2;
    constexpr double gridResY = 0.5;

    struct Grid : public GridPrototype {
        bool getLongLatOffset(int ix, int iy, double &longOffsetRadian,
                              double &latOffsetRadian) const {
            if (ix == iQueriedX) {
                longOffsetRadian = DegToRad(longOffsetQueriedX);
            } else if (ix == iQueriedX + 1) {
                longOffsetRadian = DegToRad(longOffsetQueriedXp1);
            } else {
                return false;
            }
            if (iy == iQueriedY) {
                latOffsetRadian = DegToRad(latOffsetQueriedY);
            } else if (iy == iQueriedY + 1) {
                latOffsetRadian = DegToRad(latOffsetQueriedYp1);
            } else {
                return false;
            }
            return true;
        }

        bool getZOffset(int ix, int iy, double &zOffset) const {
            if (ix == iQueriedX && iy == iQueriedY) {
                zOffset = zOffsetQueriedXY;
            } else if (ix == iQueriedX + 1 && iy == iQueriedY) {
                zOffset = zOffsetQueriedXp1Y;
            } else if (ix == iQueriedX && iy == iQueriedY + 1) {
                zOffset = zOffsetQueriedXYp1;
            } else if (ix == iQueriedX + 1 && iy == iQueriedY + 1) {
                zOffset = zOffsetQueriedXp1Yp1;
            } else {
                return false;
            }
            return true;
        }

        bool getLongLatZOffset(int ix, int iy, double &longOffsetRadian,
                               double &latOffsetRadian, double &zOffset) const {
            return getLongLatOffset(ix, iy, longOffsetRadian,
                                    latOffsetRadian) &&
                   getZOffset(ix, iy, zOffset);
        }

#ifdef DEBUG_DEFMODEL
        std::string name() const { return std::string(); }
#endif
    };

    struct GridSet : public GridSetPrototype<Grid> {

        Grid grid{};

        GridSet() {
            grid.minx = DegToRad(gridMinX);
            grid.miny = DegToRad(gridMinY);
            grid.resx = DegToRad(gridResX);
            grid.resy = DegToRad(gridResY);
            grid.width =
                1 + static_cast<int>(0.5 + (gridMaxX - gridMinX) / gridResX);
            grid.height =
                1 + static_cast<int>(0.5 + (gridMaxY - gridMinY) / gridResY);
        }

        const Grid *gridAt(double /*x */, double /* y */) { return &grid; }
    };

    struct EvaluatorIface : public EvaluatorIfacePrototype<Grid, GridSet> {
        std::unique_ptr<GridSet> open(const std::string &filename) {
            if (filename != "bla.tif")
                return nullptr;
            return std::unique_ptr<GridSet>(new GridSet());
        }

        bool isGeographicCRS(const std::string & /* crsDef */) { return true; }

#ifdef DEBUG_DEFMODEL
        void log(const std::string & /* msg */) {}
#endif
    };

    EvaluatorIface iface;

    Evaluator<Grid, GridSet, EvaluatorIface> eval(MasterFile::parse(j.dump()),
                                                  iface, 1, 1);
    double newLong;
    double newLat;
    double newZ;
    constexpr double tValid = 2018;
    constexpr double EPS = 1e-9;
    constexpr double zVal = 100;

    // Query on exact grid intersection
    {
        const double longitude = gridMinX + iQueriedX * gridResX;
        const double lat = gridMinY + iQueriedY * gridResY;
        EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                 zVal, tValid, newLong, newLat, newZ));
        EXPECT_NEAR(RadToDeg(newLong), longitude + tFactor * longOffsetQueriedX,
                    EPS);
        EXPECT_NEAR(RadToDeg(newLat), lat + tFactor * latOffsetQueriedY, EPS);
        EXPECT_EQ(newZ, zVal);
    }

    // Query between grid points
    {
        constexpr double alphaX = 0.25;
        constexpr double alphaY = 0.125;
        const double longitude =
            gridMinX + iQueriedX * gridResX + alphaX * gridResX;
        const double lat = gridMinY + iQueriedY * gridResY + alphaY * gridResY;
        EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                 zVal, tValid, newLong, newLat, newZ));
        EXPECT_NEAR(RadToDeg(newLong),
                    longitude + tFactor * (longOffsetQueriedX +
                                           alphaX * (longOffsetQueriedXp1 -
                                                     longOffsetQueriedX)),
                    EPS);
        EXPECT_NEAR(
            RadToDeg(newLat),
            lat + tFactor * (latOffsetQueriedY + alphaY * (latOffsetQueriedYp1 -
                                                           latOffsetQueriedY)),
            EPS);
        EXPECT_EQ(newZ, zVal);
    }

    // Longitude < model min
    {
        const double longitude = modelMinX - 1e-1;
        const double lat = gridMinY + iQueriedY * gridResY;
        EXPECT_FALSE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                  zVal, tValid, newLong, newLat, newZ));
    }

    // Longitude > model max
    {
        const double longitude = modelMaxX + 1e-1;
        const double lat = gridMinY + iQueriedY * gridResY;
        EXPECT_FALSE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                  zVal, tValid, newLong, newLat, newZ));
    }

    // Latitude < model min
    {
        const double longitude = gridMinX + iQueriedX * gridResX;
        const double lat = modelMinY - 1e-1;
        EXPECT_FALSE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                  zVal, tValid, newLong, newLat, newZ));
    }

    // Latitude > model max
    {
        const double longitude = gridMinX + iQueriedX * gridResX;
        const double lat = modelMaxY + 1e-1;
        EXPECT_FALSE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                  zVal, tValid, newLong, newLat, newZ));
    }

    // Before timeExtent.first
    {
        const double longitude = gridMinX + iQueriedX * gridResX;
        const double lat = gridMinY + iQueriedY * gridResY;
        EXPECT_FALSE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                  zVal, 1000, newLong, newLat, newZ));
    }

    // After timeExtent.last
    {
        const double longitude = gridMinX + iQueriedX * gridResX;
        const double lat = gridMinY + iQueriedY * gridResY;
        EXPECT_FALSE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                  zVal, 3000, newLong, newLat, newZ));
    }

    // Longitude < grid min
    {
        const double longitude = gridMinX - 1e-1;
        const double lat = gridMinY + iQueriedY * gridResY;
        EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                 zVal, tValid, newLong, newLat, newZ));
        EXPECT_NEAR(RadToDeg(newLong), longitude, EPS);
        EXPECT_NEAR(RadToDeg(newLat), lat, EPS);
        EXPECT_EQ(newZ, zVal);
    }

    // Longitude > grid max
    {
        const double longitude = gridMaxX + 1e-1;
        const double lat = gridMinY + iQueriedY * gridResY;
        EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                 zVal, tValid, newLong, newLat, newZ));
        EXPECT_NEAR(RadToDeg(newLong), longitude, EPS);
        EXPECT_NEAR(RadToDeg(newLat), lat, EPS);
        EXPECT_EQ(newZ, zVal);
    }

    // Latitude < grid min
    {
        const double longitude = gridMinX + iQueriedX * gridResX;
        const double lat = gridMinY - 1e-1;
        EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                 zVal, tValid, newLong, newLat, newZ));
        EXPECT_NEAR(RadToDeg(newLong), longitude, EPS);
        EXPECT_NEAR(RadToDeg(newLat), lat, EPS);
        EXPECT_EQ(newZ, zVal);
    }

    // Latitude > grid max
    {
        const double longitude = gridMinX + iQueriedX * gridResX;
        const double lat = gridMaxY + 1e-1;
        EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                 zVal, tValid, newLong, newLat, newZ));
        EXPECT_NEAR(RadToDeg(newLong), longitude, EPS);
        EXPECT_NEAR(RadToDeg(newLat), lat, EPS);
        EXPECT_EQ(newZ, zVal);
    }

    // Time function values to zero
    {
        const double longitude = gridMinX + iQueriedX * gridResX;
        const double lat = gridMinY + iQueriedY * gridResY;
        EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                 zVal, 2000, newLong, newLat, newZ));
        EXPECT_NEAR(RadToDeg(newLong), longitude, EPS);
        EXPECT_NEAR(RadToDeg(newLat), lat, EPS);
        EXPECT_EQ(newZ, zVal);
    }

    // Test vertical
    j["components"][0]["displacement_type"] = "vertical";
    j["vertical_offset_unit"] = "metre";
    Evaluator<Grid, GridSet, EvaluatorIface> evalVertical(
        MasterFile::parse(j.dump()), iface, 1, 1);
    {
        constexpr double alphaX = 0.25;
        constexpr double alphaY = 0.125;
        const double longitude =
            gridMinX + iQueriedX * gridResX + alphaX * gridResX;
        const double lat = gridMinY + iQueriedY * gridResY + alphaY * gridResY;
        EXPECT_TRUE(evalVertical.forward(iface, DegToRad(longitude),
                                         DegToRad(lat), zVal, tValid, newLong,
                                         newLat, newZ));
        EXPECT_NEAR(RadToDeg(newLong), longitude, EPS);
        EXPECT_NEAR(RadToDeg(newLat), lat, EPS);

        const double zBottom =
            zOffsetQueriedXY + alphaX * (zOffsetQueriedXp1Y - zOffsetQueriedXY);
        const double zTop =
            zOffsetQueriedXYp1 +
            alphaX * (zOffsetQueriedXp1Yp1 - zOffsetQueriedXYp1);
        EXPECT_NEAR(
            newZ, zVal + tFactor * (zBottom + alphaY * (zTop - zBottom)), EPS);
    }

    // Test 3d
    j["components"][0]["displacement_type"] = "3d";
    j["vertical_offset_unit"] = "metre";
    Evaluator<Grid, GridSet, EvaluatorIface> eval3d(MasterFile::parse(j.dump()),
                                                    iface, 1, 1);
    {
        constexpr double alphaX = 0.25;
        constexpr double alphaY = 0.125;
        const double longitude =
            gridMinX + iQueriedX * gridResX + alphaX * gridResX;
        const double lat = gridMinY + iQueriedY * gridResY + alphaY * gridResY;
        EXPECT_TRUE(eval3d.forward(iface, DegToRad(longitude), DegToRad(lat),
                                   zVal, tValid, newLong, newLat, newZ));
        EXPECT_NEAR(RadToDeg(newLong),
                    longitude + tFactor * (longOffsetQueriedX +
                                           alphaX * (longOffsetQueriedXp1 -
                                                     longOffsetQueriedX)),
                    EPS);
        EXPECT_NEAR(
            RadToDeg(newLat),
            lat + tFactor * (latOffsetQueriedY + alphaY * (latOffsetQueriedYp1 -
                                                           latOffsetQueriedY)),
            EPS);

        const double zBottom =
            zOffsetQueriedXY + alphaX * (zOffsetQueriedXp1Y - zOffsetQueriedXY);
        const double zTop =
            zOffsetQueriedXYp1 +
            alphaX * (zOffsetQueriedXp1Yp1 - zOffsetQueriedXYp1);
        EXPECT_NEAR(
            newZ, zVal + tFactor * (zBottom + alphaY * (zTop - zBottom)), EPS);
    }
}

// ---------------------------------------------------------------------------

inline void DeltaLongLatToEastingNorthing(double phi, double dlam, double dphi,
                                          double a, double b, double &de,
                                          double &dn) {
    const double sinphi = sin(phi);
    const double cosphi = cos(phi);
    const double a2 = a * a;
    const double b2 = b * b;
    const double X = a2 * (cosphi * cosphi) + b2 * (sinphi * sinphi);
    const double sqrtX = sqrt(X);
    de = dlam * (a2 * cosphi) / sqrtX;
    dn = dphi * a2 * b2 / (sqrtX * X);
}

// ---------------------------------------------------------------------------

TEST(defmodel, evaluator_horizontal_unit_metre) {

    json j(getMinValidContent());
    j["horizontal_offset_method"] = "addition";
    j["horizontal_offset_unit"] = "metre";
    j["vertical_offset_unit"] = "metre";
    constexpr double tFactor = 0.5;
    constexpr double gridMinX = 165.8;
    constexpr double gridMinY = -37.5;
    constexpr double gridMaxX = 166.2;
    constexpr double gridMaxY = -37.2;
    constexpr double gridResX = gridMaxX - gridMinX;
    constexpr double gridResY = gridMaxY - gridMinY;
    constexpr int extraPointX = 1;
    constexpr int extraPointY = 1;

    j["components"] = {
        {{"displacement_type", "horizontal"},
         {"uncertainty_type", "none"},
         {"extent",
          {{"type", "bbox"},
           {"parameters",
            {{"bbox",
              {gridMinX - extraPointX * gridResX,
               gridMinY - extraPointY * gridResY, gridMaxX, gridMaxY}}}}}},
         {"spatial_model",
          {
              {"type", "GeoTIFF"},
              {"interpolation_method", "XXXXXXX"},
              {"filename", "bla.tif"},
          }},
         {"time_function",
          {{"type", "piecewise"},
           {"parameters",
            {{"before_first", "zero"},
             {"after_last", "zero"},
             {"model",
              {{{"epoch", "2010-01-01T00:00:00Z"}, {"scale_factor", tFactor}},
               {{"epoch", "2020-01-01T00:00:00Z"},
                {"scale_factor", tFactor}}}}}}}}}};

    struct Grid : public GridPrototype {
        bool getEastingNorthingOffset(int ix, int iy, double &eastingOffset,
                                      double &northingOffset) const {
            ix -= extraPointX;
            iy -= extraPointY;
            if (ix == -1)
                ix = 0;
            if (iy == -1)
                iy = 0;
            if (ix == 0 && iy == 0) {
                eastingOffset = 0.4f;
                northingOffset = -0.2f;
            } else if (ix == 1 && iy == 0) {
                eastingOffset = 0.5f;
                northingOffset = -0.25f;
            } else if (ix == 0 && iy == 1) {
                eastingOffset = 0.8f;
                northingOffset = -0.4f;
            } else if (ix == 1 && iy == 1) {
                eastingOffset = 1.f;
                northingOffset = -0.3f;
            } else {
                return false;
            }
            return true;
        }

        bool getZOffset(int ix, int iy, double &zOffset) const {
            ix -= extraPointX;
            iy -= extraPointY;
            if (ix == -1)
                ix = 0;
            if (iy == -1)
                iy = 0;
            if (ix == 0 && iy == 0) {
                zOffset = 0.84f;
            } else if (ix == 1 && iy == 0) {
                zOffset = 0.75f;
            } else if (ix == 0 && iy == 1) {
                zOffset = 0.36f;
            } else if (ix == 1 && iy == 1) {
                zOffset = 0.f;
            } else {
                return false;
            }
            return true;
        }

        bool getEastingNorthingZOffset(int ix, int iy, double &eastingOffset,
                                       double &northingOffset,
                                       double &zOffset) const {
            return getEastingNorthingOffset(ix, iy, eastingOffset,
                                            northingOffset) &&
                   getZOffset(ix, iy, zOffset);
        }

#ifdef DEBUG_DEFMODEL
        std::string name() const { return std::string(); }
#endif
    };

    struct GridSet : public GridSetPrototype<Grid> {

        Grid grid{};

        GridSet() {
            grid.minx = DegToRad(gridMinX - extraPointX * gridResX);
            grid.miny = DegToRad(gridMinY - extraPointY * gridResY);
            grid.resx = DegToRad(gridResX);
            grid.resy = DegToRad(gridResY);
            grid.width = 2 + extraPointX;
            grid.height = 2 + extraPointY;
        }

        const Grid *gridAt(double /*x */, double /* y */) { return &grid; }
    };

    struct EvaluatorIface : public EvaluatorIfacePrototype<Grid, GridSet> {
        std::unique_ptr<GridSet> open(const std::string &filename) {
            if (filename != "bla.tif")
                return nullptr;
            return std::unique_ptr<GridSet>(new GridSet());
        }

        bool isGeographicCRS(const std::string & /* crsDef */) { return true; }

#ifdef DEBUG_DEFMODEL
        void log(const std::string & /* msg */) {}
#endif

        void geographicToGeocentric(double lam, double phi, double height,
                                    double a, double /*b*/, double es,
                                    double &X, double &Y, double &Z) {
            PJ_CONTEXT *ctx = proj_context_create();
            PJ *cart = proj_create(
                ctx,
                ("+proj=cart +a=" + osgeo::proj::internal::toString(a, 18) +
                 " +es=" + osgeo::proj::internal::toString(es, 18))
                    .c_str());
            PJ_LPZ lpz;
            lpz.lam = lam;
            lpz.phi = phi;
            lpz.z = height;
            PJ_XYZ xyz = cart->fwd3d(lpz, cart);
            X = xyz.x;
            Y = xyz.y;
            Z = xyz.z;
            proj_destroy(cart);
            proj_context_destroy(ctx);
        }

        void geocentricToGeographic(double X, double Y, double Z, double a,
                                    double /*b*/, double es, double &lam,
                                    double &phi, double &height) {
            PJ_CONTEXT *ctx = proj_context_create();
            PJ *cart = proj_create(
                ctx,
                ("+proj=cart +a=" + osgeo::proj::internal::toString(a, 18) +
                 " +es=" + osgeo::proj::internal::toString(es, 18))
                    .c_str());
            PJ_XYZ xyz;
            xyz.x = X;
            xyz.y = Y;
            xyz.z = Z;
            PJ_LPZ lpz = cart->inv3d(xyz, cart);
            lam = lpz.lam;
            phi = lpz.phi;
            height = lpz.z;
            proj_destroy(cart);
            proj_context_destroy(ctx);
        }
    };
    EvaluatorIface iface;

    constexpr double a = 6378137;
    constexpr double b = 6356752.314140;
    constexpr double tValid = 2018;
    constexpr double zVal = 100;

    const struct {
        double longitude;
        double lat;
        double expected_de;
        double expected_dn;
        double expected_dz;
        const char *displacement_type;
        const char *interpolation_method;
    } testPoints[] = {
        {gridMinX - extraPointX * gridResX - 1e-11,
         gridMinY - extraPointY * gridResY - 1e-11, 0.4, -0.2, 0, "horizontal",
         "bilinear"},
        {gridMinX, gridMinY, 0.4, -0.2, 0, "horizontal", "bilinear"},
        {gridMaxX, gridMinY, 0.5, -0.25, 0, "horizontal", "bilinear"},
        {gridMinX, gridMaxY, 0.8, -0.4, 0, "horizontal", "bilinear"},
        {gridMaxX, gridMaxY, 1, -0.3, 0, "horizontal", "bilinear"},
        {gridMaxX + 1e-11, gridMaxY + 1e-11, 1, -0.3, 0, "horizontal",
         "bilinear"},
        {165.9, -37.3, 0.70833334, -0.32083334, 0, "horizontal", "bilinear"},
        {165.9, -37.3, 0.70833334, -0.32083334, 0.4525, "3d", "bilinear"},

        {gridMinX, gridMinY, 0.4, -0.2, 0, "horizontal", "geocentric_bilinear"},
        {gridMaxX, gridMinY, 0.5, -0.25, 0, "horizontal",
         "geocentric_bilinear"},
        {gridMinX, gridMaxY, 0.8, -0.4, 0, "horizontal", "geocentric_bilinear"},
        {gridMaxX, gridMaxY, 1, -0.3, 0, "horizontal", "geocentric_bilinear"},
        {165.9, -37.3, 0.7083692044608846, -0.3209642339711405, 0, "horizontal",
         "geocentric_bilinear"},
        {165.9, -37.3, 0.7083692044608846, -0.3209642339711405, 0.4525, "3d",
         "geocentric_bilinear"},
    };

    for (const auto &testPoint : testPoints) {
        j["components"][0]["displacement_type"] = testPoint.displacement_type;
        j["components"][0]["spatial_model"]["interpolation_method"] =
            testPoint.interpolation_method;
        Evaluator<Grid, GridSet, EvaluatorIface> eval(
            MasterFile::parse(j.dump()), iface, a, b);

        const double longitude = testPoint.longitude;
        const double lat = testPoint.lat;
        double newLong;
        double newLat;
        double newZ;
        EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                 zVal, tValid, newLong, newLat, newZ))
            << longitude << " " << lat << " " << testPoint.displacement_type
            << testPoint.interpolation_method;
        EXPECT_NEAR(newZ - zVal, tFactor * testPoint.expected_dz, 1e-8)
            << longitude << " " << lat << " " << testPoint.displacement_type
            << testPoint.interpolation_method;

        double de;
        double dn;
        DeltaLongLatToEastingNorthing(DegToRad(lat),
                                      newLong - DegToRad(longitude),
                                      newLat - DegToRad(lat), a, b, de, dn);
        EXPECT_NEAR(de, tFactor * testPoint.expected_de, 1e-8)
            << longitude << " " << lat << " " << testPoint.displacement_type
            << testPoint.interpolation_method;
        EXPECT_NEAR(dn, tFactor * testPoint.expected_dn, 1e-8)
            << longitude << " " << lat << " " << testPoint.displacement_type
            << testPoint.interpolation_method;

        if (longitude == gridMinX && lat == gridMinY) {
            // Redo the exact same test, to test caching
            double newLong2;
            double newLat2;
            double newZ2;
            EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                     zVal, tValid, newLong2, newLat2, newZ2));
            EXPECT_EQ(newLong2, newLong);
            EXPECT_EQ(newLat2, newLat);
            EXPECT_EQ(newZ2, newZ);

            // Shift in longitude
            EXPECT_TRUE(eval.forward(iface, DegToRad(longitude - gridResX / 2),
                                     DegToRad(lat), zVal, tValid, newLong2,
                                     newLat2, newZ2));

            // Redo test at original position
            EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                     zVal, tValid, newLong2, newLat2, newZ2));
            EXPECT_EQ(newLong2, newLong);
            EXPECT_EQ(newLat2, newLat);
            EXPECT_EQ(newZ2, newZ);

            // Shift in latitude
            EXPECT_TRUE(eval.forward(iface, DegToRad(longitude),
                                     DegToRad(lat - gridResY / 2), zVal, tValid,
                                     newLong2, newLat2, newZ2));

            // Redo test at original position
            EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                     zVal, tValid, newLong2, newLat2, newZ2));
            EXPECT_EQ(newLong2, newLong);
            EXPECT_EQ(newLat2, newLat);
            EXPECT_EQ(newZ2, newZ);
        }
    }

    // Test inverse()
    {
        j["horizontal_offset_method"] = "addition";
        j["components"][0]["displacement_type"] = "3d";
        j["components"][0]["spatial_model"]["interpolation_method"] =
            "bilinear";
        Evaluator<Grid, GridSet, EvaluatorIface> eval(
            MasterFile::parse(j.dump()), iface, a, b);

        const double longitude = 165.9;
        const double lat = -37.3;
        double newLong;
        double newLat;
        double newZ;
        EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                 zVal, tValid, newLong, newLat, newZ));

        double invLongitude;
        double invLat;
        double invZ;
        EXPECT_TRUE(eval.inverse(iface, newLong, newLat, newZ, tValid,
                                 invLongitude, invLat, invZ));
        EXPECT_NEAR(RadToDeg(invLongitude), longitude, 1e-10);
        EXPECT_NEAR(RadToDeg(invLat), lat, 1e-10);
        EXPECT_NEAR(invZ, zVal, 1e-4);
    }

    // Test horizontal_offset_method = geocentric
    {
        j["horizontal_offset_method"] = "geocentric";
        j["components"][0]["displacement_type"] = "3d";
        j["components"][0]["spatial_model"]["interpolation_method"] =
            "bilinear";
        Evaluator<Grid, GridSet, EvaluatorIface> eval(
            MasterFile::parse(j.dump()), iface, a, b);

        const double longitude = gridMinX;
        const double lat = gridMinY;
        double newLong;
        double newLat;
        double newZ;
        EXPECT_TRUE(eval.forward(iface, DegToRad(longitude), DegToRad(lat),
                                 zVal, tValid, newLong, newLat, newZ));

        double de;
        double dn;
        constexpr double expected_de = 0.40000000948081327;
        constexpr double expected_dn = -0.19999999810542682;
        DeltaLongLatToEastingNorthing(DegToRad(lat),
                                      newLong - DegToRad(longitude),
                                      newLat - DegToRad(lat), a, b, de, dn);
        EXPECT_NEAR(de, tFactor * expected_de, 1e-10);
        EXPECT_NEAR(dn, tFactor * expected_dn, 1e-9);
        EXPECT_NEAR(newZ - zVal, tFactor * 0.84, 1e-4);
    }
}

// ---------------------------------------------------------------------------

TEST(defmodel, evaluator_projected_crs) {

    json j(getMinValidContent());
    j["horizontal_offset_method"] = "addition";
    j["horizontal_offset_unit"] = "metre";
    j["vertical_offset_unit"] = "metre";
    constexpr double gridMinX = 10000;
    constexpr double gridMinY = 20000;
    constexpr double gridMaxX = 30000;
    constexpr double gridMaxY = 40000;
    constexpr double gridResX = gridMaxX - gridMinX;
    constexpr double gridResY = gridMaxY - gridMinY;

    j["extent"]["parameters"] = {
        {"bbox", {gridMinX, gridMinY, gridMaxX, gridMaxY}}};
    j["components"] = {
        {{"displacement_type", "horizontal"},
         {"uncertainty_type", "none"},
         {"extent",
          {{"type", "bbox"},
           {"parameters",
            {{"bbox", {gridMinX, gridMinY, gridMaxX, gridMaxY}}}}}},
         {"spatial_model",
          {
              {"type", "GeoTIFF"},
              {"interpolation_method", "bilinear"},
              {"filename", "bla.tif"},
          }},
         {"time_function", {{"type", "constant"}}}}};

    struct Grid : public GridPrototype {
        bool getEastingNorthingOffset(int ix, int iy, double &eastingOffset,
                                      double &northingOffset) const {
            if (ix == 0 && iy == 0) {
                eastingOffset = 0.4;
                northingOffset = -0.2;
            } else if (ix == 1 && iy == 0) {
                eastingOffset = 0.5;
                northingOffset = -0.25;
            } else if (ix == 0 && iy == 1) {
                eastingOffset = 0.8;
                northingOffset = -0.4;
            } else if (ix == 1 && iy == 1) {
                eastingOffset = 1.;
                northingOffset = -0.3;
            } else {
                return false;
            }
            return true;
        }

#ifdef DEBUG_DEFMODEL
        std::string name() const { return std::string(); }
#endif
    };

    struct GridSet : public GridSetPrototype<Grid> {

        Grid grid{};

        GridSet() {
            grid.minx = gridMinX;
            grid.miny = gridMinY;
            grid.resx = gridResX;
            grid.resy = gridResY;
            grid.width = 2;
            grid.height = 2;
        }

        const Grid *gridAt(double /*x */, double /* y */) { return &grid; }
    };

    struct EvaluatorIface : public EvaluatorIfacePrototype<Grid, GridSet> {
        std::unique_ptr<GridSet> open(const std::string &filename) {
            if (filename != "bla.tif")
                return nullptr;
            return std::unique_ptr<GridSet>(new GridSet());
        }

        bool isGeographicCRS(const std::string & /* crsDef */) { return false; }

#ifdef DEBUG_DEFMODEL
        void log(const std::string & /* msg */) {}
#endif
    };
    EvaluatorIface iface;

    constexpr double a = 6378137;
    constexpr double b = 6356752.314140;
    constexpr double tValid = 2018;
    constexpr double zVal = 100;

    Evaluator<Grid, GridSet, EvaluatorIface> eval(MasterFile::parse(j.dump()),
                                                  iface, a, b);

    double newX;
    double newY;
    double newZ;
    EXPECT_TRUE(eval.forward(iface, gridMinX, gridMinY, zVal, tValid, newX,
                             newY, newZ));
    EXPECT_NEAR(newX - gridMinX, 0.4, 1e-8);
    EXPECT_NEAR(newY - gridMinY, -0.2, 1e-8);
    EXPECT_NEAR(newZ - zVal, 0, 1e-8);

    {
        json jcopy(j);
        jcopy["horizontal_offset_unit"] = "degree";
        EXPECT_THROW((Evaluator<Grid, GridSet, EvaluatorIface>(
                         MasterFile::parse(jcopy.dump()), iface, a, b)),
                     EvaluatorException);
    }

    {
        json jcopy(j);
        jcopy["horizontal_offset_method"] = "geocentric";
        EXPECT_THROW((Evaluator<Grid, GridSet, EvaluatorIface>(
                         MasterFile::parse(jcopy.dump()), iface, a, b)),
                     EvaluatorException);
    }

    {
        json jcopy(j);
        jcopy["components"][0]["spatial_model"]["interpolation_method"] =
            "geocentric_bilinear";
        EXPECT_THROW((Evaluator<Grid, GridSet, EvaluatorIface>(
                         MasterFile::parse(jcopy.dump()), iface, a, b)),
                     EvaluatorException);
    }
}

} // namespace

#ifdef _MSC_VER
#pragma warning(pop)
#endif
