Вопрос о большом круге

Вот реализация C, которую я нашел, должна быть достаточно простой для перевода на PHP:

 #define KmPerDegree 111.12000071117 #define DegreesPerKm (1.0/KmPerDegree) #define PI M_PI #define TwoPI (M_PI+M_PI) #define HalfPI M_PI_2 #define RadiansPerDegree (PI/180.0) #define DegreesPerRadian (180.0/PI) #define copysign(x,y) (((y)<0.0)?-fabs(x):fabs(x)) #define NGT1(x) (fabs(x)>1.0?copysign(1.0,x):(x)) #define ArcCos(x) (fabs(x)>1?quiet_nan():acos(x)) #define hav(x) ((1.0-cos(x))*0.5) /* haversine */ #define ahav(x) (ArcCos(NGT1(1.0-((x)*2.0)))) /* arc haversine */ #define sec(x) (1.0/cos(x)) /* secant */ #define csc(x) (1.0/sin(x)) /* cosecant */ /* ** GreatCirclePos() -- ** ** Compute ending position from course and great-circle distance. ** ** Given a starting latitude (decimal), the initial great-circle ** course and a distance along the course track, compute the ending ** position (decimal latitude and longitude). ** This is the inverse function to GreatCircleDist). */ void GreatCirclePos(dist, course, slt, slg, xlt, xlg) double dist; /* -> great-circle distance (km) */ double course; /* -> initial great-circle course (degrees) */ double slt; /* -> starting decimal latitude (-S) */ double slg; /* -> starting decimal longitude(-W) */ double *xlt; /* <- ending decimal latitude (-S) */ double *xlg; /* <- ending decimal longitude(-W) */ { double c, d, dLo, L1, L2, coL1, coL2, l; if (dist > KmPerDegree*180.0) { course -= 180.0; if (course < 0.0) course += 360.0; dist = KmPerDegree*360.0-dist; } if (course > 180.0) course -= 360.0; c = course*RadiansPerDegree; d = dist*DegreesPerKm*RadiansPerDegree; L1 = slt*RadiansPerDegree; slg *= RadiansPerDegree; coL1 = (90.0-slt)*RadiansPerDegree; coL2 = ahav(hav(c)/(sec(L1)*csc(d))+hav(d-coL1)); L2 = HalfPI-coL2; l = L2-L1; if ((dLo=(cos(L1)*cos(L2))) != 0.0) dLo = ahav((hav(d)-hav(l))/dLo); if (c < 0.0) dLo = -dLo; slg += dLo; if (slg < -PI) slg += TwoPI; else if (slg > PI) slg -= TwoPI; *xlt = L2*DegreesPerRadian; *xlg = slg*DegreesPerRadian; } /* GreatCirclePos() */ 

Источник: http://sam.ucsd.edu/sio210/propseawater/ppsw_c/gcdist.c

Чтобы ответить на мой собственный вопрос, так это для кого-то любопытного, класс PHP, преобразованный из функции C, предоставленной Чадом Берчем:

 class GreatCircle { /* * Find a point a certain distance and vector away from an initial point * converted from c function found at: http://sam.ucsd.edu/sio210/propseawater/ppsw_c/gcdist.c * * @param int distance in meters * @param double direction in degrees ie 0 = North, 90 = East, etc. * @param double lon starting longitude * @param double lat starting latitude * @return array ('lon' => $lon, 'lat' => $lat) */ public static function getPositionByDistance($distance, $direction, $lon, $lat) { $metersPerDegree = 111120.00071117; $degreesPerMeter = 1.0 / $metersPerDegree; $radiansPerDegree = pi() / 180.0; $degreesPerRadian = 180.0 / pi(); if ($distance > $metersPerDegree*180) { $direction -= 180.0; if ($direction < 0.0) { $direction += 360.0; } $distance = $metersPerDegree * 360.0 - $distance; } if ($direction > 180.0) { $direction -= 360.0; } $c = $direction * $radiansPerDegree; $d = $distance * $degreesPerMeter * $radiansPerDegree; $L1 = $lat * $radiansPerDegree; $lon *= $radiansPerDegree; $coL1 = (90.0 - $lat) * $radiansPerDegree; $coL2 = self::ahav(self::hav($c) / (self::sec($L1) * self::csc($d)) + self::hav($d - $coL1)); $L2 = (pi() / 2) - $coL2; $l = $L2 - $L1; $dLo = (cos($L1) * cos($L2)); if ($dLo != 0.0) { $dLo = self::ahav((self::hav($d) - self::hav($l)) / $dLo); } if ($c < 0.0) { $dLo = -$dLo; } $lon += $dLo; if ($lon < -pi()) { $lon += 2 * pi(); } elseif ($lon > pi()) { $lon -= 2 * pi(); } $xlat = $L2 * $degreesPerRadian; $xlon = $lon * $degreesPerRadian; return array('lon' => $xlon, 'lat' => $xlat); } /* * copy the sign */ private static function copysign($x, $y) { return ((($y) < 0.0) ? - abs($x) : abs($x)); } /* * not greater than 1 */ private static function ngt1($x) { return (abs($x) > 1.0 ? self::copysign(1.0 , $x) : ($x)); } /* * haversine */ private static function hav($x) { return ((1.0 - cos($x)) * 0.5); } /* * arc haversine */ private static function ahav($x) { return acos(self::ngt1(1.0 - ($x * 2.0))); } /* * secant */ private static function sec($x) { return (1.0 / cos($x)); } /* * cosecant */ private static function csc($x) { return (1.0 / sin($x)); } } 

Я бы подумал, что было бы труднее отбросить Haversine fomula, а затем создать свою собственную.

Сначала угол, созданный из ядра Земли, перемещая «прямую» линию на поверхности (вы считаете, что она прямая, но она изогнута).

Угол в радиантах = длина дуги / радиус. Угол = ArcLen / 6371 км

Широта должна быть легкой, просто «вертикальной» (северной / южной) составляющей вашего угла.

Lat1 + Cos (подшипник) * Угол

Различия по долготе варьируются в зависимости от широты. Так что это становится сложнее. Вы будете использовать:

Грех (подшипник) * Угол (с востоком определяется как отрицательный), чтобы найти угол в направлении долготы, но переход к реальной долготе на этой широте будет более сложным.

См. Раздел Точка назначения, заданное расстояние и опознавание от начальной точки на этом веб-сайте: http://www.movable-type.co.uk/scripts/latlong.html