Recover Lat/Log from only two sets of (d,lat,log)

recover_lat_long_2.py

import numpy
import scipy
import scipy.optimize


def latlon_system(d1, my_phi1, my_lambda1, d2, my_phi2, my_lambda2):
    """
    Non-linear objective function to solve for lat,long of unknown, given
    known positions and distances...

    Relate distance to lat,long from the Azimuthal equidistant projection
    :param d1: distance Trial one (km)
    :param my_phi1: latitude Trial one (rad)
    :param my_lambda1: longitude Trial one (rad)
    :param d2: distance Trial two (km)
    :param my_phi2: latitude trial two (rad)
    :param my_lambda2: latitude trial two (rad)
    :return:
    """

    def F(x):
        phi1 = x[0]
        lambda0 = x[1]
        t1 = numpy.sin(phi1)
        t2 = numpy.sin(my_phi1)
        t3 = numpy.cos(phi1)
        t4 = numpy.cos(my_phi1)
        t5 = numpy.cos(-my_lambda1 + lambda0)
        t6 = numpy.arccos(t1 * t2 + t3 * t4 * t5)
        t7 = numpy.sin(my_phi2)
        t8 = numpy.cos(my_phi2)
        t9 = numpy.cos(-my_lambda2 + lambda0)
        t10 = numpy.arccos(t1 * t7 + t3 * t8 * t9)
        return 6371.0 * numpy.asarray([t6, t10]) - numpy.asarray([d1, d2])

    def J(x):
        phi1 = x[0]
        lambda0 = x[1]

        t1 = numpy.cos(phi1)
        t2 = numpy.sin(my_phi1)
        t4 = numpy.sin(phi1)
        t5 = numpy.cos(my_phi1)
        t7 = -my_lambda1 + lambda0
        t8 = numpy.cos(t7)
        t12 = t5 * t1
        t15 = (t2 * t4 + t8 * t12) ** 2
        t17 = numpy.sqrt(0.1e1 - t15)
        t18 = 0.1e1 / t17
        t20 = numpy.sin(my_phi2)
        t22 = numpy.cos(my_phi2)
        t24 = -my_lambda2 + lambda0
        t25 = numpy.cos(t24)
        t29 = t22 * t1
        t32 = (t20 * t4 + t25 * t29) ** 2
        t34 = numpy.sqrt(0.1e1 - t32)
        t35 = 0.1e1 / t34
        t37 = numpy.sin(t7)
        t40 = numpy.sin(t24)
        cg1 = numpy.zeros((2, 2))
        cg1[0, 0] = -t18 * (t2 * t1 - t8 * t5 * t4)
        cg1[1, 0] = -t35 * (t20 * t1 - t25 * t22 * t4)
        cg1[0, 1] = t18 * t37 * t12
        cg1[1, 1] = t35 * t40 * t29
        return 6371.0 * cg1

    return F, J


objective, J = latlon_system(7208.479564, -.6152005524, -1.13595974, 8898.539059, -.152005524, -0.13595974)
x = scipy.optimize.root(objective, [1, 1], jac=J)
print(x)