C - формула гаверсина немного отличается

Мой вывод кажется немного неверным при вычислении расстояния между двумя координатами широты и долготы, и я не могу понять, почему. Ниже мой код (реализация формулы Хаверсина)

float calcDistance(double latHome, double lonHome, double latDest, double lonDest) {

    double pi = 3.141592653589793;
    int R = 6371; //Radius of the Earth
    latHome = (pi/180)*(latHome);
    latDest = (pi/180)*(latDest);
    double differenceLon = (pi/180)*(lonDest - lonHome);
    double differenceLat = (pi/180)*(latDest - latHome);
    double a = sin(differenceLat/2) * sin(differenceLat/2) + 
    cos(latHome) * cos(latDest) * 
    sin(differenceLon/2) * sin(differenceLon/2);
    double c = 2 * atan2(sqrt(a), sqrt(1-a));
    double distance = R * c;
    printf("%f\n", distance);
    return distance;
}

Ввод: 38,898556 -77,037852 38,897147 -77,043934

Выход: 0,526339

Должно получиться 0,5492


c haversine
person Phil O'kelly    schedule 17.09.2015    source источник
comment
Было бы полезно, если бы вы привели пример ввода, фактического вывода и ожидаемого вывода.   -  person user3386109    schedule 17.09.2015
comment
Извините, добавил пример   -  person Phil O'kelly    schedule 17.09.2015
comment
Код правильный (синтаксически), убедитесь, что вы используете правильную формулу. Там нет фанкового усечения   -  person vsoftco    schedule 17.09.2015
comment
Сама формула кажется мне неправильной. Почему вы используете арктангенс? формула представляет собой arcsin(sqrt(a)) для вашей переменной a.   -  person ely    schedule 17.09.2015
comment
Что такое выходная единица? Майлз? Метры? Стержни? Парсек?   -  person stackoverflowuser2010    schedule 17.09.2015
comment
Вы сделали расчет вручную на калькуляторе и сравнили все промежуточные значения в нем с промежуточными значениями в вашей функции? Вы распечатали входные данные, чтобы убедиться, что они соответствуют вашим ожиданиям? Откуда вы знаете, что ответ должен быть именно таким, как вы говорите?   -  person Jonathan Leffler    schedule 17.09.2015
comment
Примечание: 2*atan2(sqrt(a), sqrt(1-a)); имеет только один квадрант ответа. С таким же успехом можно использовать 2*asin(sqrt(a)) для поддержания точности.   -  person chux - Reinstate Monica    schedule 17.09.2015

Ответы (1)


arrow_upward
3
arrow_downward

Единственное, о чем я могу думать (и если фиксированное дает ожидаемый результат), это тот факт, что вы перезаписываете локальные параметры latHome и latDest (но, вероятно, следует использовать исходные при расчете differenceLon и differenceLat выше a). Используйте другое имя для них, например latHomeTmp, затем вычислите a с новыми, и это сработает, например:

double latHomeTmp = (pi/180)*(latHome);
double latDestTmp = (pi/180)*(latDest);
double a = sin(differenceLat/2.) * sin(differenceLat/2.) + 
           cos(latHomeTmp) * cos(latDestTmp) * 
           sin(differenceLon/2.) * sin(differenceLon/2.);

В общем, рекомендуется сохранять параметры, передаваемые функции, неизменяемыми, чтобы избежать подобных ситуаций (конечно, это не твердое правило, но я обычно его придерживаюсь).

person vsoftco    schedule 17.09.2015
comment
Да, проблема была в том, что широты переводились из градусов в радианы, а потом разница между широтами (уже в радианах) снова конвертировалась из градусов в радианы. Таким образом, проблема также может быть решена путем удаления (pi/180) из расчета разницы широт. - person user3386109; 17.09.2015