import math


def calcDSH(q1, e1, i1, O1, w1, q2, e2, i2, O2, w2):
    """ Calculate the Southworth and Hawking meteoroid orbit dissimilarity criterion.

    Arguments:
        q1: [double] perihelion distance of the first orbit
        e1: [double] num. eccentricity of the first orbit
        i1: [double] inclination of the first orbit (rad)
        O1: [double] longitude of ascending node of the first orbit (rad)
        w1: [double] argument of perihelion of the first orbit (rad)
        q2: [double] perihelion distance of the second orbit
        e2: [double] num. eccentricity of the second orbit
        i2: [double] inclination of the second orbit (rad)
        O2: [double] longitude of ascending node of the second orbit (rad)
        w2: [double] argument of perihelion of the second orbit (rad)

    Return:
        [double] D_SH value

    """

    rho = 1

    if (abs(O2 - O1) > math.pi):
        rho = -1


    I21 = math.acos(math.cos(i1)*math.cos(i2) + math.sin(i1)*math.sin(i2)*math.cos(O2 - O1))


    asin_val = math.cos((i2 + i1)/2.0)*math.sin((O2 - O1)/2.0)*(1/math.cos(I21/2.0))

    # Name sure the value going into asin is not beyond the bounds due to numerical reasons
    if abs(asin_val) > 1:
        asin_val = math.copysign(1.0, asin_val)

    pi21 = w2 - w1 + 2*rho*math.asin(asin_val)

    DSH2 = pow((e2 - e1), 2) + pow((q2 - q1), 2) + pow((2 * math.sin(I21/2.0)), 2) + \
        pow((e2 + e1)/2.0, 2)*pow((2 * math.sin(pi21 / 2.0)), 2)


    return math.sqrt(DSH2)