import math


def torsion_angle(c1, c2, c3, c4):
    """
    float <- torsion_angle(a, b, c, d)
    returns the torsion angle in degrees between 3D pts a,b,c,d
    """

    v1 = (c1[0]-c2[0], c1[1]-c2[1], c1[2]-c2[2])
    v2 = (c2[0]-c3[0], c2[1]-c3[1], c2[2]-c3[2])
    v3 = (c3[0]-c4[0], c3[1]-c4[1], c3[2]-c4[2])

    p = (v2[1]*v1[2] - v1[1]*v2[2],
         v1[0]*v2[2] - v2[0]*v1[2],
         v2[0]*v1[1] - v1[0]*v2[1])

    q = (v3[1]*v2[2] - v2[1]*v3[2],
         v2[0]*v3[2] - v3[0]*v2[2],
         v3[0]*v2[1] - v2[0]*v3[1])

    n = 1.0 / math.sqrt( p[0]*p[0] + p[1]*p[1] + p[2]*p[2] )
    p = (p[0]*n, p[1]*n, p[2]*n )
    
    n = 1.0 / math.sqrt( q[0]*q[0] + q[1]*q[1] + q[2]*q[2] )
    q = (q[0]*n, q[1]*n, q[2]*n )

    xtheta = p[0]*q[0] + p[1]*q[1] + p[2]*q[2]
  
    if xtheta > 1.0: xtheta = 1.0
    if xtheta < -1.0: xtheta = -1.0
    theta = math.acos(xtheta) * 57.29578
    absth = math.fabs(theta)

    if absth < 0.001:
        return 0.0
    elif math.fabs(absth - 180.0) < 0.001:
        return 180.0

    s = v1[0]*q[0] + v1[1]*q[1] + v1[2]*q[2]
    if s < 0.0:
        theta = 360.0 - theta
  
    if theta > 180.0:
        theta = theta - 360.0

    return theta