helios_vector_types.h

Go to the documentation of this file.

 
#ifndef HELIOS_VECTOR_TYPES
#define HELIOS_VECTOR_TYPES
 
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
 
typedef unsigned int uint;
 
#include <cassert>
#include <string>
#include <vector>
// #include <stdexcept>
#include <cmath>
#include <iomanip>
#include <limits>
 
 
namespace helios {
 
 
    struct int2 {
 
        int x;
        int y;
 
        /* initializes to zero */
        constexpr int2() : x(0), y(0) {
        }
 
        explicit int2(const std::vector<int> &v) {
            if (v.size() != 2) {
                throw(std::runtime_error("ERROR: vector for initialization of int2 must have length 2."));
            }
            x = v[0];
            y = v[1];
        }
        constexpr explicit int2(const int v[2]) : x(v[0]), y(v[1]) {
        }
 
        constexpr int2(int v0, int v1) : x(v0), y(v1) {
        }
 
 
        constexpr int2 operator+(const int2 &a) const noexcept;
        inline int2 &operator+=(const int2 &a) noexcept;
        inline int2 &operator-=(const int2 &a) noexcept;
        constexpr int2 operator-(const int2 &a) const noexcept;
        constexpr bool operator==(const int2 &a) const noexcept;
        constexpr bool operator!=(const int2 &a) const noexcept;
        constexpr int2 operator-() const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, const helios::int2 &vec) {
            return os << "helios::int2<" << vec.x << ", " << vec.y << ">";
        }
    };
 
 
    constexpr int2 make_int2(int x, int y) {
        return {x, y};
    }
 
 
    inline int2 make_int2(int X[2]) {
        return {X[0], X[1]};
    }
 
    constexpr int2 int2::operator+(const int2 &a) const noexcept {
        return {a.x + x, a.y + y};
    }
 
    inline int2 &int2::operator+=(const int2 &a) noexcept {
        this->x += a.x;
        this->y += a.y;
        return *this;
    }
 
    inline int2 &int2::operator-=(const int2 &a) noexcept {
        this->x -= a.x;
        this->y -= a.y;
        return *this;
    }
 
    constexpr int2 int2::operator-(const int2 &a) const noexcept {
        return {x - a.x, y - a.y};
    }
 
    constexpr bool int2::operator==(const int2 &a) const noexcept {
        return a.x == x && a.y == y;
    }
 
    constexpr bool int2::operator!=(const int2 &a) const noexcept {
        return a.x != x || a.y != y;
    }
 
    constexpr int2 int2::operator-() const noexcept {
        return {-x, -y};
    }
 
 
    struct int3 {
 
        int x;
        int y;
        int z;
 
        /* initializes to zero */
        constexpr int3() : x(0), y(0), z(0) {
        }
 
        explicit int3(const std::vector<int> &v) {
            if (v.size() != 3) {
                throw(std::runtime_error("ERROR: vector for initialization of int3 must have length 3."));
            }
            x = v[0];
            y = v[1];
            z = v[2];
        }
        constexpr explicit int3(const int v[3]) : x(v[0]), y(v[1]), z(v[2]) {
        }
 
        constexpr int3(int v0, int v1, int v2) : x(v0), y(v1), z(v2) {
        }
 
        constexpr int3 operator+(const int3 &a) const noexcept;
        inline int3 &operator+=(const int3 &a) noexcept;
        inline int3 &operator-=(const int3 &a) noexcept;
        constexpr int3 operator-(const int3 &a) const noexcept;
        constexpr bool operator==(const int3 &a) const noexcept;
        constexpr bool operator!=(const int3 &a) const noexcept;
        constexpr int3 operator-() const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, const helios::int3 &vec) {
            return os << "helios::int3<" << vec.x << ", " << vec.y << ", " << vec.z << ">";
        }
    };
 
    constexpr int3 make_int3(int X, int Y, int Z) {
        return {X, Y, Z};
    }
 
    constexpr int3 make_int3(int X[3]) {
        return {X[0], X[1], X[2]};
    }
 
    constexpr int3 int3::operator+(const int3 &a) const noexcept {
        return {a.x + x, a.y + y, a.z + z};
    }
 
    inline int3 &int3::operator+=(const int3 &a) noexcept {
        this->x += a.x;
        this->y += a.y;
        this->z += a.z;
        return *this;
    }
 
    inline int3 &int3::operator-=(const int3 &a) noexcept {
        this->x -= a.x;
        this->y -= a.y;
        this->z -= a.z;
        return *this;
    }
 
    constexpr int3 int3::operator-(const int3 &a) const noexcept {
        return {x - a.x, y - a.y, z - a.z};
    }
 
    constexpr bool int3::operator==(const int3 &a) const noexcept {
        return a.x == x && a.y == y && a.z == z;
    }
 
    constexpr bool int3::operator!=(const int3 &a) const noexcept {
        return a.x != x || a.y != y || a.z != z;
    }
 
    constexpr int3 int3::operator-() const noexcept {
        return {-x, -y, -z};
    }
 
 
    struct int4 {
 
 
        int x;
        int y;
        int z;
        int w;
 
        /* initializes to zero */
        constexpr int4() : x(0), y(0), z(0), w(0) {
        }
 
        explicit int4(const std::vector<int> &v) {
            if (v.size() != 4) {
                throw(std::runtime_error("ERROR: vector for initialization of int4 must have length 4."));
            }
            x = v[0];
            y = v[1];
            z = v[2];
            w = v[3];
        }
 
        explicit constexpr int4(const int v[4]) : x(v[0]), y(v[1]), z(v[2]), w(v[3]) {
        }
 
        constexpr int4(int v0, int v1, int v2, int v3) : x(v0), y(v1), z(v2), w(v3) {
        }
 
        constexpr int4 operator+(const int4 &a) const noexcept;
        inline int4 &operator+=(const int4 &a) noexcept;
        inline int4 &operator-=(const int4 &a) noexcept;
        constexpr int4 operator-(const int4 &a) const noexcept;
        constexpr bool operator==(const int4 &a) const noexcept;
        constexpr bool operator!=(const int4 &a) const noexcept;
        constexpr int4 operator-() const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, const helios::int4 &vec) {
            return os << "helios::int4<" << vec.x << ", " << vec.y << ", " << vec.z << ", " << vec.w << ">";
        }
    };
 
    constexpr int4 make_int4(int x, int y, int z, int w) {
        return {x, y, z, w};
    }
 
    constexpr int4 make_int4(const int X[4]) {
        return {X[0], X[1], X[2], X[3]};
    }
 
    constexpr int4 int4::operator+(const int4 &a) const noexcept {
        return {a.x + x, a.y + y, a.z + z, a.w + w};
    }
 
    inline int4 &int4::operator+=(const int4 &a) noexcept {
        this->x += a.x;
        this->y += a.y;
        this->z += a.z;
        this->w += a.w;
        return *this;
    }
 
    inline int4 &int4::operator-=(const int4 &a) noexcept {
        this->x -= a.x;
        this->y -= a.y;
        this->z -= a.z;
        this->w -= a.w;
        return *this;
    }
 
    constexpr int4 int4::operator-(const int4 &a) const noexcept {
        return {x - a.x, y - a.y, z - a.z, w - a.w};
    }
 
    constexpr bool int4::operator==(const int4 &a) const noexcept {
        return a.x == x && a.y == y && a.z == z && a.w == w;
    }
 
    constexpr bool int4::operator!=(const int4 &a) const noexcept {
        return a.x != x || a.y != y || a.z != z || a.w != w;
    }
 
    constexpr int4 int4::operator-() const noexcept {
        return {-x, -y, -z, -w};
    }
 
 
    struct uint2 {
 
 
        unsigned int x;
        unsigned int y;
 
        /* initializes to zero */
        constexpr uint2() : x(0u), y(0u) {
        }
 
        explicit uint2(const std::vector<unsigned int> &v) {
            if (v.size() != 2) {
                throw(std::runtime_error("ERROR: vector for initialization of uint2 must have length 2."));
            }
            x = v[0];
            y = v[1];
        }
        constexpr explicit uint2(const unsigned int v[2]) : x(v[0]), y(v[1]) {
        }
 
        constexpr uint2(unsigned int v0, unsigned int v1) : x(v0), y(v1) {
        }
 
 
        constexpr uint2 operator+(const uint2 &a) const noexcept;
        inline uint2 &operator+=(const uint2 &a) noexcept;
        inline uint2 &operator-=(const uint2 &a) noexcept;
        constexpr uint2 operator-(const uint2 &a) const noexcept;
        constexpr bool operator==(const uint2 &a) const noexcept;
        constexpr bool operator!=(const uint2 &a) const noexcept;
        constexpr uint2 operator-() const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, const helios::uint2 &vec) {
            return os << "helios::uint2<" << vec.x << ", " << vec.y << ">";
        }
    };
 
 
    constexpr uint2 make_uint2(unsigned int x, unsigned int y) {
        return {x, y};
    }
 
 
    constexpr uint2 make_uint2(unsigned int X[2]) {
        return {X[0], X[1]};
    }
 
    constexpr uint2 uint2::operator+(const uint2 &a) const noexcept {
        return {a.x + x, a.y + y};
    }
 
    inline uint2 &uint2::operator+=(const uint2 &a) noexcept {
        this->x += a.x;
        this->y += a.y;
        return *this;
    }
 
    inline uint2 &uint2::operator-=(const uint2 &a) noexcept {
        this->x -= a.x;
        this->y -= a.y;
        return *this;
    }
 
    constexpr uint2 uint2::operator-(const uint2 &a) const noexcept {
        return {x - a.x, y - a.y};
    }
 
    constexpr bool uint2::operator==(const uint2 &a) const noexcept {
        return a.x == x && a.y == y;
    }
 
    constexpr bool uint2::operator!=(const uint2 &a) const noexcept {
        return a.x != x || a.y != y;
    }
 
    constexpr uint2 uint2::operator-() const noexcept {
        return {-x, -y};
    }
 
 
 
    struct uint3 {
 
 
        unsigned int x;
        unsigned int y;
        unsigned int z;
 
 
        constexpr uint3() : x(0u), y(0u), z(0u) {
        }
 
        explicit uint3(const std::vector<unsigned int> &v) {
            if (v.size() != 3) {
                throw(std::runtime_error("ERROR: vector for initialization of uint3 must have length 3."));
            }
            x = v[0];
            y = v[1];
            z = v[2];
        }
        constexpr explicit uint3(const unsigned int v[3]) : x(v[0]), y(v[1]), z(v[2]) {
        }
 
        constexpr uint3(unsigned int v0, unsigned int v1, unsigned int v2) : x(v0), y(v1), z(v2) {
        }
 
        constexpr uint3 operator+(const uint3 &a) const noexcept;
        inline uint3 &operator+=(const uint3 &a) noexcept;
        inline uint3 &operator-=(const uint3 &a) noexcept;
        constexpr uint3 operator-(const uint3 &a) const noexcept;
        constexpr bool operator==(const uint3 &a) const noexcept;
        constexpr bool operator!=(const uint3 &a) const noexcept;
        constexpr uint3 operator-() const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, const helios::uint3 &vec) {
            return os << "helios::uint3<" << vec.x << ", " << vec.y << ", " << vec.z << ">";
        }
    };
 
 
    constexpr uint3 make_uint3(unsigned int X, unsigned int Y, unsigned int Z) {
        return {X, Y, Z};
    }
 
 
    constexpr uint3 make_uint3(unsigned int X[3]) {
        return {X[0], X[1], X[2]};
    }
 
    constexpr uint3 uint3::operator+(const uint3 &a) const noexcept {
        return {a.x + x, a.y + y, a.z + z};
    }
 
    inline uint3 &uint3::operator+=(const uint3 &a) noexcept {
        this->x += a.x;
        this->y += a.y;
        this->z += a.z;
        return *this;
    }
 
    inline uint3 &uint3::operator-=(const uint3 &a) noexcept {
        this->x -= a.x;
        this->y -= a.y;
        this->z -= a.z;
        return *this;
    }
 
    constexpr uint3 uint3::operator-(const uint3 &a) const noexcept {
        return {x - a.x, y - a.y, z - a.z};
    }
 
    constexpr bool uint3::operator==(const uint3 &a) const noexcept {
        return a.x == x && a.y == y && a.z == z;
    }
 
    constexpr bool uint3::operator!=(const uint3 &a) const noexcept {
        return a.x != x || a.y != y || a.z != z;
    }
 
    constexpr uint3 uint3::operator-() const noexcept {
        return {-x, -y, -z};
    }
 
 
 
    struct uint4 {
 
        unsigned int x;
        unsigned int y;
        unsigned int z;
        unsigned int w;
 
 
        constexpr uint4() : x(0u), y(0u), z(0u), w(0u) {
        }
 
        explicit uint4(const std::vector<unsigned int> &v) {
            if (v.size() != 4) {
                throw(std::runtime_error("ERROR: vector for initialization of uint4 must have length 4."));
            }
            x = v[0];
            y = v[1];
            z = v[2];
            w = v[3];
        }
 
        constexpr explicit uint4(const unsigned int v[4]) : x(v[0]), y(v[1]), z(v[2]), w(v[3]) {
        }
 
        constexpr uint4(unsigned int v0, unsigned int v1, unsigned int v2, unsigned int v3) : x(v0), y(v1), z(v2), w(v3) {
        }
 
        constexpr uint4 operator+(const uint4 &a) const noexcept;
        inline uint4 &operator+=(const uint4 &a) noexcept;
        inline uint4 &operator-=(const uint4 &a) noexcept;
        constexpr uint4 operator-(const uint4 &a) const noexcept;
        constexpr bool operator==(const uint4 &a) const noexcept;
        constexpr bool operator!=(const uint4 &a) const noexcept;
        constexpr uint4 operator-() const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, const helios::uint4 &vec) {
            return os << "helios::uint4<" << vec.x << ", " << vec.y << ", " << vec.z << ", " << vec.w << ">";
        }
    };
 
 
    constexpr uint4 make_uint4(unsigned int x, unsigned int y, unsigned int z, unsigned int w) {
        return {x, y, z, w};
    }
 
 
    constexpr uint4 make_uint4(const unsigned int X[4]) {
        return {X[0], X[1], X[2], X[3]};
    }
 
    constexpr uint4 uint4::operator+(const uint4 &a) const noexcept {
        return {a.x + x, a.y + y, a.z + z, a.w + w};
    }
 
    inline uint4 &uint4::operator+=(const uint4 &a) noexcept {
        this->x += a.x;
        this->y += a.y;
        this->z += a.z;
        this->w += a.w;
        return *this;
    }
 
    inline uint4 &uint4::operator-=(const uint4 &a) noexcept {
        this->x -= a.x;
        this->y -= a.y;
        this->z -= a.z;
        this->w -= a.w;
        return *this;
    }
 
    constexpr uint4 uint4::operator-(const uint4 &a) const noexcept {
        return {x - a.x, y - a.y, z - a.z, w - a.w};
    }
 
    constexpr bool uint4::operator==(const uint4 &a) const noexcept {
        return a.x == x && a.y == y && a.z == z && a.w == w;
    }
 
    constexpr bool uint4::operator!=(const uint4 &a) const noexcept {
        return a.x != x || a.y != y || a.z != z || a.w != w;
    }
 
    constexpr uint4 uint4::operator-() const noexcept {
        return {static_cast<unsigned int>(0u - x), static_cast<unsigned int>(0u - y), static_cast<unsigned int>(0u - z), static_cast<unsigned int>(0u - w)};
    }
 
 
    struct vec2 {
 
 
        float x;
        float y;
 
        vec2 normalize() {
            const float m2 = x * x + y * y;
            if (m2 <= 0)
                return *this;
            const float inv = 1.0f / sqrtf(m2);
            this->x *= inv;
            this->y *= inv;
            return *this;
        }
 
 
        [[nodiscard]] float magnitude() const {
            return sqrtf(x * x + y * y);
        }
 
 
        constexpr vec2() : x(0), y(0) {
        }
 
        explicit vec2(const std::vector<float> &v) {
            if (v.size() != 2) {
                throw(std::runtime_error("ERROR: vector for initialization of vec2 must have length 2."));
            }
            x = v[0];
            y = v[1];
        }
        explicit constexpr vec2(const float v[2]) : x(v[0]), y(v[1]) {
        }
 
        constexpr vec2(float v0, float v1) : x(v0), y(v1) {
        }
 
        constexpr float operator*(const vec2 &a) const noexcept;
        constexpr vec2 operator+(const vec2 &a) const noexcept;
        inline vec2 &operator+=(const vec2 &a) noexcept;
        inline vec2 &operator-=(const vec2 &a) noexcept;
        inline vec2 &operator*=(float a) noexcept;
        inline vec2 &operator/=(float a) noexcept;
        constexpr vec2 operator-(const vec2 &a) const noexcept;
        constexpr vec2 operator*(float a) const noexcept;
        inline vec2 operator/(float a) const noexcept;
        constexpr vec2 operator+(float a) const noexcept;
        constexpr vec2 operator-(float a) const noexcept;
        constexpr bool operator==(const vec2 &a) const noexcept;
        constexpr bool operator!=(const vec2 &a) const noexcept;
        constexpr vec2 operator-() const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, const helios::vec2 &vec) {
            return os << "helios::vec2<" << vec.x << ", " << vec.y << ">";
        }
    };
 
    constexpr vec2 operator*(float a, const vec2 &v) noexcept;
    constexpr vec2 operator+(float a, const vec2 &v) noexcept;
    constexpr vec2 operator-(float a, const vec2 &v) noexcept;
 
 
    constexpr vec2 make_vec2(float x, float y) {
        return {x, y};
    }
 
    inline vec2 normalize(const vec2 &v) {
        const float m2 = v.x * v.x + v.y * v.y;
        const float inv = 1.0f / sqrtf(m2);
        return {v.x * inv, v.y * inv};
    }
 
    constexpr float vec2::operator*(const vec2 &a) const noexcept {
        return a.x * x + a.y * y;
    }
 
    constexpr vec2 vec2::operator+(const vec2 &a) const noexcept {
        return {a.x + x, a.y + y};
    }
 
    inline vec2 &vec2::operator+=(const vec2 &a) noexcept {
        x += a.x;
        y += a.y;
        return *this;
    }
 
    inline vec2 &vec2::operator-=(const vec2 &a) noexcept {
        x -= a.x;
        y -= a.y;
        return *this;
    }
 
    inline vec2 &vec2::operator*=(float a) noexcept {
        x *= a;
        y *= a;
        return *this;
    }
 
    inline vec2 &vec2::operator/=(float a) noexcept {
        x /= a;
        y /= a;
        return *this;
    }
 
    constexpr vec2 vec2::operator+(const float a) const noexcept {
        return {a + x, a + y};
    }
 
    constexpr vec2 operator+(const float a, const vec2 &v) noexcept {
        return {a + v.x, a + v.y};
    }
 
    constexpr vec2 vec2::operator-(const vec2 &a) const noexcept {
        return {x - a.x, y - a.y};
    }
 
    constexpr vec2 vec2::operator-(const float a) const noexcept {
        return {x - a, y - a};
    }
 
    constexpr vec2 operator-(const float a, const vec2 &v) noexcept {
        return {a - v.x, a - v.y};
    }
 
    constexpr vec2 vec2::operator*(const float a) const noexcept {
        return {x * a, y * a};
    }
 
    constexpr vec2 operator*(const float a, const vec2 &v) noexcept {
        return {a * v.x, a * v.y};
    }
 
    inline vec2 vec2::operator/(const float a) const noexcept {
        const float inv_a = (a != 0.f ? 1.f / a : std::copysign(std::numeric_limits<float>::infinity(), a));
        return {x * inv_a, y * inv_a};
    }
 
    constexpr bool vec2::operator==(const vec2 &a) const noexcept {
        return x == a.x && y == a.y;
    }
 
    constexpr bool vec2::operator!=(const vec2 &a) const noexcept {
        return x != a.x || y != a.y;
    }
 
    constexpr vec2 vec2::operator-() const noexcept {
        return {-x, -y};
    }
 
 
    struct vec3 {
 
        float x;
        float y;
        float z;
 
        vec3 normalize() {
            const float m2 = x * x + y * y + z * z;
            if (m2 <= 0)
                return *this;
            const float inv = 1.0f / sqrtf(m2);
            this->x *= inv;
            this->y *= inv;
            this->z *= inv;
            return *this;
        }
 
 
        [[nodiscard]] float magnitude() const {
            return sqrtf(x * x + y * y + z * z);
        }
 
 
        constexpr vec3() : x(0), y(0), z(0) {
        }
 
        explicit vec3(const std::vector<float> &v) {
            if (v.size() != 3) {
                throw(std::runtime_error("ERROR: vector for initialization of vec3 must have length 3."));
            }
            x = v[0];
            y = v[1];
            z = v[2];
        }
        explicit constexpr vec3(const float v[3]) : x(v[0]), y(v[1]), z(v[2]) {
        }
 
        constexpr vec3(float v0, float v1, float v2) : x(v0), y(v1), z(v2) {
        }
 
        constexpr float operator*(const vec3 &a) const noexcept;
        constexpr vec3 operator+(const vec3 &a) const noexcept;
        inline vec3 &operator+=(const vec3 &a) noexcept;
        inline vec3 &operator-=(const vec3 &a) noexcept;
        inline vec3 &operator*=(float a) noexcept;
        inline vec3 &operator/=(float a) noexcept;
        constexpr vec3 operator-(const vec3 &a) const noexcept;
        constexpr vec3 operator*(float a) const noexcept;
        inline vec3 operator/(float a) const noexcept;
        constexpr vec3 operator+(float a) const noexcept;
        constexpr vec3 operator-(float a) const noexcept;
        constexpr bool operator==(const vec3 &a) const noexcept;
        constexpr bool operator!=(const vec3 &a) const noexcept;
        constexpr vec3 operator-() const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, const helios::vec3 &vec) {
            return os << "helios::vec3<" << vec.x << ", " << vec.y << ", " << vec.z << ">";
        }
    };
 
    constexpr vec3 operator*(float a, const vec3 &v) noexcept;
    constexpr vec3 operator+(float a, const vec3 &v) noexcept;
    constexpr vec3 operator-(float a, const vec3 &v) noexcept;
 
    constexpr vec3 make_vec3(float x, float y, float z) {
        return {x, y, z};
    }
 
    constexpr vec3 make_vec3(float X[3]) {
        return {X[0], X[1], X[2]};
    }
 
    constexpr vec3 cross(const vec3 &a, const vec3 &b) {
        return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
    }
 
    inline vec3 normalize(const vec3 &v) noexcept {
        const float m2 = v.x * v.x + v.y * v.y + v.z * v.z;
        const float inv = 1.0f / sqrtf(m2);
        return {v.x * inv, v.y * inv, v.z * inv};
    }
 
    constexpr float vec3::operator*(const vec3 &a) const noexcept {
        return a.x * x + a.y * y + a.z * z;
    }
 
    constexpr vec3 vec3::operator+(const vec3 &a) const noexcept {
        return {a.x + x, a.y + y, a.z + z};
    }
 
    inline vec3 &vec3::operator+=(const vec3 &a) noexcept {
        this->x += a.x;
        this->y += a.y;
        this->z += a.z;
        return *this;
    }
 
    inline vec3 &vec3::operator-=(const vec3 &a) noexcept {
        this->x -= a.x;
        this->y -= a.y;
        this->z -= a.z;
        return *this;
    }
 
    inline vec3 &vec3::operator*=(float a) noexcept {
        x *= a;
        y *= a;
        z *= a;
        return *this;
    }
 
    inline vec3 &vec3::operator/=(float a) noexcept {
        x /= a;
        y /= a;
        z /= a;
        return *this;
    }
 
    constexpr vec3 vec3::operator+(float a) const noexcept {
        return {x + a, y + a, z + a};
    }
 
    constexpr vec3 operator+(float a, const vec3 &v) noexcept {
        return {a + v.x, a + v.y, a + v.z};
    }
 
    constexpr vec3 vec3::operator-(const vec3 &a) const noexcept {
        return {x - a.x, y - a.y, z - a.z};
    }
 
    constexpr vec3 vec3::operator-(float a) const noexcept {
        return {x - a, y - a, z - a};
    }
 
    constexpr vec3 operator-(float a, const vec3 &v) noexcept {
        return {a - v.x, a - v.y, a - v.z};
    }
 
    constexpr vec3 vec3::operator*(float a) const noexcept {
        return {x * a, y * a, z * a};
    }
 
    constexpr vec3 operator*(float a, const vec3 &v) noexcept {
        return {a * v.x, a * v.y, a * v.z};
    }
 
    inline vec3 vec3::operator/(const float a) const noexcept {
        const float inv_a = (a != 0.f ? 1.f / a : std::copysign(std::numeric_limits<float>::infinity(), a));
        return {x * inv_a, y * inv_a, z * inv_a};
    }
 
    constexpr bool vec3::operator==(const vec3 &a) const noexcept {
        return x == a.x && y == a.y && z == a.z;
    }
 
    constexpr bool vec3::operator!=(const vec3 &a) const noexcept {
        return x != a.x || y != a.y || z != a.z;
    }
 
    constexpr vec3 vec3::operator-() const noexcept {
        return {-x, -y, -z};
    }
 
 
    struct vec4 {
 
 
        float x;
        float y;
        float z;
        float w;
 
        vec4 normalize() {
            const float m2 = x * x + y * y + z * z + w * w;
            if (m2 <= 0)
                return *this;
            const float inv = 1.0f / sqrtf(m2);
            this->x *= inv;
            this->y *= inv;
            this->z *= inv;
            this->w *= inv;
            return *this;
        }
 
 
        [[nodiscard]] float magnitude() const {
            return sqrt(x * x + y * y + z * z + w * w);
        }
 
 
        constexpr vec4() : x(0), y(0), z(0), w(0) {
        }
 
        explicit vec4(const std::vector<float> &v) {
            if (v.size() != 4) {
                throw(std::runtime_error("ERROR: vector for initialization of vec4 must have length 4."));
            }
            x = v[0];
            y = v[1];
            z = v[2];
            w = v[3];
        }
        explicit constexpr vec4(const float v[4]) : x(v[0]), y(v[1]), z(v[2]), w(v[3]) {
        }
 
        constexpr vec4(float v0, float v1, float v2, float v3) : x(v0), y(v1), z(v2), w(v3) {
        }
 
        constexpr float operator*(const vec4 &a) const noexcept;
        constexpr vec4 operator+(const vec4 &a) const noexcept;
        inline vec4 &operator+=(const vec4 &a) noexcept;
        inline vec4 &operator-=(const vec4 &a) noexcept;
        inline vec4 &operator*=(float a) noexcept;
        inline vec4 &operator/=(float a) noexcept;
        constexpr vec4 operator-(const vec4 &a) const noexcept;
        constexpr vec4 operator*(float a) const noexcept;
        inline vec4 operator/(float a) const noexcept;
        constexpr vec4 operator+(float a) const noexcept;
        constexpr vec4 operator-(float a) const noexcept;
        constexpr bool operator==(const vec4 &a) const noexcept;
        constexpr bool operator!=(const vec4 &a) const noexcept;
        constexpr vec4 operator-() const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, const helios::vec4 &vec) {
            return os << "helios::vec4<" << vec.x << ", " << vec.y << ", " << vec.z << ", " << vec.w << ">";
        }
    };
 
    constexpr vec4 operator*(float a, const vec4 &v) noexcept;
    constexpr vec4 operator+(float a, const vec4 &v) noexcept;
    constexpr vec4 operator-(float a, const vec4 &v) noexcept;
 
    constexpr vec4 make_vec4(float x, float y, float z, float w) {
        return {x, y, z, w};
    }
 
    constexpr vec4 make_vec4(float X[4]) {
        return {X[0], X[1], X[2], X[3]};
    }
 
    inline vec4 normalize(const vec4 &v) {
        const float m2 = v.x * v.x + v.y * v.y + v.z * v.z + v.w * v.w;
        const float inv = 1.0f / sqrtf(m2);
        return {v.x * inv, v.y * inv, v.z * inv, v.w * inv};
    }
 
    constexpr float vec4::operator*(const vec4 &a) const noexcept {
        return a.x * x + a.y * y + a.z * z + a.w * w;
    }
 
    constexpr vec4 vec4::operator+(const vec4 &a) const noexcept {
        return {a.x + x, a.y + y, a.z + z, a.w + w};
    }
 
    inline vec4 &vec4::operator+=(const vec4 &a) noexcept {
        this->x += a.x;
        this->y += a.y;
        this->z += a.z;
        this->w += a.w;
        return *this;
    }
 
    inline vec4 &vec4::operator-=(const vec4 &a) noexcept {
        this->x -= a.x;
        this->y -= a.y;
        this->z -= a.z;
        this->w -= a.w;
        return *this;
    }
 
    inline vec4 &vec4::operator*=(float a) noexcept {
        x *= a;
        y *= a;
        z *= a;
        w *= a;
        return *this;
    }
 
    inline vec4 &vec4::operator/=(float a) noexcept {
        x /= a;
        y /= a;
        z /= a;
        w /= a;
        return *this;
    }
 
    constexpr vec4 vec4::operator+(float a) const noexcept {
        return {x + a, y + a, z + a, w + a};
    }
 
    constexpr vec4 operator+(float a, const vec4 &v) noexcept {
        return {a + v.x, a + v.y, a + v.z, a + v.w};
    }
 
    constexpr vec4 vec4::operator-(const vec4 &a) const noexcept {
        return {x - a.x, y - a.y, z - a.z, w - a.w};
    }
 
    constexpr vec4 vec4::operator-(float a) const noexcept {
        return {x - a, y - a, z - a, w - a};
    }
 
    constexpr vec4 operator-(float a, const vec4 &v) noexcept {
        return {a - v.x, a - v.y, a - v.z, a - v.w};
    }
 
    constexpr vec4 vec4::operator*(float a) const noexcept {
        return {x * a, y * a, z * a, w * a};
    }
 
    constexpr vec4 operator*(float a, const vec4 &v) noexcept {
        return {a * v.x, a * v.y, a * v.z, a * v.w};
    }
 
    inline vec4 vec4::operator/(const float a) const noexcept {
        const float inv_a = (a != 0.f ? 1.f / a : std::copysign(std::numeric_limits<float>::infinity(), a));
        return {x * inv_a, y * inv_a, z * inv_a, w * inv_a};
    }
 
    constexpr bool vec4::operator==(const vec4 &a) const noexcept {
        return x == a.x && y == a.y && z == a.z && w == a.w;
    }
 
    constexpr bool vec4::operator!=(const vec4 &a) const noexcept {
        return x != a.x || y != a.y || z != a.z || w != a.w;
    }
 
    constexpr vec4 vec4::operator-() const noexcept {
        return {-x, -y, -z, -w};
    }
 
 
    struct RGBcolor {
 
 
 
        float r;
 
        float g;
 
        float b;
 
        constexpr RGBcolor() : r(0), g(0), b(0) {
        }
 
 
        constexpr RGBcolor(float red, float green, float blue) : r(RGBcolor::clamp(red)), g(RGBcolor::clamp(green)), b(RGBcolor::clamp(blue)) {
        }
 
 
        explicit constexpr RGBcolor(const float C[3]) : r(RGBcolor::clamp(C[0])), g(RGBcolor::clamp(C[1])), b(RGBcolor::clamp(C[2])) {
        }
 
 
        explicit RGBcolor(const std::vector<float> &C) {
            if (C.size() != 3) {
                throw(std::runtime_error("ERROR: cannot initialize RGBcolor using supplied vector - size should be 3."));
            }
            r = RGBcolor::clamp(C[0]);
            g = RGBcolor::clamp(C[1]);
            b = RGBcolor::clamp(C[2]);
        }
 
        explicit constexpr RGBcolor(const vec3 &C) : r(RGBcolor::clamp(C.x)), g(RGBcolor::clamp(C.y)), b(RGBcolor::clamp(C.z)) {
        }
 
 
        void scale(float scale_factor) {
            if (scale_factor < 0) {
                throw(std::runtime_error("ERROR (RGBcolor::scale): cannot scale RGBcolor by negative factor."));
            }
            r *= scale_factor;
            g *= scale_factor;
            b *= scale_factor;
        }
 
        friend std::ostream &operator<<(std::ostream &os, const helios::RGBcolor &c) {
            return os << "helios::RGBcolor<" << c.r << ", " << c.g << ", " << c.b << ">";
        }
 
        constexpr RGBcolor operator+(const RGBcolor &c) const noexcept;
        constexpr RGBcolor operator-(const RGBcolor &c) const noexcept;
        constexpr bool operator==(const RGBcolor &c) const noexcept;
        constexpr bool operator!=(const RGBcolor &c) const noexcept;
 
    private:
 
        static constexpr float clamp(float value) {
            return (value < 0.f) ? 0.f : (value > 1.f) ? 1.f : value;
        }
    };
 
 
    inline RGBcolor make_RGBcolor(float red, float green, float blue) {
        return {red, green, blue};
    }
 
 
    RGBcolor blend(const RGBcolor &color0, const RGBcolor &color1, float weight);
 
    constexpr RGBcolor RGBcolor::operator+(const RGBcolor &c) const noexcept {
        return {clamp(r + c.r), clamp(g + c.g), clamp(b + c.b)};
    }
 
    constexpr RGBcolor RGBcolor::operator-(const RGBcolor &c) const noexcept {
        return {clamp(r - c.r), clamp(g - c.g), clamp(b - c.b)};
    }
 
    constexpr bool RGBcolor::operator==(const RGBcolor &c) const noexcept {
        constexpr float epsilon = 1e-6f;
        auto constexpr_abs = [](float x) constexpr { return x < 0 ? -x : x; };
        return (constexpr_abs(r - c.r) < epsilon) && (constexpr_abs(g - c.g) < epsilon) && (constexpr_abs(b - c.b) < epsilon);
    }
 
    constexpr bool RGBcolor::operator!=(const RGBcolor &c) const noexcept {
        return c.r != r || c.g != g || c.b != b;
    }
 
    namespace RGB {
 
        extern RGBcolor red;
        extern RGBcolor blue;
        extern RGBcolor green;
        extern RGBcolor cyan;
        extern RGBcolor magenta;
        extern RGBcolor yellow;
        extern RGBcolor orange;
        extern RGBcolor violet;
        extern RGBcolor black;
        extern RGBcolor white;
        extern RGBcolor lime;
        extern RGBcolor silver;
        extern RGBcolor gray;
        extern RGBcolor navy;
        extern RGBcolor brown;
        extern RGBcolor khaki;
        extern RGBcolor greenyellow;
        extern RGBcolor forestgreen;
        extern RGBcolor yellowgreen;
        extern RGBcolor goldenrod;
 
    } // namespace RGB
 
 
    struct RGBAcolor {
 
 
        float r;
 
        float g;
 
        float b;
 
        float a;
 
        constexpr RGBAcolor() : r(0), g(0), b(0), a(1) {
        }
 
 
        constexpr RGBAcolor(float red, float green, float blue, float alpha) : r(RGBAcolor::clamp(red)), g(RGBAcolor::clamp(green)), b(RGBAcolor::clamp(blue)), a(RGBAcolor::clamp(alpha)) {
        }
 
 
        explicit constexpr RGBAcolor(const float C[4]) : r(RGBAcolor::clamp(C[0])), g(RGBAcolor::clamp(C[1])), b(RGBAcolor::clamp(C[2])), a(RGBAcolor::clamp(C[3])) {
        }
 
 
        explicit RGBAcolor(const std::vector<float> &C) {
            if (C.size() != 4) {
                throw(std::runtime_error("ERROR: cannot initialize RGBAcolor using supplied vector - size should be 4."));
            }
            r = RGBAcolor::clamp(C[0]);
            g = RGBAcolor::clamp(C[1]);
            b = RGBAcolor::clamp(C[2]);
            a = RGBAcolor::clamp(C[3]);
        }
 
 
        void scale(float scale_factor) {
            if (scale_factor < 0) {
                throw(std::runtime_error("ERROR (RGBAcolor::scale): cannot scale RGBAcolor by negative factor."));
            }
            r *= scale_factor;
            g *= scale_factor;
            b *= scale_factor;
            a *= scale_factor;
        }
 
        friend std::ostream &operator<<(std::ostream &os, const helios::RGBAcolor &c) {
            return os << "helios::RGBAcolor<" << c.r << ", " << c.g << ", " << c.b << ", " << c.a << ">";
        }
 
        constexpr RGBAcolor operator+(const RGBAcolor &c) const noexcept;
        constexpr RGBAcolor operator-(const RGBAcolor &c) const noexcept;
        constexpr bool operator==(const helios::RGBAcolor &c) const noexcept;
        constexpr bool operator!=(const helios::RGBAcolor &c) const noexcept;
 
    private:
 
        static constexpr float clamp(float value) {
            return (value < 0.f) ? 0.f : (value > 1.f) ? 1.f : value;
        }
    };
 
 
    inline RGBAcolor make_RGBAcolor(float r, float g, float b, float a) {
        return {r, g, b, a};
    }
 
 
    inline RGBAcolor make_RGBAcolor(RGBcolor color, float a) {
        return {color.r, color.g, color.b, a};
    }
 
 
    RGBAcolor blend(const helios::RGBAcolor &color0, const helios::RGBAcolor &color1, float weight);
 
    constexpr bool RGBAcolor::operator==(const RGBAcolor &c) const noexcept {
        constexpr float epsilon = 1e-6f;
        auto constexpr_abs = [](float x) constexpr { return x < 0 ? -x : x; };
        return (constexpr_abs(r - c.r) < epsilon) && (constexpr_abs(g - c.g) < epsilon) && (constexpr_abs(b - c.b) < epsilon && (constexpr_abs(a - c.a) < epsilon));
    }
 
    constexpr bool RGBAcolor::operator!=(const RGBAcolor &c) const noexcept {
        return c.r != r || c.g != g || c.b != b || c.a != a;
    }
 
    constexpr RGBAcolor RGBAcolor::operator+(const RGBAcolor &c) const noexcept {
        return {clamp(r + c.r), clamp(g + c.g), clamp(b + c.b), clamp(a + c.a)};
    }
 
    constexpr RGBAcolor RGBAcolor::operator-(const RGBAcolor &c) const noexcept {
        return {clamp(r - c.r), clamp(g - c.g), clamp(b - c.b), clamp(a - c.a)};
    }
 
    namespace RGBA {
 
        extern RGBAcolor red;
        extern RGBAcolor blue;
        extern RGBAcolor green;
        extern RGBAcolor cyan;
        extern RGBAcolor magenta;
        extern RGBAcolor yellow;
        extern RGBAcolor orange;
        extern RGBAcolor violet;
        extern RGBAcolor black;
        extern RGBAcolor white;
        extern RGBAcolor lime;
        extern RGBAcolor silver;
        extern RGBAcolor gray;
        extern RGBAcolor navy;
        extern RGBAcolor brown;
        extern RGBAcolor khaki;
        extern RGBAcolor greenyellow;
        extern RGBAcolor forestgreen;
        extern RGBAcolor yellowgreen;
        extern RGBAcolor goldenrod;
 
    } // namespace RGBA
 
 
    struct Date {
 
 
        int day;
        int month;
        int year;
 
        constexpr Date() : day(1), month(1), year(2000) {
        }
 
 
        Date(int day, int month, int year) : day(day), month(month), year(year) {
 
            if (day < 1 || day > 31) {
                throw(std::runtime_error("ERROR (Date constructor): Day of month is out of range (day of " + std::to_string(day) + " was given)."));
            } else if (month < 1 || month > 12) {
                throw(std::runtime_error("ERROR (Date constructor): Month of year is out of range (month of " + std::to_string(month) + " was given)."));
            } else if (year < 1000) {
                throw(std::runtime_error("ERROR (Date constructor): Year should be specified in YYYY format."));
            }
        }
 
 
        [[nodiscard]] int JulianDay() const;
 
        void incrementDay();
 
 
        [[nodiscard]] bool isLeapYear() const;
 
        constexpr bool operator==(const helios::Date &c) const noexcept;
        constexpr bool operator!=(const helios::Date &c) const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, helios::Date const &d) {
            return os << d.year << "-" << std::setfill('0') << std::setw(2) << d.month << "-" << std::setfill('0') << std::setw(2) << d.day;
        }
    };
 
 
    inline Date make_Date(int day, int month, int year) {
        if (day < 1 || day > 31) {
            throw(std::runtime_error("ERROR (make_Date): Day of month is out of range (day of " + std::to_string(day) + " was given)."));
        } else if (month < 1 || month > 12) {
            throw(std::runtime_error("ERROR (make_Date): Month of year is out of range (month of " + std::to_string(month) + " was given)."));
        } else if (year < 1000) {
            throw(std::runtime_error("ERROR (make_Date): Year should be specified in YYYY format."));
        }
 
        return {day, month, year};
    }
 
    constexpr bool Date::operator==(const Date &c) const noexcept {
        return c.day == day && c.month == month && c.year == year;
    }
 
    constexpr bool Date::operator!=(const Date &c) const noexcept {
        return c.day != day || c.month != month || c.year != year;
    }
 
 
    inline Date Julian2Calendar(int JulianDay, int year) {
 
        int month = -1;
 
        int skips_leap[] = {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335};
        int skips_nonleap[] = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334};
        const int *skips = nullptr;
 
        if (JulianDay < 1 || JulianDay > 366) {
            throw(std::runtime_error("ERROR (Julian2Calendar): Julian day of " + std::to_string(JulianDay) + " is out of range (should be >0 and <367)."));
        }
 
        if ((year - 2000) % 4 == 0) { // leap year
            skips = skips_leap;
        } else { // non-leap year
            skips = skips_nonleap;
        }
 
        for (int i = 1; i <= 12; i++) {
            if (i == 12) {
                month = 12;
            } else if (JulianDay > skips[i - 1] && JulianDay <= skips[i]) {
                month = i;
                break;
            }
        }
        assert(month > 0 && month < 13);
 
        int day = JulianDay - skips[month - 1];
 
        assert(day > 0 && day < 32);
 
        return {day, month, year};
    }
 
 
    inline int Calendar2Julian(Date date) {
 
        int skips_leap[] = {0, 31, 60, 91, 121, 152, 182, 214, 244, 274, 305, 335};
        int skips_nonleap[] = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334};
        int *skips = nullptr;
 
        if ((date.year - 2000) % 4 == 0) { // leap year
            skips = skips_leap;
        } else { // non-leap year
            skips = skips_nonleap;
        }
 
        int JD = skips[date.month - 1] + date.day;
        if (JD <= 0 || JD > 366) {
            throw("ERROR (Calendar2Julian): Julian day of " + std::to_string(JD) + " is out of range (should be >0 and <=366).");
        }
 
        return JD;
    }
 
 
    inline Date make_Date(int JulianDay, int year) {
        if (JulianDay < 1 || JulianDay > 366) {
            throw(std::runtime_error("ERROR (make_Date): Julian day is out of range (Julian day of " + std::to_string(JulianDay) + " was given)."));
        } else if (year < 1000) {
            throw(std::runtime_error("ERROR (make_Date): Year should be specified in YYYY format."));
        }
 
        return Julian2Calendar(JulianDay, year);
    }
 
 
 
    struct Time {
 
        int second;
        int minute;
        int hour;
 
        constexpr Time() : second(0), minute(0), hour(0) {
        }
 
 
        Time(int hour, int minute) : second(0), minute(minute), hour(hour) {
 
            if (minute < 0 || minute > 59) {
                throw(std::runtime_error("ERROR (Time constructor): Minute out of range (0-59)."));
            } else if (hour < 0 || hour > 23) {
                throw(std::runtime_error("ERROR (Time constructor): Hour out of range (0-23)."));
            }
        }
 
 
        Time(int hour, int minute, int second) : second(second), minute(minute), hour(hour) {
 
            if (second < 0 || second > 59) {
                throw(std::runtime_error("ERROR (Time constructor): Second out of range (0-59)."));
            } else if (minute < 0 || minute > 59) {
                throw(std::runtime_error("ERROR (Time constructor): Minute out of range (0-59)."));
            } else if (hour < 0 || hour > 23) {
                throw(std::runtime_error("ERROR (Time constructor): Hour out of range (0-23)."));
            }
        }
 
        constexpr bool operator==(const helios::Time &c) const noexcept;
        constexpr bool operator!=(const helios::Time &c) const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, helios::Time const &t) {
            return os << t.hour << ":" << std::setfill('0') << std::setw(2) << t.minute << ":" << std::setfill('0') << std::setw(2) << t.second;
        }
    };
 
 
    inline Time make_Time(int hour, int minute) {
        if (minute < 0 || minute > 59) {
            throw(std::runtime_error("ERROR (make_Time): Minute of " + std::to_string(minute) + " out of range (0-59)."));
        } else if (hour < 0 || hour > 23) {
            throw(std::runtime_error("ERROR (make_Time): Hour of " + std::to_string(hour) + " out of range (0-23)."));
        }
 
        return {hour, minute, 0};
    }
 
 
    inline Time make_Time(int hour, int minute, int second) {
        if (second < 0 || second > 59) {
            throw(std::runtime_error("ERROR (make_Time): Second of " + std::to_string(second) + " out of range (0-59)."));
        } else if (minute < 0 || minute > 59) {
            throw(std::runtime_error("ERROR (make_Time): Minute of " + std::to_string(minute) + " out of range (0-59)."));
        } else if (hour < 0 || hour > 23) {
            throw(std::runtime_error("ERROR (make_Time): Hour of " + std::to_string(hour) + " out of range (0-23)."));
        }
 
        return {hour, minute, second};
    }
 
    constexpr bool Time::operator==(const Time &c) const noexcept {
        return c.hour == hour && c.minute == minute && c.second == second;
    }
 
    constexpr bool Time::operator!=(const Time &c) const noexcept {
        return c.hour != hour || c.minute != minute || c.second != second;
    }
 
 
    struct Location {
 
        float latitude_deg;
        float longitude_deg;
        float UTC_offset;
 
        Location() {
            latitude_deg = 38.55f;
            longitude_deg = 121.76f;
            UTC_offset = 8;
        }
 
 
        Location(float latitude_deg, float longitude_deg, float UTC_offset) {
 
            this->latitude_deg = latitude_deg;
            this->longitude_deg = longitude_deg;
            this->UTC_offset = UTC_offset;
        }
 
        bool operator==(const helios::Location &c) const noexcept;
        bool operator!=(const helios::Location &c) const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, helios::Location const &t) {
            return os << "<" << t.latitude_deg << "," << t.longitude_deg << "," << t.UTC_offset << ">";
        }
    };
 
 
    inline Location make_Location(float latitude_deg, float longitude_deg, float UTC_offset) {
        return {latitude_deg, longitude_deg, UTC_offset};
    }
 
    inline bool Location::operator==(const Location &c) const noexcept {
        return c.latitude_deg == latitude_deg && c.longitude_deg == longitude_deg && c.UTC_offset == UTC_offset;
    }
 
    inline bool Location::operator!=(const Location &c) const noexcept {
        return c.latitude_deg != latitude_deg || c.longitude_deg != longitude_deg || c.UTC_offset != UTC_offset;
    }
 
 
    struct SphericalCoord {
    private:
        float elevation_private;
        float zenith_private;
 
    public:
        float radius;
 
        const float &elevation;
 
        const float &zenith;
        float azimuth;
 
 
        SphericalCoord() : elevation_private(0.5f * static_cast<float>(M_PI)), zenith_private(0.f), elevation(elevation_private), zenith(zenith_private) {
            radius = 1;
            azimuth = 0;
        }
 
        SphericalCoord(float radius, float elevation_radians, float azimuth_radians) :
            elevation_private(elevation_radians), zenith_private(0.5f * static_cast<float>(M_PI) - elevation_radians), radius(radius), elevation(elevation_private), zenith(zenith_private), azimuth(azimuth_radians) {
        }
 
        SphericalCoord(const SphericalCoord &c) : elevation_private(c.elevation_private), zenith_private(c.zenith_private), radius(c.radius), elevation(c.elevation_private), zenith(c.zenith_private), azimuth(c.azimuth) {
        }
 
        SphericalCoord &operator=(const SphericalCoord &c) noexcept {
            if (this != &c) {
                elevation_private = c.elevation_private;
                zenith_private = c.zenith_private;
                radius = c.radius;
                azimuth = c.azimuth;
            }
            return *this;
        }
 
        bool operator==(const helios::SphericalCoord &c) const noexcept;
        bool operator!=(const helios::SphericalCoord &c) const noexcept;
 
        friend std::ostream &operator<<(std::ostream &os, helios::SphericalCoord const &coord) {
            return os << "helios::SphericalCoord<" << coord.radius << ", " << coord.elevation << ", " << coord.azimuth << ">";
        }
    };
 
 
    inline SphericalCoord make_SphericalCoord(float elevation_radians, float azimuth_radians) {
        return {1, elevation_radians, azimuth_radians};
    }
 
 
    inline SphericalCoord make_SphericalCoord(float radius, float elevation_radians, float azimuth_radians) {
#ifdef HELIOS_DEBUG
        if (radius <= 0.f) {
            throw("ERROR (helios::make_SphericalCoord): radius should be positive, but " + std::to_string(radius) + " was given.");
        }
#endif
        return {radius, elevation_radians, azimuth_radians};
    }
 
    inline bool SphericalCoord::operator==(const SphericalCoord &c) const noexcept {
        return c.radius == radius && c.zenith == zenith && c.elevation == elevation && c.azimuth == azimuth;
    }
 
    inline bool SphericalCoord::operator!=(const SphericalCoord &c) const noexcept {
        return c.radius != radius || c.zenith != zenith || c.elevation != elevation || c.azimuth != azimuth;
    }
 
 
} // namespace helios
 
#endif