# How to calculate the geodesic curvature of a discrete 3D curve?

I have coordinates of a set of points that form a closed loop that lies in a 3D surface. I know the equation of the surface and I can calculate it's surface normal at any point. I found that for a unit speed curve $$\alpha(t)$$ I can calculate the geodesic curvature as $$\kappa_g(t) = \alpha''(t)\cdot (n(t)\times\alpha'(t))$$

where $$n(t)$$ is the surface normal. I am confused about how to calculate the unit-speed derivatives $$\alpha'(t)$$ and $$\alpha''(t)$$. I think I should use some finite difference scheme like $$\alpha'(t) \approx \frac{\alpha(t + h) - \alpha(t - h)}{2h}$$

So my algorithm for calculating the geodesic curvature at any point $$P$$ of this curve should be as follows:

1. Identify a reference point on the curve to measure arc-length from.
2. For the point $$P$$ calculate the arc-length by adding up the Euclidean distances between successive points from the reference point to the point P. This arc length gives me $$t$$.
3. Find the position vectors of points at $$t + h$$ and $$t - h$$ by interpolation and then use the above finite difference equation to calculate $$\alpha'(t)$$.

Is my algorithm correct? I am not clear about how to get arc-length parameterization from the points on the curve.

Are there any standard algorithms (or even Python libraries) for such a calculation?