text: /*
   * (C) Copyright 2013 ECMWF.
   *
   * This software is licensed under the terms of the Apache Licence Version 2.0
   * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
   * In applying this licence, ECMWF does not waive the privileges and immunities
   * granted to it by virtue of its status as an intergovernmental organisation
   * nor does it submit to any jurisdiction.
   */
  
  #pragma once
  
  #include <memory>
  #include <sstream>
  #include <string>
  
  #include "atlas/runtime/Exception.h"
  #include "atlas/util/Earth.h"
  #include "atlas/util/Object.h"
  #include "atlas/util/ObjectHandle.h"
  #include "atlas/util/Point.h"
  
  namespace atlas {
  namespace geometry {
  
  //------------------------------------------------------------------------------------------------------
  
  namespace detail {
  
  class GeometryBase : public util::Object {
  public:
      virtual ~GeometryBase()                                             = default;
      virtual void lonlat2xyz( const Point2&, Point3& ) const             = 0;
      virtual void xyz2lonlat( const Point3&, Point2& ) const             = 0;
      virtual double distance( const Point2& p1, const Point2& p2 ) const = 0;
      virtual double distance( const Point3& p1, const Point3& p2 ) const = 0;
      virtual double radius() const                                       = 0;
      virtual double area() const                                         = 0;
  
      Point3 xyz( const Point2& lonlat ) const {
          Point3 xyz;
          lonlat2xyz( lonlat, xyz );
          return xyz;
      }
      Point2 lonlat( const Point3& xyz ) const {
          Point2 lonlat;
          xyz2lonlat( xyz, lonlat );
          return lonlat;
      }
  };
  
  //------------------------------------------------------------------------------------------------------
  
  template <typename SphereT>
  class GeometrySphereT : public GeometryBase {
  public:
      void lonlat2xyz( const Point2& lonlat, Point3& xyz ) const override {
          SphereT::convertSphericalToCartesian( lonlat, xyz );
      }
      void xyz2lonlat( const Point3& xyz, Point2& lonlat ) const override {
          SphereT::convertCartesianToSpherical( xyz, lonlat );
      }
      double distance( const Point2& p1, const Point2& p2 ) const override { return SphereT::distance( p1, p2 ); }
      double distance( const Point3& p1, const Point3& p2 ) const override { return SphereT::distance( p1, p2 ); }
      double radius() const override { return SphereT::radius(); }
      double area() const override { return SphereT::area(); }
  };
  
  class GeometrySphere : public GeometryBase {
      using Sphere = eckit::geometry::Sphere;
  
  public:
      GeometrySphere( double radius ) : radius_( radius ) {}
      void lonlat2xyz( const Point2& lonlat, Point3& xyz ) const override {
          Sphere::convertSphericalToCartesian( radius_, lonlat, xyz );
      }
      void xyz2lonlat( const Point3& xyz, Point2& lonlat ) const override {
          Sphere::convertCartesianToSpherical( radius_, xyz, lonlat );
      }
      double distance( const Point2& p1, const Point2& p2 ) const override { return Sphere::distance( radius_, p1, p2 ); }
      double distance( const Point3& p1, const Point3& p2 ) const override { return Sphere::distance( radius_, p1, p2 ); }
      double radius() const override { return radius_; }
      double area() const override { return Sphere::area( radius_ ); }
  
  private:
      double radius_;
  };
  
  }  // namespace detail
  }  // namespace geometry
  
  //------------------------------------------------------------------------------------------------------
  
  class Geometry : public util::ObjectHandle<geometry::detail::GeometryBase> {
  public:
      using Handle::Handle;
  
      Geometry() : Handle( build<geometry::detail::GeometrySphereT<util::Earth>>() ) {}
      Geometry( const std::string& name ) : Handle( build( name ) ) {}
      Geometry( const char* name ) : Handle( build( name ) ) {}
      Geometry( double radius ) : Handle( build<geometry::detail::GeometrySphere>( radius ) ) {}
  
      template <typename SphereT>
      Geometry( const SphereT& ) : Handle( build<SphereT>() ) {}
  
      Point3 xyz( const Point2& lonlat ) const { return get()->xyz( lonlat ); }
      Point2 lonlat( const Point3& xyz ) const { return get()->lonlat( xyz ); }
      void xyz2lonlat( const Point3& xyz, Point2& lonlat ) const { get()->xyz2lonlat( xyz, lonlat ); }
      void lonlat2xyz( const Point2& lonlat, Point3& xyz ) const { get()->lonlat2xyz( lonlat, xyz ); }
      double distance( const Point2& p1, const Point2& p2 ) const { return get()->distance( p1, p2 ); }
      double distance( const Point3& p1, const Point3& p2 ) const { return get()->distance( p1, p2 ); }
      double radius() const { return get()->radius(); }
      double area() const { return get()->area(); }
  
  protected:
      template <typename GeometryT, typename... Args>
      static Implementation* build( Args... args ) {
          return new GeometryT( args... );
      }
  
      static Implementation* build( const std::string& name ) {
          // Factory without self registration
          if ( name == "Earth" ) {
              return build<geometry::detail::GeometrySphereT<util::Earth>>();
          }
          else if ( name == "UnitSphere" ) {
              return build<geometry::detail::GeometrySphereT<eckit::geometry::UnitSphere>>();
          }
          else {
              ATLAS_THROW_EXCEPTION( "name " << name << " is not a valid key for a Geometry" );
          }
      }
  };
  
  //------------------------------------------------------------------------------------------------------
  
  namespace geometry {
  using Earth = Geometry;  // Sphere with util::Earth radius by default
  class UnitSphere : public Geometry {
  public:
      UnitSphere() : Geometry( /*radius*/ 1. ) {}
  };
  
  }  // namespace geometry
  
  
  // ------------------------------------------------------------------
  // C wrapper interfaces to C++ routines
  
  extern "C" {
  Geometry::Implementation* atlas__Geometry__new_name( const char* name );
  Geometry::Implementation* atlas__Geometry__new_radius( const double radius );
  void atlas__Geometry__delete( Geometry::Implementation* This );
  void atlas__Geometry__xyz2lonlat( Geometry::Implementation* This, const double x, const double y, const double z,
                                    double& lon, double& lat );
  void atlas__Geometry__lonlat2xyz( Geometry::Implementation* This, const double lon, const double lat, double& x,
                                    double& y, double& z );
  double atlas__Geometry__distance_lonlat( Geometry::Implementation* This, const double lon1, const double lat1,
                                           const double lon2, const double lat2 );
  double atlas__Geometry__distance_xyz( Geometry::Implementation* This, const double x1, const double y1, const double z1,
                                        const double x2, const double y2, const double z2 );
  double atlas__Geometry__radius( Geometry::Implementation* This );
  double atlas__Geometry__area( Geometry::Implementation* This );
  }
  
  //------------------------------------------------------------------------------------------------------
  
  }  // namespace atlas
  