# Curve curvature in numpy

EDIT: I put together this answer off and on over a couple of hours, so I missed your latest edits indicating that you only needed curvature. Hopefully, this answer will be helpful regardless.

Other than doing some curve-fitting, our method of approximating derivatives is via finite differences. Thankfully, `numpy` has a `gradient` method that does these difference calculations for us, taking care of the details of averaging previous and next slopes for each interior point and leaving each endpoint alone, etc.

```import numpy as np

a = np.array([ [  0.  ,   0.  ],[  0.3 ,   0.  ],[  1.25,  -0.1 ],
[  2.1 ,  -0.9 ],[  2.85,  -2.3 ],[  3.8 ,  -3.95],
[  5.  ,  -5.75],[  6.4 ,  -7.8 ],[  8.05,  -9.9 ],
[  9.9 , -11.6 ],[ 12.05, -12.85],[ 14.25, -13.7 ],
[ 16.5 , -13.8 ],[ 19.25, -13.35],[ 21.3 , -12.2 ],
[ 22.8 , -10.5 ],[ 23.55,  -8.15],[ 22.95,  -6.1 ],
[ 21.35,  -3.95],[ 19.1 ,  -1.9 ]])
```

Now, we compute the derivatives of each variable and put them together (for some reason, if we just call `np.gradient(a)`, we get a list of arrays…not sure off the top of my head what’s going on there, but I’ll just work around it for now):

```dx_dt = np.gradient(a[:, 0])
velocity = np.array([ [dx_dt[i], dy_dt[i]] for i in range(dx_dt.size)])
```

This gives us the following vector for `velocity`:

```array([[ 0.3  ,  0.   ],
[ 0.625, -0.05 ],
[ 0.9  , -0.45 ],
[ 0.8  , -1.1  ],
[ 0.85 , -1.525],
[ 1.075, -1.725],
[ 1.3  , -1.925],
[ 1.525, -2.075],
[ 1.75 , -1.9  ],
[ 2.   , -1.475],
[ 2.175, -1.05 ],
[ 2.225, -0.475],
[ 2.5  ,  0.175],
[ 2.4  ,  0.8  ],
[ 1.775,  1.425],
[ 1.125,  2.025],
[ 0.075,  2.2  ],
[-1.1  ,  2.1  ],
[-1.925,  2.1  ],
[-2.25 ,  2.05 ]])
```

which makes sense when glancing at the scatterplot of `a`.

Now, for speed, we take the length of the velocity vector. However, there’s one thing that we haven’t really kept in mind here: everything is a function of `t`. Thus, `ds/dt` is really a scalar function of `t` (as opposed to a vector function of `t`), just like `dx/dt` and `dy/dt`. Thus, we will represent `ds_dt` as a `numpy` array of values at each of the one second time intervals, each value corresponding to an approximation of the speed at each second:

```ds_dt = np.sqrt(dx_dt * dx_dt + dy_dt * dy_dt)
```

This yields the following array:

```array([ 0.3       ,  0.62699681,  1.00623059,  1.36014705,  1.74588803,
2.03254766,  2.32284847,  2.57512136,  2.58311827,  2.48508048,
2.41518633,  2.27513736,  2.50611752,  2.52982213,  2.27623593,
2.31651678,  2.20127804,  2.37065392,  2.8487936 ,  3.04384625])
```

which, again, makes some sense as you look at the gaps between the dots on the scatterplot of `a`: the object picks up speed, slowing down a bit as it takes the corner, and then speeds back up even more.

Now, in order to find the unit tangent vector, we need to make a small transformation to `ds_dt` so that its size is the same as that of `velocity` (this effectively allows us to divide the vector-valued function `velocity` by the (representation of) the scalar function `ds_dt`):

```tangent = np.array([1/ds_dt] * 2).transpose() * velocity
```

This yields the following `numpy` array:

```array([[ 1.        ,  0.        ],
[ 0.99681528, -0.07974522],
[ 0.89442719, -0.4472136 ],
[ 0.5881717 , -0.80873608],
[ 0.48685826, -0.87348099],
[ 0.52889289, -0.84868859],
[ 0.55965769, -0.82872388],
[ 0.5922051 , -0.80578727],
[ 0.67747575, -0.73554511],
[ 0.80480291, -0.59354215],
[ 0.90055164, -0.43474907],
[ 0.97796293, -0.2087786 ],
[ 0.99755897,  0.06982913],
[ 0.9486833 ,  0.31622777],
[ 0.77979614,  0.62603352],
[ 0.48564293,  0.87415728],
[ 0.03407112,  0.99941941],
[-0.46400699,  0.88583154],
[-0.67572463,  0.73715414],
[-0.73919634,  0.67349   ]])
```

Note two things: 1. At each value of `t``tangent` is pointing in the same direction as `velocity`, and 2. At each value of `t``tangent` is a unit vector. Indeed:

In [12]:

```In [12]: np.sqrt(tangent[:,0] * tangent[:,0] + tangent[:,1] * tangent[:,1])
Out[12]:
array([ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,
1.,  1.,  1.,  1.,  1.,  1.,  1.])
```

Now, since we take the derivative of the tangent vector and divide by its length to get the unit normal vector, we do the same trick (isolating the components of `tangent` for convenience):

```tangent_x = tangent[:, 0]
tangent_y = tangent[:, 1]

dT_dt = np.array([ [deriv_tangent_x[i], deriv_tangent_y[i]] for i in range(deriv_tangent_x.size)])

length_dT_dt = np.sqrt(deriv_tangent_x * deriv_tangent_x + deriv_tangent_y * deriv_tangent_y)

normal = np.array([1/length_dT_dt] * 2).transpose() * dT_dt
```

This gives us the following vector for `normal`:

```array([[-0.03990439, -0.9992035 ],
[-0.22975292, -0.97324899],
[-0.48897562, -0.87229745],
[-0.69107645, -0.72278167],
[-0.8292422 , -0.55888941],
[ 0.85188045,  0.52373629],
[ 0.8278434 ,  0.56095927],
[ 0.78434982,  0.62031876],
[ 0.70769355,  0.70651953],
[ 0.59568265,  0.80321988],
[ 0.41039706,  0.91190693],
[ 0.18879684,  0.98201617],
[-0.05568352,  0.99844847],
[-0.36457012,  0.93117594],
[-0.63863584,  0.76950911],
[-0.89417603,  0.44771557],
[-0.99992445,  0.0122923 ],
[-0.93801622, -0.34659137],
[-0.79170904, -0.61089835],
[-0.70603568, -0.70817626]])
```

Note that the normal vector represents the direction in which the curve is turning. The vector above then makes sense when viewed in conjunction with the scatterplot for `a`. In particular, we go from turning down to turning up after the fifth point, and we start turning to the left (with respect to the x axis) after the 12th point.

Finally, to get the tangential and normal components of acceleration, we need the second derivatives of `s``x`, and `y` with respect to `t`, and then we can get the curvature and the rest of our components (keeping in mind that they are all scalar functions of `t`):

```d2s_dt2 = np.gradient(ds_dt)