// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 
// FrenetFramesAlongParametricCurve.pov
// POV-Ray SDL file for illustrating Frenet frames along curve
// Copyright 2016 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.6;

#include "colors.inc"

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

#macro CurveTransform(FnX, FnY, FnZ, D1FnX, D1FnY, D1FnZ, D2FnX, D2FnY, D2FnZ, Time)

  #local pP = <FnX(Time), FnY(Time), FnZ(Time)>;
  #local vD1 = <D1FnX(Time), D1FnY(Time), D1FnZ(Time)>;
  #local vD2 = <D2FnX(Time), D2FnY(Time), D2FnZ(Time)>;

  #local vT = vnormalize(vD1); // Tangent vector
  #local vB = vnormalize(vcross(vD1, vD2)); // Binormal vector
  #local vN = vcross(vT, vB); // Normal vector

  transform {
    matrix <
      vT.x, vT.y, vT.z,
      vB.x, vB.y, vB.z,
      vN.x, vN.y, vN.z,
      pP.x, pP.y, pP.z
    >
  }

#end // macro CurveTransform


#macro ShowCurve(xFn, yFn, zFn, MinT, MaxT, ResT, WrapT, Radius)

  #if (ResT = 0)
    #local dT = 0;
  #else
    #local dT = (MaxT - MinT)/ResT;
  #end // if

  #if (dT != 0)
    #local I = 0;
    #local pT = <xFn(MinT), yFn(MinT), zFn(MinT)>;
    #while (I < ResT)
      #local pF = pT;
      #local NextI = I + 1;
      #local NextT = MinT + NextI*dT;
      #local pT = <xFn(NextT), yFn(NextT), zFn(NextT)>;
      sphere  { pF, Radius }
      cylinder { pF, pT, Radius }
      #local I = NextI;
    #end // while
    #if (!WrapT)
      sphere { pT, Radius*2 }
    #end // if
  #else
    sphere { <xFn(MinT), yFn(MinT), zFn(MinT)>, Radius }
  #end // if

#end // macro ShowCurve


// Numerical approximation of the first derivative

#macro FirstDerivativeOfFunction(Fn)

  #local H = 1E-4;
  #local H2 = H*2;

  function(U) { (Fn(U + H) - Fn(U - H))/H2 }

#end // FirstDerivativeOfFunction


// Numerical approximation of the second derivative

#macro SecondDerivativeOfFunction(Fn)

  #local H = 1E-4;
  #local HH = pow(H, 2);

  function(U) { (Fn(U + H) - 2*Fn(U) + Fn(U - H))/HH }

#end // SecondDerivativeOfFunction


#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

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

//   x(t) = +cos(t) +2*cos(2*t)
//   y(t) = +2*sin(3*t)
//   z(t) = +sin(t) -2*sin(2*t)

//  x'(t) = -sin(t) -4*sin(2*t)
//  y'(t) = +6*cos(3.t)
//  z'(t) = +cos(t) -4*cos(2*t)

// x''(t) = -cos(t) -8*cos(2*t)
// y''(t) = -18*sin(3*t)
// z''(t) = -sin(t) +8*sin(2*t)

// The functions for the curve
#declare xFn = function(_t) { +cos(_t) +2*cos(2*_t) }
#declare yFn = function(_t) { +2*sin(3*_t) }
#declare zFn = function(_t) { +sin(_t) -2*sin(2*_t) }

// The first derivatives of the functions for the curve
#declare D1_xFn = function(_t) { -sin(_t) -4*sin(2*_t) }
#declare D1_yFn = function(_t) { +6*cos(3*_t) }
#declare D1_zFn = function(_t) { +cos(_t) -4*cos(2*_t) }

// The second derivatives of the functions for the curve
#declare D2_xFn = function(_t) { -cos(_t) -8*cos(2*_t) }
#declare D2_yFn = function(_t) { -18*sin(3*_t) }
#declare D2_zFn = function(_t) { -sin(_t) +8*sin(2*_t) }

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 
// Alternative way of finding the derivitives of the functions 

/*
#undef D1_xFn
#undef D1_yFn
#undef D1_zFn

#declare D1_xFn = FirstDerivativeOfFunction(xFn)
#declare D1_yFn = FirstDerivativeOfFunction(yFn)
#declare D1_zFn = FirstDerivativeOfFunction(zFn)

#undef D2_xFn
#undef D2_yFn
#undef D2_zFn

#declare D2_xFn = SecondDerivativeOfFunction(xFn)
#declare D2_yFn = SecondDerivativeOfFunction(yFn)
#declare D2_zFn = SecondDerivativeOfFunction(zFn)
*/

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 
// Show the parametric curve

#declare NrOfSteps = 4*360;
#declare Step = 2*pi/NrOfSteps;

#declare Radius = 0.08;

#declare ClosedLoop = true;

#declare Curve =
  union {
    ShowCurve(xFn, yFn, zFn, 0, 2*pi, NrOfSteps, ClosedLoop, Radius)
    pigment { color White }
  }

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 
// Arrows to show the axes in a coordinate system

#declare StemRadius = 0.1;

#declare SmallArrows =
  union {
    sphere {
      <0, 0, 0>, 2*StemRadius
      pigment { color White }
    }
    object {
      Arrow(0*x, 4*x, StemRadius)
      pigment { color Red }
    }
    object {
      Arrow(0*y, 4*y, StemRadius)
      pigment { color Green }
    }
    object {
      Arrow(0*z, 4*z, StemRadius)
      pigment { color Blue }
    }
    scale <1, 1, 1>/5
  }

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 
// Triangle to go into the Frenet frames

#declare A = 0.4;

#declare p0 = vrotate(A*z, -120*x);
#declare p1 = vrotate(A*z,    0*x);
#declare p2 = vrotate(A*z, +120*x);

#declare R = 0.02;

#declare Triangle =
  union {
    union {
      sphere { p0, 1.5*R }
      sphere { p1, 1.5*R }
      sphere { p2, 1.5*R }
      pigment { color Gray }
    }
    cylinder {
      p0 ,p1, R
      pigment { color rgb <0.3, 0.4, 0.7> }
    }
    cylinder {
      p1 ,p2, R
      pigment { color rgb <0.7, 0.3, 0.4> }
    }
    cylinder {
      p2 ,p0, R
      pigment { color rgb <0.4, 0.7, 0.3> }
    }
  }

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 
// Polygon to go into the Frenet frames

#declare Rad = 0.5;

#declare CircleFnX = function(_t) { 0 }
#declare CircleFnY = function(_t) { Rad*cos(_t) }
#declare CircleFnZ = function(_t) { Rad*sin(_t) }

#declare Polygon =
  union {
    ShowCurve(
      CircleFnX, CircleFnY, CircleFnZ,
      0, 2*pi, 6,
      true,
      0.03
    )
    pigment { color (Blue + Gray)/2 }
  }

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10 
// Put something into the Frenet frames

#declare Shape = Triangle;
// #declare Shape = SmallArrows // Try to uncomment this
// #declare Shape = Polygon // - or this

#declare ShapesAroundCurve =
  union {
    #declare I = 0;
    #while (I < NrOfSteps)
      object {
        Shape
        CurveTransform(
          xFn, yFn, zFn,
          D1_xFn, D1_yFn, D1_zFn,
          D2_xFn, D2_yFn, D2_zFn,
          I*Step
        )
      }
      #declare I = I + 4;
    #end
  }

#declare SpiralAroundCurve =
  union {
    #declare I = 0;
    #while (I < NrOfSteps)
      sphere {
        vrotate(0.1*z, 30*degrees(I*Step)*x), 0.02
        CurveTransform(
          xFn, yFn, zFn,
          D1_xFn, D1_yFn, D1_zFn,
          D2_xFn, D2_yFn, D2_zFn,
          I*Step
        )
      }
      #declare I = I + 1;
    #end
    pigment { color 2*White }
  }

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

union {
  object { Curve }
  object { ShapesAroundCurve }
//  object { SpiralAroundCurve } // Try to uncomment this
  rotate 30*y
}

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

background { color Black }

light_source {
  <100, 100, 0>
  color Gray70
  shadowless
  rotate -120*y
}

light_source {
  <100, 100, 0>
  color Gray70
  shadowless
}

light_source {
  <100, 100, 0>
  color Gray70
  shadowless
  rotate 120*y
}

camera {
  orthographic
  location 8*y
  look_at 0*y
}

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