How to calculate the spiral around a curve? I have a curve given by a set of points and want to build a spiral around it (like this).

I tried it by
STEP 1:
calculated the vector of each step by
$$v_n = \left<x_n,y_n,z_n\right> - \left<x_{n-1},y_{n-1},z_{n-1}\right>$$
STEP 2:
I assumed a vector/plane of the rotating trajectory as
$$ w_n = \left<\cos(i),1,\sin(i)\right>$$
where $i$ represents the intervals and increases by the point number.
And calculated the perpendicular vector to find the spiral points by
$$X = v_y \cdot w_z - v_z \cdot w_y $$
$$Y = v_z \cdot w_x - v_x \cdot w_z $$
$$Z = v_x \cdot w_y - v_y \cdot w_x $$
And, of course, I normalized the scale.
The spiral is formed, but I have discontinuity or deformation at sharp angles (probably because the vector direction is changed). It is not a matter of smoothness, as the spiral changes its direction to rotate in the opposite direction.
Where did I do wrong?
 A: spiral will evolve according to helical angle $$\phi_n=\phi_{n-1}+2 \pi \frac{|\vec{x}_n-\vec{x}_{n-1}|}{\lambda}$$ where $\lambda$ is the pitch of the helix and is measured relative to an "up" orientation normal vector $\vec{u}$ which is the component of the acceleration vector $\vec{a}_n$ perpendicular to the velocity vector $\vec{v}_n$ $$\vec{u}_n=\hat{a}_n-(\hat{a}_n \cdot \hat{v}_n ) \hat{v}_n$$ where $$ \vec{v}_n=\vec{x}_n-\vec{x}_{n-1}$$ and $$\vec{a}_n=\vec{v}_n-\vec{v}_{n-1}$$ the location of the helix is then given by $$\vec{h}_n=\vec{x}_n+\hat{u}_nR\cos \phi_n+(\hat{u}_n\times\hat{v}_n)R\sin \phi_n$$ where $R$ is the radius of the helix
A: In three dimensions, you need two vector-valued functions to describe a helix:
$$\begin{array}{ll}
\vec{n}(t) & \text{Axis} \\
\vec{u}(t) & \text{Direction of } \theta = 0 \\
\end{array}$$
with the two being perpendicular at each point $(t)$,
$$\vec{n}(t) \cdot \vec{u}(t) = 0$$
These two vectors and their cross product, normalized to unit length,
$$\left\lbrace\begin{aligned}
\displaystyle \hat{n}(t) &= \frac{\vec{n}(t)}{\left\lVert\vec{n}(t)\right\rVert} \\
\displaystyle \hat{u}(t) &= \frac{\vec{u}(t)}{\left\lVert\vec{u}(t)\right\rVert} \\
\displaystyle \hat{v}(t) &= \frac{\vec{n}(t) \times \vec{u}(t)}{\left\lVert\vec{n}(t) \times \vec{u}(t)\right\rVert} \\
\end{aligned}\right.$$
form an orthonormal basis at each point $t$.
If $\lambda$ is the pitch (turns per $t$ increasing by one) and $r$ is the radius of the helix, then the helix is described by vector-valued function
$$\vec{p}(t) = \vec{n}(t) + \hat{u}(t) \bigr( r \cos( 2 \pi \lambda t) \bigr) + \hat{v}(t) \bigr( r \sin( 2 \pi \lambda t) \bigr)$$

Typically, you use the points the curve passes through, to define a cubic Bézier curve between consecutive pair of points $\vec{p}_i$ and $\vec{p}_{i+1}$, so that each curve segment is defined by a vector-valued function $\vec{c}_i$:
$$\vec{c}_i(t) = (1-t)^3 \vec{p}_{i} + 3 (1-t)^2 t \vec{a}_i + 3 (1-t) t^2 \vec{b}_i + t^3 \vec{p}_{i+1}, \quad t = 0 .. 1$$
The two additional vectors per curve segment, $\vec{a}_i$ and $\vec{b}_i$, are control points.  The curve passes only through the endpoints $\vec{p}_i$ (at $t 0 0$) and $\vec{p}_{i+1}$ (at $t = 1$), but not through the control points.
Note that if you have $n$ points, $i = 0 .. n-1$, you have $n-1$ consecutive pairs, and thus $n-1$ curve segments, $i = 0 .. n-2$.
Usually the control points are calculated based on the endpoints and on the preceding and succeeding points $\vec{p}_{i-1}$ and $\vec{p}_{i+2}$.  For example, you can try
$$\begin{aligned}
\displaystyle \vec{a}_{i} &= \frac{\vec{p}_{i-1} + \vec{p}_{i+1}}{2} \\
\displaystyle \vec{b}_{i} &= \frac{\vec{p}_{i} + \vec{p}_{i+2}}{2} \\
\end{aligned} \quad \text{except} \quad
\begin{aligned}
\displaystyle \vec{a}_0 &= \frac{\vec{p}_0 + \vec{p}_1}{2} \\
\displaystyle \vec{b}_{n-2} &= \frac{\vec{p}_{n-2} + \vec{p}_{n-1}}{2} \\
\end{aligned}$$
which is probably the simplest way one can define the additional control points, and get a G1-continuous curve.
The important thing is that the two control points around each end point are in the same line, so that the axis is G1-continuous: the curve itself being continuous (G0), but also its tangent is continuous.
In this case, the axis vector is the curve tangent, i.e.
$$\displaystyle \vec{n}_{i}(t) = \frac{d \vec{c}_i(t)}{d t} =
(-3 t^2 + 6 t - 3) \vec{p}_{i} + (9 t^2 - 12 t + 3) \vec{a}_i + (-9 t^2 + 6 t) \vec{b}_i + (3 t^2) \vec{p}_{i+1}$$
The direction where $\theta = 0$ is conceptually harder.  We can either use a fixed direction (which means that if we painted a straight line on the helix when it is at rest, we'd twist the helix so that the painted line is always in that same direction with respect to the axis it twists around), or we can try to recover it somehow from the set of points we have.
For example, if the points are such that no three consecutive points $\vec{p}_i$ are collinear, we can use the normal of the plane they define:
$$\vec{u}_i = (\vec{p}_{i+1} - \vec{p}_{i}) \times (\vec{p}_{i+2} - \vec{p}_{i})$$
and then interpolate the direction somehow, for example linearly, between points:
$$\vec{u}_i(t) = (1 - t) \vec{u}_i + t \vec{u}_{i+1}$$
In all cases, we should first calculate the tangent vector $\vec{n}$ at some point $t$, and normalize that to unit length $\hat{n}$, then use the Gram-Schmidt process to orthogonalize and normalize $\vec{u}$:
$$\begin{aligned}
\hat{n} &= \frac{\vec{n}}{\sqrt{\vec{n}\cdot\vec{n}}} \\
\vec{w} &= \vec{u} - \hat{n} ( \hat{n} \cdot \vec{u} ) \\
\hat{u} &= \frac{\vec{w}}{\sqrt{\vec{w}\cdot\vec{w}}} \\
\hat{v} &= \hat{n} \times \hat{u} \\
\end{aligned}$$
where $\vec{w}$ is just a temporary vector, helping with the computation.
You do this for every point you calculate for the helix, separately.  (Alternatively, you define the vector-valued functions $\vec{n}_i(t)$ and $\vec{u}_i(t)$ so that they are perpendicular at all points $t$, but except for some very rare cases where you define the curve using sine and cosine functions from the get go, that is not worth the effort at all.)

Using cubic curves like Bézier curves above has the downside that the distance between the consecutive end points defines the local "stretch" of the helix.
In fact, it would be best if the axis was described using a curve where the curve parameter $t$ would be a product of the distance and the local compression, i.e. "helical distance" as measured in turns, and there would be an associated basis, three vectors, one of which is parallel to the tangent of the curve, at each point along the curve.  One of these is the triple clothoid, but I haven't experimented with it at all.
In practice, especially if you think about how a human could manipulate such a helix, it would be best if each point the helix axis passes through, was associated with the $\theta = 0$ direction vector also.  Then, both the axis itself and the direction vectors could be interpolated as cubic Bèziers, and orthonormalized using Gram-Schmidt at each point, with points where the two are parallel ignored and replaced with "straight" line segments (i.e. previous or next valid vector).
A: Probably your steps are too large.  You are rotating a full radian between steps.  If you use finer steps in both directions you will improve things.  If you just divide $i$ by $10$ in the $w$ equation the pitch will increase by a factor $10$.  Maybe you should interpolate $10$ steps between your $v$ points, then divide $i$ by $10$ in $w$.  The most challenging points are those where the axis of the helix bends most.  You need the steps in your helix rather small by comparison with that scale.
