Duff.h

/* Copyright (c) Queen's Printer for Ontario, 2020. */
/* Copyright (c) His Majesty the King in Right of Canada as represented by the Minister of Natural Resources, 2025. */
 
/* SPDX-License-Identifier: AGPL-3.0-or-later */
 
#pragma once
#include "LookupTable.h"
namespace fs::fuel
{
template <int Ash, int Rho, int B0, int B1, int B2, int B3>
class DuffType;
class Duff
{
public:
  virtual ~Duff() = default;
  Duff(const Duff& rhs) noexcept = delete;
  Duff(Duff&& rhs) noexcept = delete;
  Duff& operator=(const Duff& rhs) noexcept = delete;
  Duff& operator=(Duff&& rhs) noexcept = delete;
  [[nodiscard]] constexpr bool operator==(const Duff& rhs) const
  {
    // HACK: only equivalent if identical
    return this == &rhs;
  }
  [[nodiscard]] constexpr bool operator!=(const Duff& rhs) const
  {
    return !operator==(rhs);
  }
  [[nodiscard]] virtual ThresholdSize probabilityOfSurvival(MathSize mc_pct) const noexcept = 0;
  //  /**
  //   * \brief Feather moss (upper) [Frandsen table 2/3]
  //   */
  //  static const DuffType<172, 464, 139873, -3296, 4904, 568> FeatherMossUpper;
  //  /**
  //   * \brief Feather moss (lower) [Frandsen table 2/3]
  //   */
  //  static const DuffType<191, 389, 132628, -1167, 3308, -2604> FeatherMossLower;
  static const DuffType<124, 218, -88306, -608, 8095, 2735> SphagnumUpper;
  //  /**
  //   * \brief Sphagnum (lower) [Frandsen table 2/3]
  //   */
  //  static const DuffType<567, 1190, 3273347, -37655, -87849, 26684> SphagnumLower;
  static const DuffType<181, 427, 90970, -1040, 1165, -646> FeatherMoss;
  static const DuffType<261, 563, 80359, -393, -591, -340> Reindeer;
  //  /**
  //   * \brief Sedge meadow (upper) [Frandsen table 2/3]
  //   */
  //  static const DuffType<233, 694, 398477, -1800, -3727, -1874> SedgeMeadowUpper;
  //  /**
  //   * \brief Sedge meadow (lower) [Frandsen table 2/3]
  //   */
  //  static const DuffType<449, 915, 290818, -2059, -2319, -420> SedgeMeadowLower;
  static const DuffType<359, 1220, 3325604, -12220, -21024, -12619> WhiteSpruce;
  static const DuffType<94, 2220, -198198, -1169, 10414, 782> Peat;
  static const DuffType<349, 2030, 372276, -1876, -2833, -951> PeatMuck;
  //  /**
  //   * \brief Sedge meadow (Seney) [Frandsen table 2/3]
  //   */
  //  static const DuffType<354, 1830, 71813, -1413, -1253, 390> SedgeMeadowSeney;
  static const DuffType<365, 1900, 451778, -3227, -3644, -362> PineSeney;
  static const DuffType<307, 1160, 586921, -2737, -5413, -1246> SprucePine;
  //  /**
  //   * \brief Grass/sedge marsh [Frandsen table 2/3]
  //   */
  //  static const DuffType<352, 1200, 2362934, -8423, -25097, -4902> GrassSedgeMarsh;
  //  /**
  //   * \brief Southern pine duff [Frandsen table 2/3]
  //   */
  //  static const DuffType<680, 1120, 586921, -2737, -5413, -1246> SouthernPine;
  //  /**
  //   * \brief Hardwood swamp (upper) [Frandsen table 2/3]
  //   */
  //  static const DuffType<182, 1380, 336907, -2946, -3002, -4040> HardwoodSwamp;
  // coefficients aren't defined in the table these came from
  // static const DuffType Pocosin;
  // static const DuffType SwampForest;
  // static const DuffType Flatwoods;
protected:
  Duff() = default;
};
template <int Ash, int Rho, int B0, int B1, int B2, int B3>
class DuffType final
  : public Duff
{
public:
  DuffType() = default;
  ~DuffType() override = default;
  DuffType(const DuffType& rhs) noexcept = delete;
  DuffType(DuffType&& rhs) noexcept = delete;
  DuffType& operator=(const DuffType& rhs) noexcept = delete;
  DuffType& operator=(DuffType&& rhs) noexcept = delete;
  [[nodiscard]] ThresholdSize probabilityOfSurvival(const MathSize mc_pct) const noexcept override
  {
    return probability_of_survival_(mc_pct);
  }
  [[nodiscard]] static constexpr MathSize ash()
  {
    return Ash / 10.0;
  }
  [[nodiscard]] static constexpr MathSize rho()
  {
    return Rho / 10.0;
  }
  [[nodiscard]] static constexpr MathSize b0()
  {
    return B0 / 10000.0;
  }
  [[nodiscard]] static constexpr MathSize b1()
  {
    return B1 / 10000.0;
  }
  [[nodiscard]] static constexpr MathSize b2()
  {
    return B2 / 10000.0;
  }
  [[nodiscard]] static constexpr MathSize b3()
  {
    return B3 / 10000.0;
  }
private:
  static constexpr auto ConstantPart = b0() + b2() * ash() + b3() * rho();
  [[nodiscard]] static ThresholdSize
    duffFunction(const MathSize mc_pct) noexcept
  {
    const auto d = 1 + exp(-(b1() * mc_pct + ConstantPart));
    if (0 == d)
    {
      return 1.0;
    }
    return 1.0 / d;
  }
  util::LookupTable<&duffFunction> probability_of_survival_{};
};
}