Algorithm from Craig 1957.

    From the original Craig paper, we can set up a pair of equations
    and solve for d13C and d18O simultaneously:

        d45 * r45 = r13 * d13
                  + 0.5 * r17 * d18
        d46 = r13 * ((r17**2 + r17 - r18) / a) * d13
            + 1 - 0.5 * r17 * ((r13**2 + r13 - r18) / a) * d18
        where a = r18 + r13 * r17 and b = 1 + r13 + r17

def delta13c_craig(r45sam, r46sam, d13cstd, r45std, r46std,
                   ks='Craig', d18ostd=23.5):
    """
    Algorithm from Craig 1957.

    From the original Craig paper, we can set up a pair of equations
    and solve for d13C and d18O simultaneously:

        d45 * r45 = r13 * d13
                  + 0.5 * r17 * d18
        d46 = r13 * ((r17**2 + r17 - r18) / a) * d13
            + 1 - 0.5 * r17 * ((r13**2 + r13 - r18) / a) * d18
        where a = r18 + r13 * r17 and b = 1 + r13 + r17
    """
    # the constants for the calculations
    # originally r13, r17, r18 = 1123.72e-5, 759.9e-6, 415.8e-5
    k = delta13c_constants()[ks]

    # TODO: not clear why need to multiply by 2?
    r13, r18 = k['S13'], 2 * k['S18']
    r17 = 2 * (k['K'] * k['S18'] ** k['A'])
    a = (r18 + r13 * r17) * (1. + r13 + r17)

    # the coefficients for the calculations
    eqn_mat = np.array([[r13, 0.5 * r17],
                        [r13 * ((r17 ** 2 + r17 - r18) / a),
                         1 - 0.5 * r17 * ((r13 ** 2 + r13 - r18) / a)]])

    # precalculate the d45 and d46 of the standard versus PDB
    r45d45std = (eqn_mat[0, 0] * d13cstd + eqn_mat[0, 1] * d18ostd)
    d46std = eqn_mat[1, 0] * d13cstd + eqn_mat[1, 1] * d18ostd

    # calculate the d45 and d46 of our sample versus PDB
    # in r45d45, r45 of PDB = r13 + r17 of PDB
    r45d45 = 1000. * (r45sam / r45std - 1.) * \
        (r13 + r17 + 0.001 * r45d45std) + r45d45std
    d46 = 1000. * (r46sam / r46std - 1.) * (1. + 0.001 * d46std) + d46std

    # solve the system of equations
    x = np.linalg.solve(eqn_mat, np.array([r45d45, d46]))
    return x[0]