### Author Topic: Dead reckoning function... (geodesic measurements)  (Read 6203 times)

#### Planetary Researcher

• Ra (Power Member)
• Posts: 134
##### Dead reckoning function... (geodesic measurements)
« on: March 05, 2005, 04:11:06 PM »
Does anyone out there have a mathematical function which will take:

initial latitude (or co-latitude)
initial longitude
direction of travel (azimuthal angle east from north)
distance of travel
radius of the sphere whose surface you're on

and output the new latitude and longitude you'd find yourself at? Needs to be accurate for small distances, relative to the radius of the sphere.

#### thare

• GIS Support Team
• Isis (Extreme Power Member)
• Posts: 1217
##### Dead reckoning function... (geodesic measurements)
« Reply #1 on: March 05, 2005, 04:12:57 PM »
It appears you are searching for some of Vincenty's formulas. Searching the web, I found a version written in Python. See the next post. Othewise, I have done this done this within ArcView 3.x calling the ReturnPointAlongGeodesic function. In ArcView 3.x and its' Avenue programming language, you can return the resultant X, Y point basically in one call:

aEndPoint = aSpheroid.ReturnPointAlongGeodesic(aStartPoint, aDistance, anAzimuth)

This states, given defined Spheroid call the routine ReturnPointAlongGeodesic sending it the Start Point (X, Y), the Distance, and the Azimuth and return the End Point. Gotta love that easy to read syntax. To bad it is basically a dead language. My extension Geodetic Profiler has some sample code you can use which calls this routine:
http://webgis.wr.usgs.gov/pgis_discussion/viewtopic.php?t=49

This can also be done in ArcMap using the geodcoord function. Search for that routine at http://arcobjectsonline.esri.com/ for some sample code that calls that routine.

There is also a downloadable Excel spreadsheet and online forms to calculate Vincenty Direct and Vincenty Inverse formulas:
http://www.ga.gov.au/nmd/geodesy/datums/calcs.jsp

Trent

#### thare

• GIS Support Team
• Isis (Extreme Power Member)
• Posts: 1217
##### Dead reckoning function... (geodesic measurements)
« Reply #2 on: March 05, 2005, 04:13:46 PM »
The original python code for Vincenty's direct formula was found at (comments removed to fit here):
http://wegener.mechanik.tu-darmstadt.de/GMT-Help/Archiv/att-8710/Geodetic_py
Code: [Select]
`# Vincenty's Direct formulae def  vinc_pt( f, a, phi1, lembda1, alpha12, s ) :    """    Returns the lat and long of projected point and reverse azimuth    given a reference point and a distance and azimuth to project.    lats, longs and azimuths are passed in decimal degrees    Returns ( phi2,  lambda2,  alpha21 ) as a tuple    """    piD4 = math.atan( 1.0 )    two_pi = piD4 * 8.0    phi1    = phi1    * piD4 / 45.0    lembda1 = lembda1 * piD4 / 45.0    alpha12 = alpha12 * piD4 / 45.0    if ( alpha12 < 0.0 ) :       alpha12 = alpha12 + two_pi    if ( alpha12 > two_pi ) :       alpha12 = alpha12 - two_pi    b = a * (1.0 - f)    TanU1 = (1-f) * math.tan(phi1)    U1 = math.atan( TanU1 )    sigma1 = math.atan2( TanU1, math.cos(alpha12) )    Sinalpha = math.cos(U1) * math.sin(alpha12)    cosalpha_sq = 1.0 - Sinalpha * Sinalpha        u2 = cosalpha_sq * (a * a - b * b ) / (b * b)    A = 1.0 + (u2 / 16384) * (4096 + u2 * (-768 + u2 * \       (320 - 175 * u2) ) )    B = (u2 / 1024) * (256 + u2 * (-128 + u2 * (74 - 47 * u2) ) )        # Starting with the approx    sigma = (s / (b * A))    last_sigma = 2.0 * sigma + 2.0   # something impossible        # Iterate the following 3 eqs unitl no sig change in sigma    # two_sigma_m , delta_sigma    while ( abs( (last_sigma - sigma) / sigma) > 1.0e-9 ) :       two_sigma_m = 2 * sigma1 + sigma       delta_sigma = B * math.sin(sigma) * ( math.cos(two_sigma_m) \             + (B/4) * (math.cos(sigma) * \             (-1 + 2 * math.pow( math.cos(two_sigma_m), 2 ) -  \             (B/6) * math.cos(two_sigma_m) * \             (-3 + 4 * math.pow(math.sin(sigma), 2 )) *  \             (-3 + 4 * math.pow( math.cos (two_sigma_m), 2 ))))) \             last_sigma = sigma       sigma = (s / (b * A)) + delta_sigma        phi2 = math.atan2 ( (math.sin(U1) * math.cos(sigma) + math.cos(U1) * math.sin(sigma) * math.cos(alpha12) ), \       ((1-f) * math.sqrt( math.pow(Sinalpha, 2) +  \       pow(math.sin(U1) * math.sin(sigma) - math.cos(U1) * math.cos(sigma) * math.cos(alpha12), 2))))        lembda = math.atan2( (math.sin(sigma) * math.sin(alpha12 )), (math.cos(U1) * math.cos(sigma) -  \            math.sin(U1) *  math.sin(sigma) * math.cos(alpha12)))        C = (f/16) * cosalpha_sq * (4 + f * (4 - 3 * cosalpha_sq ))    omega = lembda - (1-C) * f * Sinalpha *  \       (sigma + C * math.sin(sigma) * (math.cos(two_sigma_m) + \       C * math.cos(sigma) * (-1 + 2 * math.pow(math.cos(two_sigma_m),2) )))    lembda2 = lembda1 + omega    alpha21 = math.atan2 ( Sinalpha, (-math.sin(U1) * math.sin(sigma) +  \       math.cos(U1) * math.cos(sigma) * math.cos(alpha12)))    alpha21 = alpha21 + two_pi / 2.0    if ( alpha21 < 0.0 ) :       alpha21 = alpha21 + two_pi    if ( alpha21 > two_pi ) :       alpha21 = alpha21 - two_pi    phi2       = phi2       * 45.0 / piD4    lembda2    = lembda2    * 45.0 / piD4    alpha21    = alpha21    * 45.0 / piD4    return phi2,  lembda2,  alpha21    # END`