minimum_curvature_inner Is Unstable

You, again are using the inverse cosine to calculate the angle between vectors:

dogleg = np.arccos(cos_inc - (sin_inc * cos_azi))

Use the angle code in the last issue to calculate.

A numerically stable way to calculate minimum curvature is this:

def arc2chord(t1, t2, arclen):
    """
    Calculates the relative vector between two survey stations,
    given tangents at the ends of the arc and the arc length between
    the tangents.
    Assumes survey values are correct
    and if arrays of values, the arrays must all be the same length.
    """
    is_scalar = np.ndim(arclen) == 0
    arclen = np.atleast_1d(arclen)
    cnt = len(arclen)
    t1 = np.atleast_2d(t1)
    t2 = np.atleast_2d(t2)

    t_add = t1 + t2 # add the tangent vectors; a vector that points to the end of the arc from the start
    lsqrd_t_add = np.einsum('ij,ij->i', t_add, t_add) # the length squared of the vector sum... same as: np.sum(t12*t12, axis=1)
    anti_parallel = lsqrd_t_add == 0 # test for anti-parallel tangent vectors, the singuar case
    lsqrd_t_add[anti_parallel] = 1.0 # set so we prevents div-by-zero when unitizing the direction vector
    len_t_add = np.sqrt(lsqrd_t_add) # the length of the addition vector
    norm_t_add = np.divide(t_add, len_t_add[:,None]) # normalize the addition vector to unit vector pointing to the end of the arc

    t_sub = t2 - t1 # subtract the tangent vectors; the chord on a unit circle
    lsqrd_t_sub = np.einsum('ij,ij->i', t_sub, t_sub) # the length squared of the vector subtraction
    len_t_sub = np.sqrt(lsqrd_t_sub) # the length of the subtraction vector

    alpha = 2.0 * np.arctan(np.divide(len_t_sub, len_t_add)) # the unoriented angle between the tangent vectors; the arc length on a unit circle

    geom_test = len_t_sub < alpha # do the degenerate circle geometry test, the straight hole test
    arc_2_chord = np.ones(cnt) # if degenerte, we are at unity
    arc_2_chord[geom_test] = np.divide(len_t_sub[geom_test], alpha[geom_test]) # where not unity, calc ratio

    relative_pos = (arclen * arc_2_chord)[:,None] * norm_t_add # arc-2-chord. For robust numeric evaluation the order of operations here are important
    relative_pos[anti_parallel] = np.array([np.nan, np.nan, np.nan]) # set any singuar cases to nan because they have vanished

    return (relative_pos[0], alpha[0]) if is_scalar else (relative_pos, alpha)

Zabamund / wellpathpy

minimum_curvature_inner Is Unstable #37