// 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 "nspline.h"
#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;


void initcircle(nspline &s)
{
  int degree=s.nb_CP();
  ofstream of("cp.txt");
  for (int i=0;i<degree;i++)
  {
    double delta=0;
//    if ( i> ( ( degree-1 ) /2 ) ) delta=1;
    s.u(i) =i+delta;
    double angle=i*n_pi/ (degree-1.);
    double random=2* (randr()-0.5) /100;
    random=0;
//  if (i<(degree-1)/2) random=0.05; else random=0;
//  if (i>(degree-1)/2) random=-0.05;
    if (i== (degree-1) /2) random=-1;
    npoint val;
    val[0]=cos(angle) * (1+random);
    val[1]=sin(angle) * (1+random);
    val[2]=0.;
    val[3]=1.;
    of << s.u(i) << " " << val[0] << " " << val[1]  << endl;
    s.CP(i) =val;
  }
  of.close();
}



void initcircle2(double pos[], double val[][2], int degree)
{
 ofstream of("cp.txt",ios::app);
 for (int i=0;i<degree;i++)
 {
  pos[i]=(i)/(degree-1.);
  double angle=i*n_pi/(degree-1.)+n_pi/4;
  double random=2*(randr()-0.5)/100;
  random=0;
//  if (i<(degree-1)/2) random=0.05; else random=0;
//  if (i>(degree-1)/2) random=-0.05;
  if (i==(degree-1)/2) random=-1;
  val[i][0]=cos(angle)*(1+random)/2+0.25;
  val[i][1]=sin(angle)*(1+random);
  of << pos[i] << " " << val[i][0] << " " << val[i][1]  << endl;
 }
 of.close();
}


void testlagrange(void)
{
 double ret[2];
 int deg=11;
 double val[100][2];
 double pos[100];

 initcircle2(pos,val,deg);
 double eps=(pos[deg-2]-pos[0])/(100*(nbpts-1));
 for (double u=pos[0];u<=(pos[nbpts-1]+eps/2);u+=eps)
 {

 lagrangeInterpolatingPolynomial (pos, val, deg, u, ret);
 cout << u << " "<< ret[0] << " "<< ret[1] << " " ;
 for (int i=0;i<deg;++i)
  cout << lagrangeBasis(pos,deg, i,u) << " " ;
 cout << endl;
 }

}

void testspline(void)
{
  npoint ret,ret2;
  int nbpts=5;
  nspline s(nbpts);
  initcircle(s);
  s.compute_natural();
//cardinalspline(val,der,nbpts,0.0);
// FDspline(val,der,pos,nbpts);
// initcircle2(pos2,val2,nbpts);
// naturalspline(val2,der2,nbpts);
//cardinalspline(val2,der2,nbpts,0.25);
// FDspline(val2,der2,pos2,nbpts);

// cardinalspline(val,der,nbpts,1.00);
// FDspline(val,der,pos,nbpts);
  double eps= (s.u(nbpts-1)-s.u(0)) / (10* (nbpts-1));
  for (double u=s.u(0);u<= (s.u(nbpts-1) +eps/2);u+=eps)
  {

    s.P(u,ret);
    ret.perspective_divide();
// evalspline(pos, val2, der2, u,nbpts, ret2);
    cout << u << " "<< ret[0] << " "<< ret[1] ; // << " " << ret2[0] << " "<< ret2[1] << " " << (ret[0]*sin(pi/4)-ret[1]*cos(pi/4))/2 +0.25<< " "<< ret[0]*cos(pi/4)+ret[1]*sin(pi/4) << " ";
    cout << endl;
  }
}

int main(void)
{
  testspline();
  return 0;
}