[转]php实现大地主题解算(正算)亦即根据已知点经纬度方向角距离求另一点经纬度

近期在做一个地图小应用用到根据已知点经纬度方向角距离求另一点经纬度，中间用到了百度地图api，却没找到实现这个的方法，百度一下知道这叫大地主题解算的正算，也找到一个vb的代码，由于我用的是php，遂将其改为php版本，经过百度地图测试发现精度还是很高的，1km的距离误差只有几米，附上代码

function rad($d){
    $Pi = pi();
    return $d*$Pi/180.0;
}
//函数里的四个参数$STARTLAT, $STARTLONG, $ANGLE1, $DISTANCE分别是起点纬度，起点经度，方向角，距离，放到代码里直接调用可
function computationx($STARTLAT, $STARTLONG, $ANGLE1, $DISTANCE){
    $Pi = pi();
    // return $Pi;
    $B1 = $STARTLAT;
    $L1 = $STARTLONG;
    $A1 = $ANGLE1;
    $S = $DISTANCE;
    $a = 6378245.0;
    $b = 6356752.3142;
    $c = $a*$a/$b;
    $alpha = ($a - $b) /$a;
    $e = sqrt($a *$a - $b *$b) /$a;
    $e2 = sqrt($a *$a - $b *$b) / $b;   
    $B1 = rad($B1);
    $L1 = rad($L1);
    $A1 = rad($A1);
    $W = sqrt(1 - $e *$e * (sin($B1) *sin($B1)));
    $V = $W * ($a / $b);
    //Dim W1 As Double
    $E1 = $e ;//第一偏心率
    //计算起点的归化纬度
    $W1 = $W; 
    //Sqr(1 - e1 * e1 * Sin(B1 ) * Sin(B1 ))
    $sinu1 = sin($B1) * sqrt(1 - $E1 * $E1) / $W1;
    $cosu1 = cos($B1) / $W1;
    //计算辅助函数值
    $sinA0 = $cosu1 * sin($A1);
    $cotq1 = $cosu1 * cos($A1);
    $sin2q1 = 2 * $cotq1 / ($cotq1*$cotq1 + 1);
    $cos2q1 = ($cotq1 *$cotq1 - 1) / ($cotq1*$cotq1 + 1);
    //计算系数AA,BB,CC及AAlpha, BBeta的值。
    $cos2A0 = 1 - $sinA0 *$sinA0;
    $e2 = sqrt($a *$a - $b *$b) / $b;
    $k2 = $e2 * $e2 * $cos2A0;
    //Dim aa, BB, CC, EE22, AAlpha, BBeta As Double
    $aa = $b * (1 + $k2 / 4 - 3 * $k2 * $k2 / 64 + 5 * $k2 * $k2 * $k2 / 256);
    $BB = $b * ($k2 / 8 - $k2 * $k2 / 32 + 15 * $k2 * $k2 * $k2 / 1024);
    $CC = $b * ($k2 * $k2 / 128 - 3 * $k2 * $k2 * $k2 / 512);
    $e2 = $E1 * $E1;
    $AAlpha = ($e2 / 2 + $e2 * $e2 / 8 + $e2 * $e2 * $e2 / 16) - ($e2 * $e2 / 16 + $e2 * $e2 * $e2 / 16) * $cos2A0 + (3 * $e2 * $e2* $e2 / 128) * $cos2A0 * $cos2A0;
    $BBeta = ($e2 * $e2 / 32 + $e2 * $e2 * $e2 / 32) * $cos2A0 - ($e2 * $e2 * $e2 / 64) * $cos2A0 * $cos2A0;
    //计算球面长度
    $q0 = ($S - ($BB + $CC * $cos2q1) * $sin2q1) / $aa;
    $sin2q1q0 = $sin2q1 * cos(2 * $q0) + $cos2q1 * sin(2 * $q0);
    $cos2q1q0 = $cos2q1 * cos(2 * $q0) - $sin2q1 * sin(2 * $q0);
    $q = $q0 + ($BB + 5 * $CC * $cos2q1q0) * $sin2q1q0 / $aa;
    //'// 计算经度差改正数
    $theta = ($AAlpha * $q + $BBeta * ($sin2q1q0 - $sin2q1)) * $sinA0;
    //计算终点大地坐标及大地方位角
    $sinu2 = $sinu1 * cos($q) + $cosu1 * cos($A1) * sin($q);
    $B2 = atan($sinu2 / (sqrt(1 - $E1 * $E1) * sqrt(1 - $sinu2 * $sinu2))) * 180 / $Pi;
    $lamuda = atan(sin($A1) * sin($q) / ($cosu1 * cos($q) - $sinu1 * sin($q) * cos($A1))) * 180 / $Pi;
    if (sin($A1) > 0){
                   if (sin($A1) * sin($q) / ($cosu1 * cos($q) - $sinu1 * sin($q) * cos($A1)) > 0){
                        $lamuda = abs($lamuda);
                    }else{
                        $lamuda = 180 - abs($lamuda);}
                }else{
                    if (sin($A1) * sin($q) / ($cosu1 * cos($q) - $sinu1 * sin($q) * cos($A1)) > 0){
                        $lamuda = abs($lamuda) - 180;}
                    else{
                        $lamuda = -abs($lamuda);}
                }
                $L2 = $L1 * 180 / $Pi + $lamuda - $theta * 180 / $Pi;
                $A2 = atan($cosu1 * sin($A1) / ($cosu1 * cos($q) * cos($A1) - $sinu1 * sin($q))) * 180 / $Pi;
                if (sin($A1) > 0){
                    if ($cosu1 * sin($A1) / ($cosu1 * cos($q) * cos($A1) - $sinu1 * sin($q)) > 0){
                        $A2 = 180 + abs($A2);}else{
                        $A2 = 360 - abs($A2);}
                }else{
                    if ($cosu1 * sin($A1) / ($cosu1 * cos($q) * cos($A1) - $sinu1 * sin($q)) > 0){
                        $A2 = abs($A2);}else{
                        $A2 = 180 - abs($A2);}}//A2代表北偏西的角度
    $location=array("lng"=>$L2,"lat"=>$B2);
    return $location;//返回坐标数组，格式如:("Bill"=>"35","Steve"=>"37","Peter"=>"43")
}

```