Compile Time International Standard Atmosphere

isa.cpp

#include <array>
#include <limits>
#include <type_traits>

#include <cmath>

namespace isa {

template<typename Params, int Offset=0>
class ISA : public Params {
    static constexpr std::size_t Layers = std::size(Params::Tgrads);
    using real_t = typename decltype(Params::layers)::value_type;

    static_assert(std::is_same_v<decltype(Params::layers), std::array<real_t, Layers + 1> const>);
    static_assert(std::is_same_v<decltype(Params::Tgrads), std::array<real_t, Layers> const>);
    static_assert(std::is_same_v<decltype(Params::R), real_t const>);
    static_assert(std::is_same_v<decltype(Params::g), real_t const>);
    static_assert(std::is_same_v<decltype(Params::T0), real_t const>);
    static_assert(std::is_same_v<decltype(Params::p0), real_t const>);
    static_assert(std::is_floating_point_v<real_t>);

    class Layer {
        std::size_t l_;
        real_t T0_;
        real_t p0_;

    public:
        constexpr Layer(std::size_t l, real_t T0, real_t p0) noexcept
            : l_{l}
            , T0_{T0}
            , p0_{p0}
        {}

        [[nodiscard]] inline constexpr real_t T(real_t h) const noexcept {
            return T0_ + Params::Tgrads[l_] * (h - Params::layers[l_]);
        }

        [[nodiscard]] inline constexpr real_t p(real_t h) const noexcept {
            if (std::fabs(Params::Tgrads[l_]) < std::numeric_limits<real_t>::epsilon())
                return p0_ * std::exp(
                    -Params::g * (h - Params::layers[l_]) /
                    (Params::R * T0_)
                );
            else
                return p0_ * std::pow(
                    (real_t(1.) + Params::Tgrads[l_] * (h - Params::layers[l_]) / T0_),
                    -Params::g / (Params::R * Params::Tgrads[l_])
                );
        }
    };

    template<std::size_t... I>
    [[nodiscard]] static constexpr auto make_layers_(std::index_sequence<I...>) noexcept {
        return std::array<Layer, Layers>{make_layer_<I>({0, Params::T0 + Offset, Params::p0})...};
    }

    template<std::size_t N, std::size_t I=0>
    [[nodiscard]] static constexpr auto make_layer_(Layer layer) noexcept {
        if constexpr (I < N)
            return make_layer_<N, I + 1>({
                I + 1,
                layer.T(Params::layers[I + 1]),
                layer.p(Params::layers[I + 1])
            });
        return layer;
    }

    template<std::size_t I=0>
    [[nodiscard]] static constexpr auto const& get_layer_([[maybe_unused]] real_t h) noexcept {
        if constexpr (I < Layers - 1) {
            if (h < Params::layers[I + 1])
                return layers_[I];
            return get_layer_<I + 1>(h);
        }
        return layers_[I];
    }

    static constexpr auto layers_{make_layers_(std::make_index_sequence<Layers>())};

public:
    [[nodiscard]] static constexpr real_t T(real_t h) noexcept {
        if (h < Params::layers[0] || h > Params::layers[Layers])
            return std::numeric_limits<real_t>::quiet_NaN();
        auto const& layer = get_layer_(h);
        return layer.T(h);
    }

    [[nodiscard]] static constexpr real_t p(real_t h) noexcept {
        if (h < Params::layers[0] || h > Params::layers[Layers])
            return std::numeric_limits<real_t>::quiet_NaN();
        auto const& layer = get_layer_(h);
        return layer.p(h);
    }

    [[nodiscard]] static constexpr real_t rho(real_t h) noexcept {
        if (h < Params::layers[0] || h > Params::layers[Layers])
            return std::numeric_limits<real_t>::quiet_NaN();
        auto const& layer = get_layer_(h);
        return layer.p(h) / (Params::R * layer.T(h));
    }

    template<int O>
    using offset = ISA<Params, O>;
};


struct EarthParams {
    static constexpr auto layers = std::array{
        00000.0,  // ground
        11000.0,  // troposphere limit
        20000.0,  // tropopause limit
        32000.0,  
        47000.0,  // stratosphere limit
        51000.0,  // stratopause limit
        71000.0,
        84552.0,  // mesosphere limit
        86000.0   // mesopause limit
    };
    static constexpr auto Tgrads = std::array{
        -0.0065,  // troposphere
        +0000.0,  // tropopause
        +0.0010,  // stratosphere
        +0.0028,  // stratosphere
        +0000.0,  // stratopause
        -0.0028,  // mesosphere
        -0.0020,  // mesosphere
        +0000.0   // mesopause
    };

    static constexpr auto R = 287.05287;
    static constexpr auto g = 9.80665;
    static constexpr auto T0 = 288.15;    // 15 ºC
    static constexpr auto p0 = 101325.0;  // 1 atm
};
using EarthISA = ISA<EarthParams>;

}