//#####################################################################
// Copyright 2002-2007, Robert Bridson, Ronald Fedkiw, Eran Guendelman, Geoffrey Irving, Craig Schroeder, Tamar Shinar.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class SPHERE
//##################################################################### 
#ifndef __SPHERE__
#define __SPHERE__

#include <Core/Basic_Geometry/BASIC_GEOMETRY_FORWARD.h>
#include <Core/Math_Tools/cube.h>
#include <Core/Math_Tools/pow.h>
#include <Core/Math_Tools/RANGE.h>
#include <Core/Matrices/SYMMETRIC_MATRIX.h>
#include <Core/Meshes/MESHES_FORWARD.h>
#include <numbers>
namespace PhysBAM{
using std::numbers::pi;
using std::numbers::pi_v;

template<class TV>
class SPHERE
{
    typedef typename TV::SCALAR T;
    enum WORKAROUND {d=TV::m};
public:
    typedef TV VECTOR_T;

    TV center;
    T radius;

    SPHERE()
        :radius(1)
    {}

    SPHERE(const TV& center,const T radius)
        :center(center),radius(radius)
    {}

    TV Normal(const TV& location) const
    {return (location-center).Normalized();}

    SYMMETRIC_MATRIX<T,TV::m> Hessian(const TV& X) const
    {TV z=X-center;T w=z.Normalize();return ((T)1-Outer_Product(z))/w;}

    bool Inside(const TV& location,const T thickness_over_two) const
    {return (location-center).Magnitude_Squared() <= sqr(radius-thickness_over_two);}

    bool Lazy_Inside(const TV& location) const
    {return (location-center).Magnitude_Squared() <= sqr(radius);}
    
    bool Outside(const TV& location,const T thickness_over_two) const
    {return (location-center).Magnitude_Squared() >= sqr(radius+thickness_over_two);}

    bool Lazy_Outside(const TV& location) const
    {return (location-center).Magnitude_Squared() >= sqr(radius);}

    bool Boundary(const TV& location,const T thickness_over_two) const
    {return !Inside(location,thickness_over_two) && !Outside(location,thickness_over_two);}

    TV Surface(const TV& location) const  
    {return radius*(location-center).Normalized()+center;}

    T Signed_Distance(const TV& location) const
    {return (location-center).Magnitude()-radius;}

    static T Size_Constant()
    {
        static const constexpr T unit_sphere_size[4]={0,2,pi_v<T>,pi_v<T>*4/3};
        return unit_sphere_size[d];
    }

    T Size() const
    {return Size_Constant()*pow<d>(radius);}

    RANGE<TV> Bounding_Box() const
    {return RANGE<TV>(center).Thickened(radius);}

    T Volume() const
    {return pi_v<T>*4/3*cube(radius);}
//#####################################################################
    void Sector_Volumes(const TV& origin,T volumes[1<<d],const T thickness_over_two=0) const;
    static std::string Name();
    T Circular_Segment_Area(const T h) const;
    template<class RW> void Read(UNTYPED_ISTREAM<RW>& input);
    template<class RW> void Write(UNTYPED_OSTREAM<RW>& output) const;
    void Fill_Uniform(RANDOM_NUMBERS<T>& rand,const RANGE<TV>& box);
//#####################################################################
};   
template<class T> class RAY_HIT;
template<class TV> class RAY;

template<class T,class TV> SPHERE<TV> Scaled(const SPHERE<TV>& b,const TV& center,T scale);
template<class TV> SPHERE<TV> Transformed(const SPHERE<TV>& b,const FRAME<TV>& f);

template<class T> TRIANGULATED_AREA<T>*
Generate_Triangles(const SPHERE<VECTOR<T,2> >& circle,int levels=4);
template<class T> TRIANGULATED_SURFACE<T>*
Generate_Triangles(const SPHERE<VECTOR<T,3> >& sphere,int levels=4);
template<class T> TRIANGULATED_SURFACE<T>*
Tessellate_Boundary(const SPHERE<VECTOR<T,3> >& sphere,int levels=4,bool exp_scale=true)
{
    PHYSBAM_ASSERT(exp_scale);
    return Generate_Triangles(sphere,levels);
}
template<class T> SEGMENTED_CURVE<VECTOR<T,2> >*
Tessellate_Boundary(const SPHERE<VECTOR<T,2> >& sphere,int levels=4,bool exp_scale=true);
template<class TV> RAY_HIT<typename TV::SCALAR>
Intersects(const RAY<TV>& ray,const SPHERE<TV>& sphere,typename TV::SCALAR isocontour=0);
template<class TV> VECTOR<RAY_HIT<typename TV::SCALAR>,2>
Intersection_Interval(const RAY<TV>& ray,const SPHERE<TV>& sphere,
    typename TV::SCALAR isocontour=0);
template<class TV,class T> void
All_Intersections(ARRAY<RAY_HIT<T> >& hits,const RAY<TV>& ray,const SPHERE<TV>& sphere,const T isocontour=0);
}
#endif
