// 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 <fstream>

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

#include "linear_algebra.h"

// function to init profiles
void init_profile(nbspline &profile);
void init_trajectory(nbspline &traj);


// fonction that computes the final surface
void extrude(nbspline &profile,nbspline &traj,nbsplinesurface &surf)
{

}




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

  nbspline profile(20,3); // need an even numbre of CPs
  nbspline trajectory(10,3);
  nbsplinesurface ext_surface; // skin of the profiles

  init_profile(profile);
  init_trajectory(trajectory);

  // display the profiles
  profile.Display(data);
  trajectory.Display(data);
  extrude(profile,trajectory,ext_surface); // function call to compute the skin surface

  ext_surface.Display(data);
  display.init_data(data);
  display.display();
  return 0;
}

void init_profile(nbspline &profile)
{
  // profil in the  x-y plane, basepoint = frame origin
  
  int degree=profile.degree();
  int nCP=profile.nb_CP();


  //nodal sequence
  for (int i=0;i<degree+1;++i)
  {
    profile.u(i)=0;
    profile.u(i+nCP)=nCP-degree;
  }
  for (int i=1;i<nCP-degree;++i)
  {
    profile.u(i+degree)=i;
  }

  double c=10.; //chord
  double t=0.15; //thickness
  // control points
  for (int i=0;i<nCP/2;++i)
  {
    double x=c*(nCP/2-i)/(nCP/2);
    double y=(t*c/0.2)*(0.2969*sqrt(x/c)-0.1260*(x/c)-0.3516*(x*x/(c*c))+0.2843*(x*x*x/(c*c*c))-0.1015*(x*x*x*x/(c*c*c*c)));
    profile.CP(i)=npoint(x,y,0.0,1.0);
    profile.CP(nCP-1-i)=npoint(x,-y,0.0,1.0);
  }
}

void init_trajectory(nbspline& traj)
{
    int degree=traj.degree();
  int nCP=traj.nb_CP();
  //nodal sequence
  for (int i=0;i<degree+1;++i)
  {
    traj.u(i)=0;
    traj.u(i+nCP)=nCP-degree;
  }
  for (int i=1;i<nCP-degree;++i)
  {
    traj.u(i+degree)=i;
  }
  double c=10.;
  for (int i=0;i<nCP;++i)
  {
    double z= (c*i)/(nCP-1);
    double y=z*z/c;
    traj.CP(i)=npoint(0.,y,z,1.0);
  }
}