Implementing Vincenty's Formula in PHP

Question

I've been attempting to implement Vincenty's formulae with the following:

    /* Implemented using Vincenty's formulae from http://en.wikipedia.org/wiki/Vincenty%27s_formulae,
 * answers "Direct Problem".
 * $latlng is a ('lat'=>x1, 'lng'=>y1) array
 * $distance is in miles
 * $angle is in degrees
 */
function addDistance($latlng, $distance, $bearing) {
    //variables
    $bearing = deg2rad($bearing);
    $iterations = 20; //avoid too-early termination while avoiding the non-convergant case

    //knowns
    $f = EARTH_SPHEROID_FLATTENING; //1/298.257223563
    $a = EARTH_RADIUS_EQUATOR_MILES; //3963.185 mi
    $phi1 = deg2rad($latlng['lat']);
    $l1 = deg2rad($latlng['lng']);
    $b = (1 - $f) * $a;

    //first block
    $tanU1 = (1-$f)*tan($phi1);
    $U1 = atan($tanU1);
    $sigma1 = atan($tanU1 / cos($bearing));
    $sinalpha = cos($U1)*sin($bearing);
    $cos2alpha = (1 - $sinalpha) * (1 + $sinalpha);
    $usquared = $cos2alpha * (($a*$a - $b*$b) / 2);
    $A = 1 + ($usquared)/16384 * (4096+$usquared*(-768+$usquared*(320 - 175*$usquared)));
    $B = ($usquared / 1024)*(256*$usquared*(-128 + $usquared * (74 - 47*$usquared)));

    //the loop - determining our value
    $sigma = $distance / ($b * $A);
    for($i = 0; $i < $iterations; ++$i) {
        $twosigmam = 2*$sigma1 + $sigma;
        $delta_sigma = $B * sin($sigma) * (cos($twosigmam)+(1/4)*$B*(cos(-1 + 2*cos(cos($twosigmam))) - (1/6)*$B*cos($twosigmam)*(-3+4*sin(sin($sigma)))*(-3+4*cos(cos($twosigmam)))));
        $sigma = $distance / ($b * $A) + $delta_sigma;
    }

    //second block
    $phi2 = atan((sin($U1)*cos($sigma)+cos($U1)*sin($sigma)*cos($bearing)) / ((1-$f) * sqrt(sin($sinalpha) + pow(sin($U1)*sin($sigma) - cos($U1)*cos($sigma)*cos($bearing), 2))));
    $lambda = atan((sin($sigma) * sin($bearing)) / (cos($U1)*cos($sigma) - sin($U1)*sin($sigma)*cos($bearing)));
    $C = ($f / 16)* $cos2alpha * (4+$f*(4-3*$cos2alpha));
    $L = $lambda - (1 - $C) * $f * $sinalpha * ($sigma + $C*sin($sigma)*(cos($twosigmam)+$C*cos($sigma)*(-1+2*cos(cos($twosigmam)))));
    $alpha2 = atan($sinalpha / (-sin($U1)*sin($sigma) + cos($U1)*cos($sigma)*cos($bearing)));

    //and return our results
    return array('lat' => rad2deg($phi2), 'lng' => rad2deg($lambda));
}

    var_dump(addDistance(array('lat' => 93.129, 'lng' => -43.221), 20, 135);

The issue is that the results are not reasonable - I'm getting variances of up to 20 latitude and longitude keeping the distance at 20. Is it not in units of elliptical distance on the sphere? Am I misunderstanding something, or is my implementation flawed?

I'm aware that http://www.movable-type.co.uk/scripts/latlong.html has a simpler implementation of distance+bearing, but I'm interested in the higher accuracy of Vincenty's formula. — freeone3000, Jul 10 '12 at 15:32
I usually do this in sql (after cutting results with a simple boundary box). I'll see if I still have the php function for you. — somedev, Jul 10 '12 at 15:37

AakashM · Accepted Answer · 2012-07-11T08:36:21.283

There are a number of errors in transcription from the wikipedia page Direct Problem section:

Your u² expression has 2 in the denominator where it should have b²;
Your A and B expressions are inconsistent about whether the initial fraction factor needs to be parenthesised to correctly express a / b * c as (a/b) * c - what happens without parentheses is a php syntax issue which I don't know the answer to, but you should favour clarity;
You should be iterating "until there is no significant change in sigma", which may or may not happen in your fixed number of iterations;
There are errors in your DELTA_sigma formula:
- on the wikipedia page, the first term inside the square bracket [ is cos sigma (-1 etc, whereas you have cos (-1 etc, which is very different;
- in the same formula and also later, note that cos² x means (cos x)(cos x), not cos cos x!
Your phi_2 formula has a sin($sinalpha) where it should have a sin($sinalpha)*sin($sinalpha);

I think that's all.

Wow. Wasn't aware I messed up that badly. Fixing the transcription errors fixed the problems, thanks. "No significant change in sigma" is nice, but that proves difficult to compute due to the subjectivity of "sigma" and floating-point error, so I've settled for a fixed number of iterations as an alternative. — freeone3000, Jul 12 '12 at 09:47

score 0 · Answer 2 · answered Jul 10 '12 at 15:43

0

Have you tried this:

https://github.com/treffynnon/Geographic-Calculations-in-PHP

answered Jul 10 '12 at 15:43

somedev

1,053
1
9
27

It's interesting, but it solves the "Inverse Problem", instead of the Direct Problem. I understand the process behind the math, vaguely, it's just the implementation I fear I've messed up. – freeone3000 Jul 10 '12 at 16:38

Implementing Vincenty's Formula in PHP

2 Answers2

Linked