// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 
// POV-Ray SDL file for even spacing of points along natural cubic splines
// Copyright 2010 by Tor Olav Kristensen
// Email: t o r . o l a v . k [_A_T_] g m a i l . c o m
// http://subcube.com
// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 

#include "colors.inc"
#include "functions.inc"

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 

#macro Arrow(pStart, pEnd, Radius)

  #local A = Radius*3;
  #local B = Radius*2;
  #local C = Radius;
  #local vDirection = vnormalize(pEnd - pStart); 
  
  union {
    difference {
      cone { <0, 0, 0>, 0, -(A + B)*vDirection, Radius + C }
      cone { -A*vDirection, Radius, -(A + 2*B)*vDirection, Radius + 2*C }
      translate pEnd
    }
    cylinder { pStart, pEnd - A*vDirection, Radius }
  }

#end // macro Arrow



#macro Size1A(Array1A)

  dimension_size(Array1A, 1)

#end // macro Size1A



#macro ExtractComponent1A(Vectors1A, v0)

  #local SizeU = Size1A(Vectors1A);
  #local Component1A = array[SizeU];

  #local I = 0;
  #while (I < SizeU)
    #local Component1A[I] = vdot(Vectors1A[I], v0);
    #local I = I + 1;
  #end // while

  Component1A

#end // macro ExtractComponent1A

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 

#macro TangentFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn)

  #local TempFn =
    function(d1_X, d1_Y, d1_Z) {
      d1_X/
      f_r(d1_X, d1_Y, d1_Z)
    }

  function(u_) {
    TempFn(
      d1_X_Fn(u_),
      d1_Y_Fn(u_),
      d1_Z_Fn(u_)
    )
  }

#end // macro TangentFunctionX_1A



#macro TangentFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn)

  #local TempFn =
    function(d1_X, d1_Y, d1_Z) {
      d1_Y/
      f_r(d1_X, d1_Y, d1_Z)
    }

  function(u_) {
    TempFn(
      d1_X_Fn(u_),
      d1_Y_Fn(u_),
      d1_Z_Fn(u_)
    )
  }

#end // macro TangentFunctionY_1A



#macro TangentFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn)

  #local TempFn =
    function(d1_X, d1_Y, d1_Z) {
      d1_Z/
      f_r(d1_X, d1_Y, d1_Z)
    }

  function(u_) {
    TempFn(
      d1_X_Fn(u_),
      d1_Y_Fn(u_),
      d1_Z_Fn(u_)
    )
  }

#end // macro TangentFunctionZ_1A



#macro BinormalFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn)

  #local TempFn =
    function(d1_X, d1_Y, d1_Z, d2_X, d2_Y, d2_Z) {
      (d1_Y*d2_Z - d1_Z*d2_Y)/
      f_r(
        d1_Y*d2_Z - d1_Z*d2_Y,
        d1_Z*d2_X - d1_X*d2_Z,
        d1_X*d2_Y - d1_Y*d2_X
      )
    }

  function(u_) {
    TempFn(
      d1_X_Fn(u_),
      d1_Y_Fn(u_),
      d1_Z_Fn(u_),
      d2_X_Fn(u_),
      d2_Y_Fn(u_),
      d2_Z_Fn(u_)
    )
  }

#end // macro BinormalFunctionX_1A



#macro BinormalFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn)

  #local TempFn =
    function(d1_X, d1_Y, d1_Z, d2_X, d2_Y, d2_Z) {
      (d1_Z*d2_X - d1_X*d2_Z)/
      f_r(
        d1_Y*d2_Z - d1_Z*d2_Y,
        d1_Z*d2_X - d1_X*d2_Z,
        d1_X*d2_Y - d1_Y*d2_X
      )
    }

  function(u_) {
    TempFn(
      d1_X_Fn(u_),
      d1_Y_Fn(u_),
      d1_Z_Fn(u_),
      d2_X_Fn(u_),
      d2_Y_Fn(u_),
      d2_Z_Fn(u_)
    )
  }

#end // macro BinormalFunctionY_1A



#macro BinormalFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn)

  #local TempFn =
    function(d1_X, d1_Y, d1_Z, d2_X, d2_Y, d2_Z) {
      (d1_X*d2_Y - d1_Y*d2_X)/
      f_r(
        d1_Y*d2_Z - d1_Z*d2_Y,
        d1_Z*d2_X - d1_X*d2_Z,
        d1_X*d2_Y - d1_Y*d2_X
      )
    }

  function(u_) {
    TempFn(
      d1_X_Fn(u_),
      d1_Y_Fn(u_),
      d1_Z_Fn(u_),
      d2_X_Fn(u_),
      d2_Y_Fn(u_),
      d2_Z_Fn(u_)
    )
  }

#end // macro BinormalFunctionZ_1A



#macro NormalFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn)

  #local B_FnY_1A = BinormalFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn);
  #local B_FnZ_1A = BinormalFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn);
  #local T_FnY_1A = TangentFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn);
  #local T_FnZ_1A = TangentFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn);

  function(u_) {
    B_FnY_1A(u_)*T_FnZ_1A(u_) - B_FnZ_1A(u_)*T_FnY_1A(u_)
  }

#end // macro NormalFunctionX_1A



#macro NormalFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn)

  #local B_FnX_1A = BinormalFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn);
  #local B_FnZ_1A = BinormalFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn);
  #local T_FnX_1A = TangentFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn);
  #local T_FnZ_1A = TangentFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn);

  function(u_) {
    B_FnZ_1A(u_)*T_FnX_1A(u_) - B_FnX_1A(u_)*T_FnZ_1A(u_)
  }

#end // macro NormalFunctionY_1A



#macro NormalFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn)

  #local B_FnX_1A = BinormalFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn);
  #local B_FnY_1A = BinormalFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn);
  #local T_FnX_1A = TangentFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn);
  #local T_FnY_1A = TangentFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn);

  function(u_) {
    B_FnX_1A(u_)*T_FnY_1A(u_) - B_FnY_1A(u_)*T_FnX_1A(u_)
  }

#end // macro NormalFunctionZ_1A

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 

#macro IntegrateTrapezoidFunction(Fn, A, B, M)

  #local D = (B - A)/M;

  function(n) {
    select(n - 1, 0, ((Fn(A) + Fn(A + n*D))/2 + sum(i, 1, n - 1, Fn(A + i*D)))*D)
  }

#end // macro IntegrateTrapezoidFunction

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 

#macro NaturalCubicSplineCoefficients1A3A(tt, yy)

  #local Size = Size1A(tt);
  #local N = Size - 1;
  #local hh = array[Size];
  #local uu = array[Size];
  #local bb = array[Size];
  #local vv = array[Size];
  #local zz = array[Size];
  #local kk = array[Size][3];

  #local I = 0;
  #while (I <= N - 1)
    #local hh[I] = tt[I+1] - tt[I];
    #local kk[I][1] = 6*(yy[I+1] - yy[I])/hh[I];
    #local I = I + 1;
  #end // while

  #local uu[1] = 2*(hh[0] + hh[1]);
  #local vv[1] = kk[1][1] - kk[0][1];
  #local I = 2;
  #while (I <= N - 1)
    #local uu[I] = 2*(hh[I] + hh[I-1]) - pow(hh[I-1], 2)/uu[I-1];
    #local vv[I] = kk[I][1] - kk[I-1][1] - hh[I-1]*vv[I-1]/uu[I-1];
    #local I = I + 1;
  #end // while

  #local zz[N] = 0;
  #local I = N - 1;
  #while (I >= 1)
    #local zz[I] = (vv[I] - hh[I]*zz[I+1])/uu[I];
    #local I = I - 1;
  #end // while
  #local zz[0] = 0;

  #local I = 0;
  #while (I <= N - 1)
    #local kk[I][0] = (zz[I+1] - zz[I])/(6*hh[I]);
    #local kk[I][1] = zz[I]/2;
    #local kk[I][2] = -hh[I]/6*(zz[I+1] + 2*zz[I]) + (yy[I+1] - yy [I])/hh[I];
    #local I = I + 1;
  #end // while

  kk

#end // macro NaturalCubicSplineCoefficients1A3A




#macro NaturalCubicSplineFunction1A(tt, yy, kk)

  #local Size = Size1A(tt);
  #local N = Size - 1;
  #local ff = array[Size];

  #local I = 0;
  #while (I <= N - 1)
    #local ff[I] =
      function(t_) {
        yy[I] + (t_ - tt[I])*(kk[I][2] + (t_ - tt[I])*(kk[I][1] + (t_ - tt[I])*kk[I][0]))
      }
    #local I = I + 1;
  #end // while

  function(t_) {
    select(t_ - tt[0], 1, 0, 0)*ff[0](t_) +
    #local I = 0;
    #while (I <= N - 1)
      select(t_ - tt[I], 0, select(t_ - tt[I+1], 1, 0))*ff[I](t_) +
      #local I = I + 1;
    #end // while
    select(t_ - tt[N], 0, 1, 1)*ff[N-1](t_)
  }

#end // macro NaturalCubicSplineFunction1A



#macro Der1_NaturalCubicSplineFunction1A(tt, kk)

  #local Size = Size1A(tt);
  #local N = Size - 1;
  #local ff = array[Size];

  #local I = 0;
  #while (I <= N - 1)
    #local ff[I] =
      function(t_) {
        kk[I][2] + (t_ - tt[I])*(2*kk[I][1] + (t_ - tt[I])*3*kk[I][0])
      }
    #local I = I + 1;
  #end // while

  function(t_) {
    select(t_ - tt[0], 1, 0, 0)*ff[0](t_) +
    #local I = 0;
    #while (I <= N - 1)
      select(t_ - tt[I], 0, select(t_ - tt[I+1], 1, 0))*ff[I](t_) +
      #local I = I + 1;
    #end // while
    select(t_ - tt[N], 0, 1, 1)*ff[N-1](t_)
  }

#end // macro Der1_NaturalCubicSplineFunction1A
        



#macro Der2_NaturalCubicSplineFunction1A(tt, kk)

  #local Size = Size1A(tt);
  #local N = Size - 1;
  #local ff = array[Size];

  #local I = 0;
  #while (I <= N - 1)
    #local ff[I] =
      function(t_) {
        2*kk[I][1] + (t_ - tt[I])*6*kk[I][0]
      }
    #local I = I + 1;
  #end // while

  function(t_) {
    select(t_ - tt[0], 1, 0, 0)*ff[0](t_) +
    #local I = 0;
    #while (I <= N - 1)
      select(t_ - tt[I], 0, select(t_ - tt[I+1], 1, 0))*ff[I](t_) +
      #local I = I + 1;
    #end // while
    select(t_ - tt[N], 0, 1, 1)*ff[N-1](t_)
  }

#end // macro Der2_NaturalCubicSplineFunction1A

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 

#declare Times =
  array[5] {
     -2,
     -1,
      0,
      0.5,
      1
  }

#declare Points =
  array[5] {
    <-2, 0, 40>,
    <0, 0, -2>,
    <1, 0, 2>,
    <0, 6, 2>,
    <8, 2, 2>
  }
    
#declare NaturalSpline =
  spline {
    natural_spline
    Times[0], Points[0]
    Times[1], Points[1]
    Times[2], Points[2]
    Times[3], Points[3]
    Times[4], Points[4]
  }


#declare Start = -2;
#declare End = 1;
#declare Span = End - Start;


#declare Intervals = 200;

#declare I = 0;
#while (I < Intervals)
  #declare T = Start + I/Intervals*Span;
  sphere {
    NaturalSpline(T), 0.04
    pigment { color White }
  }
  #declare I = I + 1;
#end // while

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 

#declare CompX = ExtractComponent1A(Points, x);
#declare CompY = ExtractComponent1A(Points, y);
#declare CompZ = ExtractComponent1A(Points, z);

#declare CoeffsX = NaturalCubicSplineCoefficients1A3A(Times, CompX)
#declare CoeffsY = NaturalCubicSplineCoefficients1A3A(Times, CompY)
#declare CoeffsZ = NaturalCubicSplineCoefficients1A3A(Times, CompZ)

#declare FnX = NaturalCubicSplineFunction1A(Times, CompX, CoeffsX);
#declare FnY = NaturalCubicSplineFunction1A(Times, CompY, CoeffsY);
#declare FnZ = NaturalCubicSplineFunction1A(Times, CompZ, CoeffsZ);

#declare Der1_FnX = Der1_NaturalCubicSplineFunction1A(Times, CoeffsX);
#declare Der1_FnY = Der1_NaturalCubicSplineFunction1A(Times, CoeffsY);
#declare Der1_FnZ = Der1_NaturalCubicSplineFunction1A(Times, CoeffsZ);

#declare Der2_FnX = Der2_NaturalCubicSplineFunction1A(Times, CoeffsX);
#declare Der2_FnY = Der2_NaturalCubicSplineFunction1A(Times, CoeffsY);
#declare Der2_FnZ = Der2_NaturalCubicSplineFunction1A(Times, CoeffsZ);


#declare T_FnX = TangentFunctionX_1A(Der1_FnX, Der1_FnY, Der1_FnZ);
#declare T_FnY = TangentFunctionY_1A(Der1_FnX, Der1_FnY, Der1_FnZ);
#declare T_FnZ = TangentFunctionZ_1A(Der1_FnX, Der1_FnY, Der1_FnZ);

#declare B_FnX = BinormalFunctionX_1A(Der1_FnX, Der1_FnY, Der1_FnZ, Der2_FnX, Der2_FnY, Der2_FnZ);
#declare B_FnY = BinormalFunctionY_1A(Der1_FnX, Der1_FnY, Der1_FnZ, Der2_FnX, Der2_FnY, Der2_FnZ);
#declare B_FnZ = BinormalFunctionZ_1A(Der1_FnX, Der1_FnY, Der1_FnZ, Der2_FnX, Der2_FnY, Der2_FnZ);

#declare N_FnX = NormalFunctionX_1A(Der1_FnX, Der1_FnY, Der1_FnZ, Der2_FnX, Der2_FnY, Der2_FnZ);
#declare N_FnY = NormalFunctionY_1A(Der1_FnX, Der1_FnY, Der1_FnZ, Der2_FnX, Der2_FnY, Der2_FnZ);
#declare N_FnZ = NormalFunctionZ_1A(Der1_FnX, Der1_FnY, Der1_FnZ, Der2_FnX, Der2_FnY, Der2_FnZ);


#declare Arrows =
  union {
    object { Arrow(0*x, 4*x, 0.05) pigment { color Red } }
    object { Arrow(0*y, 4*y, 0.05) pigment { color Green } }
    object { Arrow(0*z, 4*z, 0.05) pigment { color Blue } }
  }

Arrows


#declare Intervals = 300;

#declare LengthFn =
  IntegrateTrapezoidFunction(
    function(t_) { f_r(Der1_FnX(t_), Der1_FnY(t_), Der1_FnZ(t_)) }
    Start, End,
    Intervals
  )

#declare TotalLength = LengthFn(Intervals);


#declare Knots = Intervals + 1;
#declare tt_ = array[Knots]
#declare yy_ = array[Knots]

#declare I = 0;
#while (I < Knots)
  #declare tt_[I] = Start + I/Intervals*Span;
  #declare yy_[I] = Start + LengthFn(I)/TotalLength*Span;
  #declare I = I + 1;
#end // while

#declare kk_ = NaturalCubicSplineCoefficients1A3A(yy_, tt_)
#declare InverseSplineFn = NaturalCubicSplineFunction1A(yy_, tt_, kk_)


#declare SmallArrows =
  object {
    Arrows
    scale <1, 1, 1>/6
  }

#declare Intervals = 400;
#declare Knots = Intervals + 1;

#declare I = 0;
#while (I < Knots)
  #declare T_ = Start + I/Intervals*Span;
  #declare T = InverseSplineFn(T_);
//  #declare T = T_; // Uncomment this line to see it without correction
    object {
      SmallArrows
      matrix <
        T_FnX(T), T_FnY(T), T_FnZ(T),
        B_FnX(T), B_FnY(T), B_FnZ(T),
        N_FnX(T), N_FnY(T), N_FnZ(T),
          FnX(T),   FnY(T),   FnZ(T)
      >
    }
  #declare I = I + 1;
#end // while

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 

light_source { <-1, 6, -2>*5 color White shadowless }

camera {
  location <1, 1, -1>*8
  look_at <2, 1, 2>
}

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10