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];
}