text: import math
  import numpy as np
  from mathutils import MatrixFactory
  
  
  # a curve consists of a function, a domain [min,max], and a resolution (divisions).
  # getPoints() returns equally spaced points along the curve
  # getFrames() returns { "tangents" : list, "normals" : list, "binormals" : list } at same set of points
  
  class Curve(object):
  
      # divisions = number of line segments. len(points) = divisions+1.
      def __init__(self, function, tMin=0, tMax=1, divisions=4, arcLengthDivisions=256):
          self.f = function
          self.tMin = tMin
          self.tMax = tMax
          self.divisions = divisions
          
          # reparameterize curve by arclength
          self.arcLengthDivisions = arcLengthDivisions
          # list of length of curve from tMin to t
          # approximate using # of points from arcLengthDivisions
          self.arcLengths = None
          self.totalArcLength = None
          self.calculateArcLengths()
          
      def calculateArcLengths(self):
          # note: points stored as numpy arrays for simpler calculation syntax
          arcLengthPoints = []
          deltaT = (self.tMax - self.tMin) / (self.arcLengthDivisions-1) 
          
          for index in range(self.arcLengthDivisions):
              t = self.tMin + index * deltaT
              point = self.getPointAt(t)
              arcLengthPoints.append( np.array(point) )
              
          self.arcLengths = [0]
          for index in range(1, self.arcLengthDivisions):
              arcLength = np.linalg.norm( arcLengthPoints[index] - arcLengthPoints[index-1] )
              self.arcLengths.append( self.arcLengths[index-1] + arcLength )
          
          self.totalArcLength = self.arcLengths[self.arcLengthDivisions-1]
          
      # u - percentage of total curve distance traversed, in [0,1]
      # t - corresponding value of variable from original parameterization
      def convert_u_to_t(self, u):
          arcLength = u * self.totalArcLength
          # index = np.searchsorted(list, value) means that list[index-1] < value <= list[index]
          index = np.searchsorted(self.arcLengths, arcLength)
          # arcLengthA < arcLength < arcLengthB
          arcLengthA = self.arcLengths[index-1]
          arcLengthB = self.arcLengths[index]
          segmentLength = arcLengthB - arcLengthA
          # determine how far along the segment this arclength would be
          percent = (arcLength - arcLengthA) / segmentLength
          # 
          deltaT = (self.tMax - self.tMin) / (self.arcLengthDivisions-1)
          tA = self.tMin + (index-1) * deltaT
          t = tA + percent * deltaT
          return t
  
      # calculates a point according to original parameterization
      def getPointAt(self, t):
          return self.f(t)
                  
      # calculates a point according to arclength parameterization
      # u - percentage of total curve distance traversed, in [0,1]
      def getPoint(self, u):
          t = self.convert_u_to_t(u)
          return self.getPointAt(t)
          
      # calculates a list of equally spaced points along the curve
      # size of list specified by divisions parameter in constructor
      def getPoints(self):
          points = []
          deltaU = 1 / (self.divisions - 1)
          for index in range(self.divisions):
              u = index * deltaU
              point = self.getPoint(u)
              points.append( point )        
          return points
  
      # calculates a tangent vector according to original parameterization
      def getTangentAt(self, t):
          h = 0.00001
          tangent = np.divide( np.subtract(self.getPointAt(t+h), self.getPointAt(t)), h )
          tangent = np.divide( tangent, np.linalg.norm(tangent) )
          return tangent
      
      # calculates a tangent vector according to arclength parameterization
      # u - percentage of total curve distance traversed, in [0,1]
      def getTangent(self, u):
          t = self.convert_u_to_t(u)
          return self.getTangentAt(t)
      
      
      # returns { "tangents" : list, "normals" : list, "binormals" : list }
      #  length of each list = divisions
      def getFrames(self):
          tangents  = []
          normals   = []
          binormals = []
          
          # calculate tangents ------------------------
          deltaU = 1 / (self.divisions - 1)
          for index in range( self.divisions ):
              u = index * deltaU
              tangent = self.getTangent(u)
              tangents.append( tangent )
  
          # calculate normals ------------------------
          # reference: https://cs.indiana.edu/ftp/techreports/TR425.pdf
          
          # calculate initial normal vector
          arbitraryVector = [1,1,1] # if tangents[0] != c*[1,1,1] else [1,1,-1]
          normal = np.cross(tangents[0], arbitraryVector)
          normals.append(normal)
          
          # calculate the remaining normal vectors
          for n in range(1, self.divisions):
          
              # default normal vector equals the previous normal vector
              # (will be used if the tangent vector has not changed direction)
              normal =  normals[n-1]
              
              # if two vectors are parallel, 
              #  then their cross product is the zero vector (which has magnitude zero)
              crossProduct = np.cross(tangents[n-1], tangents[n])
              # therefore, if this magnitude is nonzero, 
              #  the tangent vectors have changed direction
              #  and we need to apply the same change to the normal vector.
              magnitude = np.linalg.norm(crossProduct)
              if (magnitude > 0.0001):
                  # normalize the cross product vector
                  crossProduct = np.divide(crossProduct, magnitude)
                  
                  # calculate the angle between the previous and current tangent vector
                  theta = math.acos( np.dot(tangents[n-1], tangents[n]) )
  
                  # calculate the transformation the rotates previous tangent vector
                  #  into the current tangent vector
                  rotationMatrix = MatrixFactory.makeRotationAxisAngle(crossProduct,theta)
                  
                  # apply this transformation to the normal vector
                  normal = list( rotationMatrix[0:3,0:3] @ normal )
                  
              # assign next normal vector
              normals.append(normal)
          
          # calculate binormals ------------------------
          for index in range( self.divisions ):
              binormal = np.cross( tangents[index], normals[index] )
              binormals.append( binormal )
  
          # calculations complete ------------------------
          return { "tangents" : tangents, "normals" : normals, "binormals" : binormals }
  
      