// Gnurbs - A curve and surface library
// Copyright (C) 2008-2026 Eric Bechet
//
// See the LICENSE file for contributions and license information.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//

#include "ncurve.h"
#include "nsurface.h"
#include <cmath>
#include <limits>

double ncurve::FPeps=std::numeric_limits<double>::epsilon();
double ncurve::FPeps2=sqrt(std::numeric_limits<double>::epsilon());
double ncurve::FPeps3=pow(std::numeric_limits<double>::epsilon(),1./3.);
double ncurve::FPeps4=pow(std::numeric_limits<double>::epsilon(),1./4.);


void ncurve::update_FPeps(bool mode)
{
  if (mode)
  {
    FPeps2=FPeps;
    FPeps3=FPeps;
    FPeps4=FPeps;
  }
  else
  {
    FPeps2=sqrt(FPeps);
    FPeps3=pow(FPeps,1./3.);
    FPeps4=pow(FPeps,1./4.);
  }

}


void nparametriccurve::dPdu(double u_,npoint& ret) const  // first derivative - generic numerical implementation
{
  double h=FPeps3*fabs(u_);

  npoint P0;
  P(u_,P0);
  npoint3 p0(P0);
  double norm=p0.norm();
  double hh=FPeps3*norm;
  double minu=min_u();
  double maxu=max_u();
  if (hh>h) h=hh;
  if (h==0.) h=1e-10;
  volatile double up=u_ + h;
  volatile double um=u_ - h;
  if (!is_periodic() && u_>=minu && u_<=maxu) // check that we are inside bounds, only if u is actually in bounds and not periodic.
  {
    if (um<minu) um=minu;
    if (up>maxu) up=maxu;
  }
  double eps=up-um;
  npoint P1,P2;
  P(um,P1);
  P(up,P2);
  npoint3 p1(P1);
  npoint3 p2(P2);
  npoint3 r((p2-p1)/eps);
  ret=r;
  ret.w()=0;
}
  
void nparametriccurve::d2Pdu2(double u_,npoint& ret) const // second derivative, generic numerical implementation
{
  double h=FPeps4*fabs(u_)*100.;
  npoint P0;
  P(u_, P0);
  npoint3 p0(P0);
  double norm=p0.norm();
  double hh=FPeps4*norm*100.;
  double minu=min_u();
  double maxu=max_u();
  if (hh>h) h=hh;
  if (h==0.) h=1e-10;
  volatile double up=u_ + h;
  volatile double um=u_ - h;
  if (!is_periodic() && u_>=minu && u_<=maxu) // check that we are inside bounds, only if u is actually in bounds and not periodic.
  {
    if (um<minu)
    {
      um=minu;
      u_=um+h;
      up=um+2*h;
    }
    if (up>maxu)
    {
      up=maxu;
      u_=up-h;
      um=up-2*h;
    }
  }
  double eps=(up-um)/2;
  npoint P1,P2,P3;
  P(um,P1);
  P(u_,P2);
  P(up,P3);
  npoint3 p1(P1);
  npoint3 p2(P2);
  npoint3 p3(P3);
  npoint3 r((p1-p2*2.0+p3)/(eps*eps));
  ret=r;
  ret.w()=0;
}

void nparametriccurve::fundamental_forms(double u, double& M1, double& M2, npoint3* J) const
{
  npoint P;
  npoint3 P_u,T;
  npoint3 P_uu,N;
  
  dPdu(u,P);P[3]=1.0; P_u=P;
  if (J) *J=P_u;
  d2Pdu2(u,P);P[3]=1.0; P_uu=P;
  M1 = P_u.dotprod(P_u);
  T=P_u;
//  P_u.print(std::cout);
//  P_uu.print(std::cout);
  double l=T.normalize();
//        T.print(std::cout);
  if (l>1e-10)
  {
    npoint3 ori;
    ori.crossprod(T,P_uu);
    double norm=ori.normalize();
    if (norm>1e-10)
    {
      N.crossprod(ori,T);
      N.normalize();
//      N.print(std::cout);
      M2=P_uu.dotprod(N);
    }
    else M2=0.;
  }
  else M2=0.;
}


void nparametriccurve::principal_curvatures(double u,double &kappa,double *m1, double *m2,double thr) const //Calcule les directions principales à partir du point pt sur la surface
{
  double M1;
  double M2;
  fundamental_forms(u,M1,M2);
  if (m1) *m1=M1;
  if (m2) *m2=M2;
  if (M1>thr)
    kappa=M2/M1;
  else
    kappa=0.;
}
  


void nparametriccurve::_Display(double from,double to, data_container& data) const    // default version : uses P(u)
{
  properties p,p1;
  p=data.getproppoints();
  p1=p;
  if (p.nb_uid<2) p.nb_uid=2;
  for (int i=0;i< nb_CP();i++)
  {
    p.uids[1]=i;
    data.setproppoints(p);
    data.add_point(CP(i));
  }
  data.setproppoints(p1);
  
  std::vector< std::pair<npoint3,double> > pts;

  npoint res1,res2;
  P(from,res1);
  P(to,res2);
  pts.push_back(std::pair<npoint3,double> (res1,from));
  int n=0;
  if (nb_CP())
    n=ceil(log(nb_CP()) /log(2.));
  if (n<3) n=3;
//  std::cout << n << " " << nb_CP() << std::endl;
  _Discretize(pts, from, res1, to,res2, 0, 0, 0, 0.005 ,n);
  pts.push_back(std::pair<npoint3,double> (res2,to));
  for (int i=0;i<pts.size()-1;++i)
  {
    line l;
    l.pts[0]=pts[i].first;
    l.pts[1]=pts[i+1].first;
    data.add_line(l);
  }
}

void nparametriccurve::Display(data_container& data) const    // default version : uses P(u)
{
  _Display(min_u(),max_u(),data);
}

void nparametriccurve::Discretize(nmesh& mesh, double du, double ds, double eps, double epsr, int n) const
{
  std::vector< std::pair<npoint3,double> > pts;
  npoint res1,res2;
  P(min_u(),res1);
  P(max_u(),res2);
  int nn=1;
  if (n)
    nn= (log(n) /log(2.));
  if (nn<0) nn=0;
  _Discretize(pts,min_u(),res1,max_u(),res2,du,ds,eps,epsr,nn);
  int deb=mesh.add_node(nmesh::nnode(-1,res1,npoint3(min_u(),0.,0.)));
  int arr;
  for (int i=0;i<pts.size();++i)
  {
    arr=mesh.add_node(nmesh::nnode(-1,pts[i].first,npoint3(pts[i].second,0.,0.)));
    nmesh::nelement el(nmesh::LINE);
    el.pts[0]=deb;
    el.pts[1]=arr;
    mesh.add_element(el);
    deb=arr;
  }
  arr=mesh.add_node(nmesh::nnode(-1,res2,npoint3(max_u(),0.,0.)));
  nmesh::nelement el(nmesh::LINE);
  el.pts[0]=deb;
  el.pts[1]=arr;
  mesh.add_element(el);
}

void nparametriccurve::Discretize(nmesh& mesh, nmesh& bnd, double du, double ds, double eps, double epsr, int n) const
{
  std::vector< std::pair<npoint3,double> > pts;
  npoint res1,res2;
  P(bnd.get_node(0).xyz[0],res1);
  P(bnd.get_node(1).xyz[0],res2);
  int nn=1;
  if (n)
    nn= (log(n) /log(2.));
  if (nn<0) nn=0;
  _Discretize(pts,bnd.get_node(0).xyz[0],res1,bnd.get_node(1).xyz[0],res2,du,ds,eps,epsr,nn);
  int deb=mesh.add_node(nmesh::nnode(bnd.get_node(0).num,res1,npoint3(bnd.get_node(0).xyz[0],0.,0.)));
  int arr;
  for (int i=0;i<pts.size();++i)
  {
    arr=mesh.add_node(nmesh::nnode(-1,pts[i].first,npoint3(pts[i].second,0.,0.)));
    nmesh::nelement el(nmesh::LINE);
    el.pts[0]=deb;
    el.pts[1]=arr;
    mesh.add_element(el);
    deb=arr;
  }
  arr=mesh.add_node(nmesh::nnode(bnd.get_node(1).num,res2,npoint3(bnd.get_node(1).xyz[0],0.,0.)));
  nmesh::nelement el(nmesh::LINE);
  el.pts[0]=deb;
  el.pts[1]=arr;
  mesh.add_element(el);
}

void nparametriccurve::_Discretize(std::vector< std::pair< npoint3, double > >& pts, double from, npoint3 Posf, double to, npoint3 Post, double du, double ds, double eps, double epsr, int n) const
{
  if ((from==0.) && (to==0.))
  {
    from=min_u();
    to=max_u();
  }
  if ((du==0.) && (ds==0.) && (eps==0.) && (epsr==0.) && (n==0)) return;
  double mid= (from+to) /2;
  npoint Poswm;
  P(mid,Poswm);
  npoint3 Posm(Poswm);
//  npoint3 Der=(Post-Posf)/(to-from);
//  npoint3 Curv=2*(Post+Posf-2*Posm)/(to-from);
  double loceps= (Posm-0.5* (Post+Posf)).norm();
  double locds= ((Posm-Post).norm() + (Posf-Posm).norm());
  double locepsr=loceps/locds;
  double locdu=fabs(to-from);
  if (((locdu>du) && (du>0.)) || ((locds>ds) && (ds>0.)) || ((loceps>eps) && (eps>0.)) || ((locepsr>epsr) && (epsr>0.)) || (n>0))
  {
    if (n>-10)
    {
      _Discretize(pts,from,Posf,mid,Posm,du,ds,eps,epsr,n-1);
      pts.push_back(std::pair< npoint3, double > (Posm,mid));
      _Discretize(pts,mid,Posm,to,Post,du,ds,eps,epsr,n-1);
    }
  }
}


void nembeddedcurve::Puv(double u,npoint& ret) const     // point along the curve for parameter u (parametric space)
{
  crv.P(u,ret);
}


void nembeddedcurve::P(double u,npoint& ret) const     // point along the curve for parameter u
{
  npoint P;
  crv.P(u,P);
  surf.P(P[0]/P[3],P[1]/P[3],ret);
}



void nembeddedcurve::relative_curvatures(double t,double kappa[2], double m1[2][2], double m2[2][2]) const  //Computes normal and geodesic curvatures of an embedded curve
{
  npoint P,P_u,P_v,P_uu,P_vv,P_uv;
  npoint P_t,P_tt;
  npoint3 Pt,dPdt,d2Pdt2;
  npoint3 Puv,dPdu,dPdv,d2Pdu2,d2Pdudv,d2Pdv2;
  double M1[2][2];
  double M2[2][2];
  crv.d2P(t,P,P_t,P_tt);
  P[3]=1.0;P_t[3]=1.0;P_tt[3]=1.0;
  Pt=P;dPdt=P_t;d2Pdt2=P_tt;
  surf.d2P(Pt[0],Pt[1],P,P_u,P_v,P_uu,P_vv,P_uv);
  P[3]=1.0;P_u[3]=1.0;P_v[3]=1.0;P_uu[3]=1.0;P_vv[3]=1.0;P_uv[3]=1.0;
  Puv=P;dPdu=P_u;dPdv=P_v;d2Pdu2=P_uu;d2Pdv2=P_vv;d2Pdudv=P_uv;
  npoint3 dgdt(dPdt[0]*dPdu+dPdt[1]*dPdv);
  npoint3 d2gdt2(d2Pdt2[0]*dPdu+d2Pdt2[1]*dPdv+dPdt[0]*dPdt[0]*d2Pdu2+2*dPdt[0]*dPdt[1]*d2Pdudv+dPdt[1]*dPdt[1]*d2Pdv2);
  npoint3 n;
  n.crossprod(dPdu,dPdv);n.normalize();
  if(m1)
  {
   m1[0][0]=M1[0][0] = dPdu.dotprod(dPdu);
   m1[1][0] = m1[0][1] = M1[1][0] = M1[0][1] = dPdu.dotprod(dPdv);
   m1[1][1]= M1[1][1] = dPdv.dotprod(dPdv);
  }
  else
  {
   M1[0][0] = dPdu.dotprod(dPdu);
   M1[1][0] = M1[0][1] = dPdu.dotprod(dPdv);
   M1[1][1] = dPdv.dotprod(dPdv);
  }  
  if(m2)
  {
    m2[0][0] = M2[0][0]=n.dotprod(d2Pdu2);
    m2[0][1] = m2[1][0] = M2[0][1] = M2[1][0] = n.dotprod(d2Pdudv);
    m2[1][1] = M2[1][1] =n.dotprod(d2Pdv2);
  }
  else
  {
    M2[0][0]=n.dotprod(d2Pdu2);
    M2[0][1] = M2[1][0] = n.dotprod(d2Pdudv);
    M2[1][1] =n.dotprod(d2Pdv2);
  }
  
  double V[2]={dPdt[0],dPdt[1]};
  double ds2=multi(V,M1,V);

  double ds=sqrt(ds2);
  npoint3 ndgdt;
  ndgdt.crossprod(n,dgdt);
  kappa[0]=d2gdt2.dotprod(n)/(ds2);
  kappa[1]=d2gdt2.dotprod(ndgdt)/(ds*ds2);
/*  
  double tot=0.;
  double norm=0.;
  double geo=0.;
  double nn=multi(V,M2,V);  
principal_curvatures(t,tot);
  norm=nn/ds2;
  geo=sqrt(tot*tot-norm*norm);
  std::cout << "  Ktot=" << tot << " Knorm=" << norm << " Kg=" << geo << std::endl;
  std::cout << "* Ktot=" << sqrt(kappa[0]*kappa[0]+kappa[1]*kappa[1]) << " Knorm=" << kappa[0] << " Kg=" << kappa[1] << std::endl;
  */  
/*  
  {
    double kappa[2];
    npoint3 J[2];
  npoint P;
  crv.P(t,P);
  npoint3 pp(P);
  crv.dPdu(t,P);
  double V[2]={P[0],P[1]};
  double u=pp[0];
  double v=pp[1];
  double (*M1)[2];
  double (*M2)[2];
  if (m1) M1=m1; else M1=new double[2][2];
  if (m2) M2=m2; else M2=new double[2][2];
  surf.fundamental_forms(u,v,M1,J,M2);
  double a=multi(V,M1,V);
  if (a>1e-10)
  {
    double b=multi(V,M2,V);
    kappa[0]=b/a;
  }
  else kappa[0]=0;
  double totk;
  principal_curvatures(t,totk);
  double diff=totk*totk-kappa[0]*kappa[0]; // should be >=0 if no numerical errors
  if (diff<0) 
  {
    kappa[1]=0.; return; 
  }
  else 
    kappa[1]=sqrt(fabs(diff));
  // now define the sign of kappa[1] w/r of the surface parametrization.
    
 std::cout << "old way " << t << "  K= " << totk << " k_n= " << kappa[0] << " k_g= " << kappa[1] << std::endl;
  }*/
}