utilities.h

Go to the documentation of this file.

/*
  Copyright (C) 2018-2024 by the authors of the World Builder code.

  This file is part of the World Builder.

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU Lesser General Public License as published
   by the Free Software Foundation, either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public License
   along with this program.  If not, see <https://www.gnu.org/licenses/>.
*/

#ifndef WORLD_BUILDER_UTILITIES_H
#define WORLD_BUILDER_UTILITIES_H


#include "world_builder/nan.h"
#include "world_builder/coordinate_systems/interface.h"
#include "world_builder/objects/natural_coordinate.h"
#include "world_builder/objects/bezier_curve.h"
#include <iostream>


namespace WorldBuilder
{

  namespace CoordinateSystems
  {
    class Interface;
  } // namespace CoordinateSystems
  namespace Utilities
  {

    inline bool approx(double a, double b, double error_factor=1e4)
    {
      return std::abs(a-b)<std::abs(std::min(a,b))*std::numeric_limits<double>::epsilon()*
             error_factor;
    }

    bool
    polygon_contains_point(const std::vector<Point<2> > &point_list,
                           const Point<2> &point);

    bool
    polygon_contains_point_implementation(const std::vector<Point<2> > &point_list,
                                          const Point<2> &point);


    double
    fraction_from_ellipse_center (const Point<2> &ellipse_center,
                                  const double semi_major_axis,
                                  const double eccentricity,
                                  const double rotation_angle,
                                  const Point<2> &point);


    double
    signed_distance_to_polygon(const std::vector<Point<2> > &point_list_,
                               const Point<2> &point_);


    std::array<double,3>
    cartesian_to_spherical_coordinates(const Point<3> &position);

    Point<3>
    spherical_to_cartesian_coordinates(const std::array<double,3> &scoord);

    std::array<double,3>
    cartesian_to_ellipsoidal_coordinates(const Point<3> &position,
                                         const double semi_major_axis_a,
                                         const double eccentricity);

    Point<3>
    ellipsoidal_to_cartesian_coordinates(const std::array<double,3> &phi_theta_d,
                                         const double semi_major_axis_a,
                                         const double eccentricity);

    CoordinateSystem
    string_to_coordinate_system (const std::string & /*coordinate_system*/);


    template<unsigned int dim>
    std::array<double,dim> convert_point_to_array(const Point<dim> &point);

    double
    string_to_double(const std::string &string);

    int
    string_to_int(const std::string &string);


    unsigned int
    string_to_unsigned_int(const std::string &string);


    Point<3> cross_product(const Point<3> &a, const Point<3> &b);

    enum class InterpolationType
    {
      None,
      Linear,
      MonotoneSpline,
      ContinuousMonotoneSpline,
      Invalid,
    };

    class interpolation
    {
      public:
        void set_points(const std::vector<double> &y);


        inline
        double operator() (const double x) const
        {
          if (x >= 0. && x <= static_cast<double>(mx_size_min))
            {
              const size_t idx = static_cast<size_t>(x);
              const double h = x-static_cast<double>(idx);
              return ((m[idx][0]*h + m[idx][1])*h + m[idx][2])*h + m[idx][3];
            }
          const size_t idx = std::min(static_cast<size_t>(std::max( static_cast<int>(x), static_cast<int>(0))),mx_size_min);
          const double h = x-static_cast<double>(idx);
          return (m[idx][1]*h + m[idx][2])*h + m[idx][3];
        }


        inline
        double operator() (const double x, const size_t idx, const double h) const
        {
          return (x >= 0. && x <= static_cast<double>(mx_size_min))
                 ?
                 ((m[idx][0]*h + m[idx][1])*h + m[idx][2])*h + m[idx][3]
                 :
                 (m[idx][1]*h + m[idx][2])*h + m[idx][3];
        }


        inline
        double value_inside (const size_t idx, const double h) const
        {
          WBAssert(idx <= mx_size_min, "Internal error: using value_inside outside the range of 0 to " << mx_size_min << ", but value was outside of this range: " << idx << ".");
          WBAssert(h >= 0 && h <= 1., "Internal error: using value_inside outside the range of 0 to " << mx_size_min << ", but value was outside of this range: " << h << ".");
          return ((m[idx][0]*h + m[idx][1])*h + m[idx][2])*h + m[idx][3];
        }


        inline
        double value_outside (const size_t idx, const double h) const
        {
          WBAssert(idx <= mx_size_min, "Internal error: using value_inside outside the range of 0 to " << mx_size_min << ", but value was outside of this range: " << idx << ".");
          WBAssert(!(static_cast<double>(idx) + h >= 0 && static_cast<double>(idx) + h <= 1.), "Internal error: using value_inside outside the range of 0 to " << mx_size_min << ", but value was outside of this range: " << static_cast<double>(idx) + h << " (h=" << h << ", idx = " << idx << ").");
          return (m[idx][1]*h + m[idx][2])*h + m[idx][3];
        }


        size_t mx_size_min;

        std::vector<std::array<double,4>> m; //m_a, m_b, m_c, m_y;

      private:
    };

    struct PointDistanceFromCurvedPlanes
    {
      PointDistanceFromCurvedPlanes(CoordinateSystem coordinate_system)
        :
        distance_from_plane(NaN::DSNAN),
        distance_along_plane(NaN::DSNAN),
        fraction_of_section(NaN::DSNAN),
        fraction_of_segment(NaN::DSNAN),
        section(NaN::ISNAN),
        segment(NaN::ISNAN),
        average_angle(NaN::DSNAN),
        depth_reference_surface(NaN::DSNAN),
        closest_trench_point(Point<3>(coordinate_system))
      {}

      double distance_from_plane;

      double distance_along_plane;

      double fraction_of_section;

      double fraction_of_segment;

      size_t section;

      size_t segment;

      double average_angle;

      double depth_reference_surface;

      Point<3> closest_trench_point;
    };

    PointDistanceFromCurvedPlanes distance_point_from_curved_planes(const Point<3> &check_point,
                                                                    const Objects::NaturalCoordinate &check_point_natural,
                                                                    const Point<2> &reference_point,
                                                                    const std::vector<Point<2> > &point_list,
                                                                    const std::vector<std::vector<double> > &plane_segment_lengths,
                                                                    const std::vector<std::vector<Point<2> > > &plane_segment_angles,
                                                                    const double start_radius,
                                                                    const std::unique_ptr<CoordinateSystems::Interface> &coordinate_system,
                                                                    const bool only_positive,
                                                                    const Objects::BezierCurve &bezier_curve);



    double wrap_angle(const double angle);

    double interpolate_angle_across_zero(const double angle_1,
                                         const double angle_2,
                                         const double fraction);

    std::array<double,3>
    euler_angles_from_rotation_matrix(const std::array<std::array<double,3>,3> &rotation_matrix);

    std::array<std::array<double,3>,3>
    euler_angles_to_rotation_matrix(double phi1, double theta, double phi2);

    std::string
    read_and_distribute_file_content(const std::string &filename);

    std::vector<double>
    calculate_ridge_distance_and_spreading(std::vector<std::vector<Point<2>>> mid_oceanic_ridges,
                                           std::vector<std::vector<double>> mid_oceanic_spreading_velocities,
                                           const std::unique_ptr<WorldBuilder::CoordinateSystems::Interface> &coordinate_system,
                                           const Objects::NaturalCoordinate &position_in_natural_coordinates_at_min_depth,
                                           const std::vector<std::vector<double>> &subducting_plate_velocities,
                                           const std::vector<double> &ridge_migration_times);

    // todo_effective
    std::vector<double>
    calculate_effective_trench_and_plate_ages(std::vector<double> ridge_parameters, double distance_along_plane);

    /*
     * Returns the result of the multiplication of two 3*3 matrix,
     * used in applying the random uniform distribution rotation matrix
     * to a given orientation (rotation matrix)
     */
    std::array<std::array<double,3>,3>
    multiply_3x3_matrices(const std::array<std::array<double,3>,3> mat1, std::array<std::array<double,3>,3> const mat2);

  } // namespace Utilities
} // namespace WorldBuilder


#endif