title: Implementing Vincenty's Formula in PHP
body: <p>I've been attempting to implement <a href="http://en.wikipedia.org/wiki/Vincenty%27s_formulae" rel="nofollow">Vincenty's formulae</a> with the following:</p>
  
  <pre><code>    /* Implemented using Vincenty's formulae from http://en.wikipedia.org/wiki/Vincenty%27s_formulae,
   * answers "Direct Problem".
   * $latlng is a ('lat'=&gt;x1, 'lng'=&gt;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 &lt; $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' =&gt; rad2deg($phi2), 'lng' =&gt; rad2deg($lambda));
  }
  
      var_dump(addDistance(array('lat' =&gt; 93.129, 'lng' =&gt; -43.221), 20, 135);
  </code></pre>
  
  <p>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? </p>