Solution I like the best is with atan2.
Here is my c code for this:
//x1, y1, z1 coordinates of point 1
//x2, y2, z2 coordinates of point 2 and so on
//b1, b2, b3 vectors described in text
//b1
b_a[0] = -(x1 - x2);
b_a[1] = -(y1 - y2);
b_a[2] = -(z1 - z2);
//b2
b_c[0] = x2 - x3;
b_c[1] = y2 - y3;
b_c[2] = z2 - z3;
//b3
c_d[0] = x4 - x3;
c_d[1] = y4 - y3;
c_d[2] = z4 - z3;
double n1[3];
double n2[3];
double m[3];
double x, y;
VectorNormalisation(b_c);
VectorNormalisation(b_a);
VectorNormalisation(c_d);
CrossProduct(b_a, b_c, n1);
CrossProduct(b_c, c_d, n2);
CrossProduct(n1, b_c, m);
x = DotProduct(n1, n2);
y = DotProduct(m, n2);
angle = 180.0 / PI * atan2(y, x);
////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////
//Functions used in the code:
double sqr(double x){
return x*x;
}
void VectorNormalisation(double *v)
{
double lenght = sqrt(sqr(v[0]) + sqr(v[1] + sqr(v[2])));
v[0] /= lenght;
v[1] /= lenght;
v[2] /= lenght;
}
double DotProduct(double *v, double *w)
{
return (v[0] * w[0] + v[1] * w[1] + v[2] * w[2]);
}
void CrossProduct(double *v, double *w, double *cross)
{
//
cross[0] = w[1] * v[2] - w[2] * v[1];
cross[1] = w[2] * v[0] - w[0] * v[2];
cross[2] = w[0] * v[1] - w[1] * v[0];
}