Math for simple 3D coordinate rotation (python) I'm rotating points from one coordinate system to the other, but drawing a blank how to do this. 
I've done some reading about Euler angles but after staring at this GIF for a while I just get dizzy.

above: from here
What I am doing:
The new $z$ axis is perpendicular the plane containing the origin, $p_1$ and $p_2$, and the new $x$ axis goes through $p_1$. There is no change in origin.
Looking for $p_{calc}$ in the rotated cartesian coordinates. Here is a plot of some normals. Math answer that's simple enough for me to re-write in python is fine. 
edit: This PDF seems to be helpful.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

p1    = np.array([-1489., -4913.,  4345.])
p2    = np.array([ 2633., -3268.,  5249.])
pcalc = np.array([-3210., -4390.,  3930.])

def normit(v):
    return v / np.sqrt((v**2).sum())

n1, n2, ncalc = [normit(p) for p in [p1, p2, pcalc]]

new_zaxis  = normit(np.cross(n1, n2))
new_xaxis  = n1
zero       = np.zeros(3) 

fig = plt.figure(figsize=[10, 8])
ax  = fig.add_subplot(1, 1, 1, projection='3d')

x, y, z = zip(zero, new_xaxis)
plt.plot(x, y, z, '-k', linewidth=3)

x, y, z = zip(zero, new_zaxis)
plt.plot(x, y, z, '--k', linewidth=3)

x, y, z = zip(zero, n2)
plt.plot(x, y, z, '-r', linewidth=1)

x, y, z = zip(zero, ncalc)
plt.plot(x, y, z, '--g', linewidth=1)

ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
ax.set_zlim(-1, 1)

plt.show()

 A: note: A nicer looking and correct answer will still get accepted, thanks!
I've read on page 27 here that a 3x3 transform matrix can be just the nine dot products - thank you U. Auckland's prof. Kelly!


above x2: screenshots from here.
Here is a very ugly implementation which seems to work. 
new_yaxis = -np.cross(new_xaxis, new_zaxis)

# new axes:
nnx, nny, nnz = new_xaxis, new_yaxis, new_zaxis
# old axes:
nox, noy, noz = np.array([1, 0, 0, 0, 1, 0, 0, 0, 1], dtype=float).reshape(3, -1)

# ulgiest rotation matrix you can imagine
top = [np.dot(nnx, n) for n in [nox, noy, noz]]
mid = [np.dot(nny, n) for n in [nox, noy, noz]]
bot = [np.dot(nnz, n) for n in [nox, noy, noz]]

def newit(vec):
    xn = sum([p*q for p,q in zip(top, vec)])
    yn = sum([p*q for p,q in zip(mid, vec)])
    zn = sum([p*q for p,q in zip(bot, vec)])
    return np.hstack((xn, yn, zn))

Let's see what happens...
nnx:         array([-0.22139284, -0.73049229,  0.64603887])
newit(nnx):  array([ 1.,  0.,  0.])

nny:         array([ 0.88747002,  0.1236673 ,  0.44396325])
newit(nny):  array([ 0.,  1.,  0.])

nnz:         array([-0.40420561,  0.67163042,  0.62091095])
newit(nnz:   array([ 0.,  0.,  1.])

OK then, this seems to be the right way to go.
A: thank you uhoh for the math. here is the same using np inner
new_yaxis = -np.cross(new_xaxis, new_zaxis)


# new axes:
new_axes  = np.array([new_xaxis, new_yaxis, new_zaxis])
# old axes:
old_axes = np.array([1, 0, 0, 0, 1, 0, 0, 0, 1], dtype=float).reshape(3, -1)

rotation_matrix = np.inner(new_axes,old_axes)
#vec can be vectors Nx3
def newit(vec):
  return np.inner(vec,rotation_matrix)

And your test
new = newit(new_axes)
#new_axes:
[[-0.22139284 -0.73049229  0.64603887]
 [ 0.88747002  0.1236673   0.44396325]
 [-0.40420561  0.67163042  0.62091095]]
#new:
[[ 1.00000000e+00 -2.06024710e-17  6.10077552e-17]
 [-2.06024710e-17  1.00000000e+00  4.13972145e-17]
 [ 6.10077552e-17  4.13972145e-17  1.00000000e+00]]
```

