How do I extract a curve in xyz from curvature/torsion using Frenet Serret equations?

Question

Somewhat of a continuation of this, opening this scab because I have the same question but the solution was not covered in post and I have been banging my head against the wall for a week trying to figure this out. I appreciate any and all help!

Same basic issue; trying to solve Frenet-Serret equations to find the curve(xyz) given curvature and torsion, instead of taking the curve (xyz) and deriving curvature and torsions.

Let $\bf T$, $\bf N$ and $\bf B$ be the Frenet frame vectors. And let us define the matrix $\bf Q$ to be

$\bf Q = \left [\begin{array}{c}\bf T \\\bf N \\\bf B\end{array}\right ] $

The Frenet-Serret equations can be defined as:

$\left ( \frac{d}{ds} \right ) {\bf Q} = {\bf Q}' = {\bf F} \cdot {\bf Q}$

where

${\bf F}= \left [ \begin{array}{} 0 & \kappa & 0 \\ -\kappa & 0 & \tau \\ 0 & -\tau & 0 \end{array} \right ]$

I have the curvature/torsion for the elliptical helix defined below;

s = np.linspace(0, 20, 1000)

# Elliptical Helix
mod = 0.5
x = mod * np.cos(s)
y = np.sin(s)
z = s / (2 * np.pi)

xyz = np.vstack((x, y, z)).T
dxyz = np.apply_along_axis(np.gradient, axis=0, arr=xyz)
ddxyz = np.apply_along_axis(np.gradient, axis=0, arr=dxyz)
dS = np.sqrt(np.sum(dxyz ** 2, axis=1))

T = np.divide(dxyz, magn(dxyz, 3))
dT = np.gradient(T)[0]
N = np.divide(dT, magn(dT, 3))
B = np.cross(T, N)
k = magn(np.cross(dxyz, ddxyz), 1)/(magn(dxyz, 1)**3)
t = np.sum(-B*N, axis=1)

where (from dipy):

def magn(xyz, n=1):
    ''' magnitude of vector
    '''
    mag = np.sum(xyz**2, axis=1)**0.5
    imag = np.where(mag == 0)
    mag[imag] = np.finfo(float).eps

    if n > 1:
        return np.tile(mag, (n, 1)).T
    return mag.reshape(len(mag), 1)

I know the reduction from xyz to $\bf T$, $\bf N$, $\bf B$, $\bf k$ and $\bf t$ was done correctly because when I discretely sum over $\bf T_{0-i}$ and multiply by $\bf \Delta s_{0-i}$ I get xyz and xyz from $\int{ \bf{TdS}}$ superimposed;

xyz_exact = []
for i in range(len(xyz)):
    xyz_exact.append(np.sum(T[0:i] * dS[0:i, None], axis=0) + np.array([1, 0, 0]))
xyz_exact = np.array(xyz_exact)

x_exact = xyz_exact.T[0]
y_exact = xyz_exact.T[1]
z_exact = xyz_exact.T[2]


fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(x, y, z)
ax.plot(x_exact, y_exact, z_exact)
plt.show()

I now want to be able to do the same but summing over a $\bf T$ obtained from curvature $\bf k$ and torsion $\bf t$ defined above.

1. I assemble the initial ${\bf Q}_i = \left [ \begin{array}{c} {\bf T}_i \\ {\bf N}_i \\ {\bf B}_i \end{array} \right ]$ from the calculated $\bf T_0$, $\bf N_0$, $\bf B_0$

   T_0 = T[0]
   N_0 = N[0]
   B_0 = B[0]
   

2. Iteratively find the derivatives by means of ${\bf Q}_i' = {\bf F}(s) {\bf Q}_i$ in the following for loop also update the quantities by means of ${\bf Q}_i= {\bf Q}_{i-1} + {\bf Q}_{i-1}'$;

   N_calc = []
   T_calc = []
   B_calc = []
   
   N_calc.append(N_0)
   B_calc.append(B_0)
   T_calc.append(T_0)
   
   k_calc = k.T[0]
   t_calc = t
   
   T_d = []
   N_d = []
   B_d = []
   
   for i in range(len(s)):
       T_step = k_calc[i] * N_calc[i]
       N_step = (-k_calc[i] * T_calc[i] + t_calc[i] * B_calc[i])
       B_step = -t_calc[i] * N_calc[i]
   
       T_d.append(T_step)
       N_d.append(N_step)
       B_d.append(B_step)
   
       T_calc.append(mag(T_calc[-1] + T_d[-1]))
       N_calc.append(mag(N_calc[-1] + N_d[-1]))
       B_calc.append(mag(B_calc[-1] + B_d[-1]))
   

3. And attempt to plot the resulting xyz by summing over the calculated $\bf T$;

   xyz_calc = []
   for i in range(len(xyz)):
       xyz_calc.append(np.sum(T_calc[0:i] * dS[0:i, None], axis=0) + np.array([0.5, 0, 0]))
   xyz_calc = np.array(xyz_calc)
   
   x_calc = xyz_calc.T[0]
   y_calc = xyz_calc.T[1]
   z_calc = xyz_calc.T[2]
   
   fig = plt.figure()
   ax = fig.gca(projection='3d')
   ax.plot(x, y, z)
   ax.plot(x_calc, y_calc, z_calc)
   
   plt.show()
   

By using this (flawed) approach, my reconstruction appears here in green and another perspective. Its clear that the reconstruction starts ok but does not grow. Closer examination of the z-component of the calculated $\bf T $ and the exact $\bf T$ from above show why, whereas the x and y components appear to trace each other (minus the small error probably due to euclidean calculation of $\bf dS$).

And so I have no idea why. Would love any and all help on this. For reference, here is the whole code;

import numpy as np
import matplotlib.pyplot as plt

def magn(xyz, n=1):
    ''' magnitude of vector
    '''
    mag = np.sum(xyz**2, axis=1)**0.5
    imag = np.where(mag == 0)
    mag[imag] = np.finfo(float).eps

    if n > 1:
        return np.tile(mag, (n, 1)).T
    return mag.reshape(len(mag), 1)

s = np.linspace(0, 20, 1000)

# Elliptical Helix
mod = 0.5
x = mod * np.cos(s)
y = np.sin(s)
z = s / (2 * np.pi)

xyz = np.vstack((x, y, z)).T
dxyz = np.apply_along_axis(np.gradient, axis=0, arr=xyz)
ddxyz = np.apply_along_axis(np.gradient, axis=0, arr=dxyz)
dS = np.sqrt(np.sum(dxyz ** 2, axis=1))

T = np.divide(dxyz, magn(dxyz, 3))
dT = np.gradient(T)[0]
N = np.divide(dT, magn(dT, 3))
B = np.cross(T, N)
k = magn(np.cross(dxyz, ddxyz), 1)/(magn(dxyz, 1)**3)
t = np.sum(-B*N, axis=1)

xyz_exact = []
for i in range(len(xyz)):
    xyz_exact.append(np.sum(T[0:i] * dS[0:i, None], axis=0) + np.array([0.5, 0, 0]))
xyz_exact = np.array(xyz_exact)

x_exact = xyz_exact.T[0]
y_exact = xyz_exact.T[1]
z_exact = xyz_exact.T[2]

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(x, y, z)
ax.plot(x_exact, y_exact, z_exact)
plt.show()


s_0 = 0
k_0 = 0
t_0 = 0
T_0 = np.array([1, 0, 0])
N_0 = np.array([0, 1, 0])
B_0 = np.cross(T_0, N_0)

# Rebuild Frame
N_calc = []
T_calc = []
B_calc = []

N_calc.append(N[0])
B_calc.append(B[0])
T_calc.append(T[0])

k_calc = k.T[0]
t_calc = t

T_d = []
N_d = []
B_d = []

for i in range(len(s)):
    T_step = k_calc[i] * N_calc[i] * dS[i]
    N_step = (-k_calc[i] * T_calc[i] + t_calc[i] * B_calc[i]) * dS[i]
    B_step = -t_calc[i] * N_calc[i] * dS[i]

    T_d.append(T_step)
    N_d.append(N_step)
    B_d.append(B_step)

    T_calc.append(mag(T_calc[-1] + T_d[-1]))
    N_calc.append(mag(N_calc[-1] + N_d[-1]))
    B_calc.append(mag(B_calc[-1] + B_d[-1]))

T_calc = np.array(T_calc)
B_calc = np.array(B_calc)
N_calc = np.array(N_calc)

xyz_calc = []
for i in range(len(xyz)):
    xyz_calc.append(np.sum(T_calc[0:i] * dS[0:i, None], axis=0) + np.array([0.5, 0, 0]))
xyz_calc = np.array(xyz_calc)

x_calc = xyz_calc.T[0]
y_calc = xyz_calc.T[1]
z_calc = xyz_calc.T[2]

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(x, y, z)
ax.plot(x_calc, y_calc, z_calc)
plt.show()

plt.figure()
plt.plot(T.T[0])
plt.plot(T_calc.T[0])
plt.show()

plt.figure()
plt.plot(T.T[1])
plt.plot(T_calc.T[1])
plt.show()

plt.figure()
plt.plot(T.T[2])
plt.plot(T_calc.T[2])
plt.show()