// 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 <iostream>
#include <vector>

#include "nvtkdisplay.h"
#include "nbsplinesurface.h"

// function that computes a geodesic on a surface from a point to another point
void geodesic(nparametricsurface &surf,npoint3 &dep_pt,npoint3 &arr_pt,std::vector<npoint3> &pts)
{
  // to code
  pts.push_back(dep_pt); // dep point is on the geodesic !
  pts.push_back(arr_pt); // arr point is on the geodesic !
  std::cout << pts.size() << " points on the geodesic" << std::endl;
}


//function to initialize the surface
void initsurface(nbsplinesurface &surf,double DX,double DY,double DZ);

int main(void)
{
  data_container data;
  nvtkdisplay display;

  std::vector<npoint3> glist; // points of the geodesic
  nbsplinesurface surface(6,3,6,3); // la surface
  initsurface(surface,0,0,0);

  npoint3 dep_point(0.25,0.5,0);// departure (u1,v1)
  npoint3 arr_point(2.5,2.9,0);// arrival (u2,v2)

  geodesic(surface,dep_point,arr_point,glist); // appel a la fonction calculant la geodesique

  if (glist.size()!=1)
    for (int i=1;i<glist.size();++i)
    {
      line l;
      npoint p;
      surface.P(glist[i-1].x(),glist[i-1].y(),p);
      l.pts[0]=p;
      surface.P(glist[i].x(),glist[i].y(),p);
      l.pts[1]=p;
      data.add_line(l);
    }
  else
  {
    npoint p;
    surface.P(glist[0].x(),glist[0].y(),p);
    data.add_point(p);
  }
  surface.Display(data);
  display.init_data(data);
  display.display();
  return 0;
}



//function to initialize the surface (implementation)
void initsurface(nbsplinesurface &surf,double DX,double DY,double DZ)
{
  for (int i=0;i<surf.nb_CP_u();++i)
  {
    for (int j=0;j<surf.nb_CP_v();++j)
    {
      double R=1.5;
      double x=i*2/((double)surf.nb_CP_u()-1)-1;
      double y=j*2/((double)surf.nb_CP_v()-1)-1;
      double rr=x+y*y;
      double z=sqrt(R*R-rr)-0.9;
      z=z*z*2;
      npoint p(x+DX,y+DY,z+DZ,1.0);
      surf.CP(i,j)=p;
    }
  }
  int endsu=surf.degree_u()+1;
  int du=surf.degree_u();
  int cpsu=surf.nb_CP_u()-1;
  for (int i=0;i<endsu;++i)
  {
    surf.u(i)=0.0;
    surf.u(cpsu+du+1-i)=(cpsu+du+2-2*endsu+1);
  }
  for (int i=1;i<=(cpsu+du+2-2*endsu);++i)
  {
    surf.u(i+endsu-1)=i;
  }
  int endsv=surf.degree_v()+1;
  int dv=surf.degree_v();
  int cpsv=surf.nb_CP_v()-1;

  for (int i=0;i<endsv;++i)
  {
    surf.v(i)=0.0;
    surf.v(cpsv+dv+1-i)=(cpsv+dv+2-2*endsv+1);
  }
  for (int i=1;i<=(cpsv+dv+2-2*endsv);++i)
  {
    surf.v(i+endsv-1)=i;
  }

}