следуя этой статье, чтобы вычислить точку из широты и долготы с азимутом и расстоянием, когда результат бега в NaN?
статья ниже:
Предлагается метод расчета конечной точки для эллипсоидных моделей земли с использованием формулы Винсенти:
Преобразуйте широту начальной точки 'lat1' (в диапазоне от -90 до 90) в радианы. lat1 = lat1 * PI / 180 Преобразует долготу начальной точки lon1 (в диапазоне от -180 до 180) в радианы. lon1 = lon1 * PI / 180 Преобразование азимута brg (в диапазоне от 0 до 360) в радианы. brg = brg * PI / 180 Эллипсоидные модели земли
Примечание: переменная «плоский» ниже представляет собой полярное уплощение Земли, используемое в различных моделях эллипсоидов. Для обычно используемого WGS-84 пусть flat = 298,257223563.
Учитывая расстояние s в метрах, большую полуось «a» в метрах, малую полуось «b» в метрах и полярное сглаживание «плоское». Рассчитайте пункт назначения по формуле Винсенти. Используются сокращенные имена переменных.
f = 1/flat
sb=sin(brg)
cb=cos(brg)
tu1=(1-f)*tan(lat1)
cu1=1/sqrt((1+tu1*tu1))
su1=tu1*cu1
s2=atan2(tu1, cb)
sa = cu1*sb
csa=1-sa*sa
us=csa*(a*a - b*b)/(b*b)
A=1+us/16384*(4096+us*(-768+us*(320-175*us)))
B = us/1024*(256+us*(-128+us*(74-47*us)))
s1=s/(b*A)
s1p = 2*PI
Пока выполняется условие, выполните следующие действия.
while (abs(s1-s1p) > 1e-12)
cs1m=cos(2*s2+s1)
ss1=sin(s1)
cs1=cos(s1)
ds1=B*ss1*(cs1m+B/4*(cs1*(-1+2*cs1m*cs1m)- B/6*cs1m*(-3+4*ss1*ss1)*(-3+4*cs1m*cs1m)))
s1p=s1
s1=s/(b*A)+ds1
Продолжить расчет после цикла.
t=su1*ss1-cu1*cs1*cb
lat2=atan2(su1*cs1+cu1*ss1*cb, (1-f)*sqrt(sa*sa + t*t))
l2=atan2(ss1*sb, cu1*cs1-su1*ss1*cb)
c=f/16*csa*(4+f*(4-3*csa))
l=l2-(1-c)*f*sa* (s1+c*ss1*(cs1m+c*cs1*(-1+2*cs1m*cs1m)))
d=atan2(sa, -t)
finalBrg=d+2*PI
backBrg=d+PI
lon2 = lon1+l;
Convert lat2, lon2, finalBrg and backBrg to degrees
lat2 = lat2 * 180/PI
lon2 = lon2 * 180/PI
finalBrg = finalBrg * 180/PI
backBrg = backBrg * 180/PI
Если lon2 находится за пределами диапазона от -180 до 180, добавьте или вычтите 360, чтобы вернуть его в этот диапазон. Если finalBrg или backBrg вне диапазона от 0 до 360, добавьте или вычтите 360, чтобы вернуть их в этот диапазон. Примечание: указанные выше переменные a, b и flat имеют следующие отношения:
b = a - (a/flat)
flat = a / (a - b)
поэтому перенос этой формулы на С # я получил:
double Latitude = 50.390202;
double Longitude = -3.9204310000000078;
double Bearing = 225;
double lat1 = Latitude * (Math.PI / 180.0);
double lon1 = Longitude * (Math.PI / 180.0);
double brg = Bearing * (Math.PI / 180.0);
double s = 1000;
double a = 6378137.0;
double b = 6356752.314245;
double f = 1 / 298.257223563;
double sb = Math.Sin(brg);
double cb = Math.Cos(brg);
double tu1 = (1-f) * Math.Tan(lat1);
double cu1 = 1 / Math.Sqrt((1+tu1*tu1));
double su1 = tu1 * cu1;
double s2 = Math.Atan2(tu1, cb);
double sa = cu1 * sb;
double csa = 1 - sa * sa;
double us = csa * (a * a - b * b) / (b * b);
double A = 1 + us / 16384 * (4096 + us * (- 768 + us * (320 - 175 * us)));
double B = us / 1024 * (256 + us * (-128 + us * (74 - 47 * us)));
double s1 = s / (b * A);
double s1p = 2 * Math.PI;
while (Math.Abs(s1-s1p) > 1e-12)
{
cs1m = Math.Cos(2 * s2 + s1);
ss1 = Math.Sin(s1);
cs1 = Math.Cos(s1);
double ds1 = B * ss1 * (cs1m + B / 4 * (cs1 * (-1 + 2 * cs1m * cs1m) - B / 6 * cs1m * (-3 + 4 * ss1 * ss1) * (-3 + 4 * cs1m * cs1m)));
s1p = s1;
s1 = s / (b * A) + ds1;
}
double t = su1 * ss1 - cu1 * cs1 * cb;
double lat2 = Math.Atan2(su1 * cs1 + cu1 * ss1 * cb, (1 - f) * Math.Sqrt(sa * sa + t * t));
double l2 = Math.Atan2(ss1 * sb, cu1 * cs1 - su1 * ss1 * cb);
double c = f / 16 * csa * (4 + f * (4 - 3 * csa));
double l = l2 - (1 - c) * f * sa * (s1 + c * ss1 * (cs1m + c * cs1 * (-1 + 2 * cs1m * cs1m)));
double d = Math.Atan2(sa, -t);
double finalBrg = d + 2 * Math.PI;
double backBrg = d + Math.PI;
double lon2 = lon1 + l;
lat2 = lat2 * 180 / Math.PI;
lon2 = lon2 * 180/ Math.PI;
finalBrg = finalBrg * 180/ Math.PI;
backBrg = backBrg * 180 / Math.PI;
if (lon2 < -180)
{
lon2 = lon2 + 360;
}
if (lon2 > 180)
{
lon2 = lon2 - 360;
}
if (finalBrg < 0)
{
finalBrg = finalBrg + 360;
}
if (finalBrg > 360)
{
finalBrg = finalBrg - 360;
}
Console.WriteLine("{0} {1}", lat2, lon2);
но когда запускается результат в NaN как для lat2, так и для lon2?
ожидаемый результат должен быть:
Широта: 50 ° 23′02 ″ с.ш. 50.38384479 Долгота: 3 ° 55′49 ″ з.д. -3,93037298 Конечный азимут: 224 ° 59′32 ″ 224.99234101 Задний азимут: 44 ° 59′32 ″ 44.99234101
делает ли С # calc что-то другое, если я сделал ошибку, которую не вижу :(
Спасибо
lat2
иlon2
- вы не можете определить, что пошло не так, просто посмотрев значения переменных. Также разбейте код на отдельные функции, которые вы можете протестировать независимо. Во всем этом невозможно разобраться, не говоря уже о тестировании. - person Panagiotis Kanavos   schedule 04.04.2016