Hi,
I think what you need to do is first convert these cartesian vectors to rotation matrix and then to Euler system. It is a simple Geometry. Below is the brief code for the reference.
% Lets take first 2 points and find Spherical coordinates.
P1 = [-426159, 501913, 845131];
P2 = [-48717.2, 499318, 847679];
v = P1-P2;
% Let's define si and theta in such a way that.
v = [r*cos(si)*cos(theta), r*sin(theta), r*sin(si)*cos(theta)]
r = norm(v);
si = atan2(v(3),v(1));
theta = atan2(v(2),sqrt(v(1).^2+v(3).^2));
j = [cos(si)*cos(theta), sin(theta), sin(si)*cos(theta)];
% Correspond to j vector you can also find orthonormal vector to j
i = [sin(si), 0, -cos(si)];
k = [cos(si)*sin(theta), -cos(theta), sin(si)*sin(theta)];
% Rotation matrix;
m = [i',j',k'];
% You can use MATLAB inbuilt function to convert rotation matrix to Euler system
eul = rotm2eul(m);
This is the typical geometrical way of doing it, There could be more ways, For reference check out vrrotvec, vrrotvec2mat and cart2sph These function can also make your calculation easy.
Cheers
