// 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 "nsurface.h"
#include "nbeziercurve.h"
#include <cmath>

void nparametricsurface::dPdu(double u_,double v_,npoint& ret) const  // first derivative - generic numerical implementation
{
    double eps=(max_u()-min_u())*1e-5;
    npoint P1,P2;
    if (u_<=min_u()+eps/2) u_+=eps/2;
    if (u_>=max_u()-eps/2) u_-=eps/2;
    P(u_-eps/2,v_,P1);
    P(u_+eps/2,v_,P2);
    npoint3 p1(P1);
    npoint3 p2(P2);
    npoint3 r((p2-p1)/eps);
    ret=r;
    ret.w()=0.;
}


void nparametricsurface::d2Pdu2(double u_,double v_,npoint& ret) const  // second derivative - generic numerical implementation
{
    double eps=(max_u()-min_u())*1e-5;
    npoint P1,P2,P3;
    if (u_<=min_u()+eps) u_+=eps;
    if (u_>=max_u()-eps) u_-=eps;
    P(u_-eps,v_,P1);
    P(u_,v_,P2);
    P(u_+eps,v_,P3);
    npoint3 p1(P1);
    npoint3 p2(P2);
    npoint3 p3(P3);
    npoint3 r((p1+p3-p2*2.0)/(eps*eps));
    ret=r;
    ret.w()=0.;
}

void nparametricsurface::dPdv(double u_,double v_,npoint& ret) const  // first derivative - generic numerical implementation
{
    double eps=(max_v()-min_v())*1e-5;
    npoint P1,P2;
    if (v_<=min_v()+eps/2) v_+=eps/2;
    if (v_>=max_v()-eps/2) v_-=eps/2;
    P(u_,v_-eps/2,P1);
    P(u_,v_+eps/2,P2);
    npoint3 p1(P1);
    npoint3 p2(P2);
    npoint3 r((p2-p1)/eps);
    ret=r;
    ret.w()=0.;
}

void nparametricsurface::d2Pdv2(double u_,double v_,npoint& ret) const  // second derivative - generic numerical implementation
{
    double eps=(max_v()-min_v())*1e-5;
    npoint P1,P2,P3;
    if (v_<=min_v()+eps) v_+=eps;
    if (v_>=max_v()-eps) v_-=eps;
    P(u_,v_-eps,P1);
    P(u_,v_,P2);
    P(u_,v_+eps,P3);
    npoint3 p1(P1);
    npoint3 p2(P2);
    npoint3 p3(P3);
    npoint3 r((p1+p3-p2*2.0)/(eps*eps));
    ret=r;
    ret.w()=0.;
}


void nparametricsurface::d2Pdudv(double u_,double v_,npoint& ret) const  // second derivative - generic numerical implementation
{
    double epsu=(max_u()-min_u())*1e-5;
    double epsv=(max_v()-min_v())*1e-5;
    npoint P1,P2,P3,P4;
    if (u_<=min_u()+epsu/2) u_+=epsu/2;
    if (u_>=max_u()-epsu/2) u_-=epsu/2;
    if (v_<=min_v()+epsv/2) v_+=epsv/2;
    if (v_>=max_v()-epsv/2) v_-=epsv/2;
    P(u_-epsu/2,v_-epsv/2,P1);
    P(u_-epsu/2,v_+epsv/2,P2);
    P(u_+epsu/2,v_-epsv/2,P3);
    P(u_+epsu/2,v_+epsv/2,P4);
    npoint3 p1(P1);
    npoint3 p2(P2);
    npoint3 p3(P3);
    npoint3 p4(P4);
    npoint3 r((p4-p2-p3+p1)/(epsu*epsv));
    ret=r;
    ret.w()=0.;
}


void nparametricsurface::Display(data_container& data) const    // default : uses P(u,res)
{
  for (int i=0;i< nb_CP();i++)
    data.add_point(CP(i));
  const int nbptsu=10*nb_CP_u() +100;
  const int nbptsv=10*nb_CP_v() +100;
  std::vector<npoint> dep(nbptsv+1);
  std::vector<npoint> arr(nbptsv+1);
  bool first=true;
  double epsu= (max_u()-min_u()) / (nbptsu);
  double epsv= (max_v()-min_v()) / (nbptsv);
  for (double u=min_u();u<= (max_u() +epsu/2);u+=epsu)
  {
    if (first)
    {
      int i=0;
      for (double v=min_v();v<= (max_v() +epsv/2);v+=epsv)
      {
        P(u,v,arr[i++]);
      }
      first=false;
    }
    else
    {
      dep=arr;
      int i=0;
      for (double v=min_v();v<= (max_v() +epsv/2);v+=epsv)
      {
        P(u,v,arr[i++]);
      }
      quad q;
      for (int j=0;j< (i-1);++j)
      {
        q.pts[0]=dep[j];
        q.pts[1]=arr[j];
        q.pts[2]=arr[j+1];
        q.pts[3]=dep[j+1];
        data.add_quad(q);
      }
//      return;
    }
  }
}