// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 
// POV-Ray SDL file for rendering natural cubic splines
// Made for easy translation to JustBasic
// Copyright 2014 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 

#version 3.7;

#include "colors.inc"

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

#declare Size = 5;
#declare N = Size - 1;

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

#declare Points =
  array[Size] {
    <-2, 0, 40>,
    <0, 0, -2>,
    <1, 0, 2>,
    <0, 6, 2>,
    <8, 2, 2>
  }

#declare Start = Times[0];
#declare End = Times[N];
#declare Span = End - Start;

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 
// The built in way
 
#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 Radius = 0.1;

#declare Intervals = 100;

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

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 
// Yet another way
// By using global variables - Not Good !

#declare X = array[Size];
#declare Y = array[Size];
#declare Z = array[Size];

#for (I, 0, N)
  #declare X[I] = Points[I].x;
  #declare Y[I] = Points[I].y;
  #declare Z[I] = Points[I].z;
#end // for


#declare tt = array[Size]; // Needed for calculating coefficients
#declare yy = array[Size]; // Needed for calculating coefficients

#declare hh = array[Size]; // Temporary use (The last value is not used)
#declare uu = array[Size]; // Temporary use (The first and last values are not used)
#declare vv = array[Size]; // Temporary use (The first and last values are not used)
#declare zz = array[Size]; // Temporary use

#declare kk = array[Size][3]; // The calculated coefficients


#macro CalculateCoefficients()

  #for (I, 0, N - 1)
    #declare hh[I] = tt[I+1] - tt[I];
    #declare kk[I][1] = 6*(yy[I+1] - yy[I])/hh[I];
  #end // for

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

  #declare zz[N] = 0;
  #for (I, N - 1, 1, -1)
    #declare zz[I] = (vv[I] - hh[I]*zz[I+1])/uu[I];
  #end // for
  #declare zz[0] = 0;

  #for (I, 0, N - 1)
    #declare kk[I][0] = (zz[I+1] - zz[I])/(6*hh[I]);
    #declare kk[I][1] = zz[I]/2;
    #declare kk[I][2] = -hh[I]/6*(zz[I+1] + 2*zz[I]) + (yy[I+1] - yy [I])/hh[I];
  #end // for

#end // macro CalculateCoefficients


// Copy Times to tt
#for (I, 0, N)
  #declare tt[I] = Times[I];
#end // for


// Calculate coefficients for the X polynomials

// Copy X to yy
#for (I, 0, N)
  #declare yy[I] = X[I];
#end // for

CalculateCoefficients()

#declare Coeffs_X = array[Size][3];

// Copy kk to Coeffs_X
#for (I, 0, N - 1)
  #declare Coeffs_X[I][0] = kk[I][0];
  #declare Coeffs_X[I][1] = kk[I][1];
  #declare Coeffs_X[I][2] = kk[I][2];
#end // for


// Calculate coefficients for the Y polynomials

// Copy Y to yy
#for (I, 0, N)
  #declare yy[I] = Y[I];
#end // for

CalculateCoefficients()

#declare Coeffs_Y = array[Size][3];

// Copy kk to Coeffs_Y
#for (I, 0, N - 1)
  #declare Coeffs_Y[I][0] = kk[I][0];
  #declare Coeffs_Y[I][1] = kk[I][1];
  #declare Coeffs_Y[I][2] = kk[I][2];
#end // for


// Calculate coefficients for the Z polynomials

// Copy Z to yy
#for (I, 0, N)
  #declare yy[I] = Z[I];
#end // for

CalculateCoefficients()

#declare Coeffs_Z = array[Size][3];

// Copy kk to Coeffs_Z
#for (I, 0, N - 1)
  #declare Coeffs_Z[I][0] = kk[I][0];
  #declare Coeffs_Z[I][1] = kk[I][1];
  #declare Coeffs_Z[I][2] = kk[I][2];
#end // for

// ===== 

#macro CubicPolynomial(T, A, B, C, D)

  (D + T*(C + T*(B + T*A))) // Equals (A*T^3 + B*T^2 + C*T + D)

#end // macro CubicPolynomial


#macro CalculateCoordinates(T, I)

    #declare dT = T - Times[I];
    #declare ComponentX =
      CubicPolynomial(
        dT,
        Coeffs_X[I][0],
        Coeffs_X[I][1],
        Coeffs_X[I][2],
        X[I]
      );
    #declare ComponentY =
      CubicPolynomial(
        dT,
        Coeffs_Y[I][0],
        Coeffs_Y[I][1],
        Coeffs_Y[I][2],
        Y[I]
      );
    #declare ComponentZ =
      CubicPolynomial(
        dT,
        Coeffs_Z[I][0],
        Coeffs_Z[I][1],
        Coeffs_Z[I][2],
        Z[I]
      );

#end // macro CalculateCoordinates


#declare Colors =
  array[Size] {
    color Cyan,
    color Green,
    color Orange,
    color Yellow,
    color Blue // Not used
  }

#declare Intervals = 400;

#declare Radius = 0.08;

#declare J = 0;
#while (J < Intervals)
  #declare T = Start + J/Intervals*Span;
  #declare ComponentX = 0;
  #declare ComponentY = 0;
  #declare ComponentZ = 0;

  #if (T < Times[0])
    CalculateCoordinates(T, 0)
    #declare Color = (Colors[0]/2 + White);
  #end // if

  #for (I, 0, N - 1)
    #if (Times[I] <= T & T < Times[I+1])
      CalculateCoordinates(T, I)
      #declare Color = Colors[I];
    #end // if
  #end // for

  #if (Times[N] <= T)
    CalculateCoordinates(T, N - 1)
    #declare Color = (Colors[N-1]/2 + White);
  #end // if

  sphere {
    <ComponentX, ComponentY, ComponentZ>, Radius
    pigment { color Color }
  }
  #declare J = J + 1;
#end // while

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

background { color <0.05, 0.10, 0.10> }

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

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

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