Компас, указывающий не на север, или как мне пригодилась тригонометрия

Фундаментальный труд

Эээ… ошибок в определении направления — туча.

Начну с того, что в точке вблизи западной границы вашей картинки «компас» должен показывать на запад (Земля-то круглая!), а ваши вычисления покажут на восток. Аналогичная ситуация случается у полюсов — так, если кратчайший путь до Мекки будет проходить через северный полюс, ваша программа опять обманет с направлением.

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

Через векторы такое делается же, пусть М — точка, где находится Мекка, А — текущая точка, тогда стрелка компаса должна будет под тем же углом, что и вектор МА
МА = (Хм-Ха, Ум-Уа) (координаты мекки — координаты текущей точки)
Угол же вектора МА тоже надо от чего-то мерять, будем от оси У, тогда возьмем вектор а{0,1}
MA * a = |MA|*|a|*cos(alpha) = Xма*1 + Yма * 0 (скалярное произведение)
cos(alpha) = Xma/|MA|
alpha = arccos(Xma/(Xma^2+Yma^2))

Как-то так, формулы не проверял. Ну и то, что земля не плоскость и к любой точке можно попасть из противоположных сторон создает свои проблемы.

А теперь, как такие задачи надо решать. Во-первых, тригонометрию надо забыть. В многих ее формулах полно сингулярностей, от которых желательно избавляться, чтобы не получилось как в историях про знаменитые глюки американских военных самолетов.
Точнее, полностью избавиться от нее, конечно же, не получится — но формулы будут использоваться только простейшие.

Во-вторых, никаких предположений вида «просто возьмем и нарисуем так». Все рисунки должны иметь смысл (вы вот проверяли, как будет выглядеть гипотенуза в вашем треугольнике на поверхности планеты?)

Считать можно двумя способами: в векторах и в кватернионах. Второй способ мне в школе не рассказывали, даже на «олимпиадном» спецкурсе, поэтому я буду использовать первый как более привычный.

Возьмем систему геоцентрическую систему отсчета: начало координат — в центре Земли, оси x и y — в плоскости экватора, ось z смотрит на северный полюс. Также для первого раза примем Землю шарообразной — все одно получится точнее, чем у вас.

Сначала построим векторы a и b — это будут радиус-векторы текущего положения и Мекки. (Тут первый раз будет применена «тригонометрия», в соответствии с известными формулами перехода от сферической системы координат к декартовой).

Кратчайший путь — это дуга большого круга, большой круг — это сечение шара плоскостью, проходящей через центр. Искомое направление — это касательная к кратчайшему пути, она, очевидно, лежит в той же самой плоскости сечения. Если построить чертеж в этой самой плоскости, то будет видно, что искомое направление — это проекция вектора b на плоскость, касающуюся Земли в точке a (отмечу, что такая простая формула получается только для шара).

Получаем вектор c (нужное направление) = — (b x a) x a / R³, где R — радиус Земли (тут я вычислил проекцию вектора на плоскость через двойное векторное умножение на нормаль к плоскости, после чего нормировал результат). В принципе, в эту формулу можно подставить формулы для векторов a и b — тогда R сократится. Только вот понятности в получившейся формуле будет ноль. Так что лучше считать последовательно, положив R=1.

Теперь надо это направление перевести в угол. Для этого надо ввести вторую систему координат, местную. Ось аппликат в такой системе будет совпадать с вектором a, ось абсцисс будет касательной к параллели, а ось y — касательной к меридиану. Найдем ось абсцисс местной системы координат. Так как она — касательная к параллели, она должна лежать в ее плоскости, а нормаль этой плоскости равна k (т.е. оси аппликат геоцентрической системы). Кроме того, она должна лежать в плоскости, касающейся Земли в точке a. Прямая, которая лежит сразу в двух плоскостях, есть пересечение этих плоскостей. Ну а направляющий вектор такой прямой есть векторное произведение их нормалей.

Мы нашли ось абсцисс местной системы координат — i₁ = a x k / R. Ну, точнее, эта формула справедлива с точностью до знака — знак надо искать по «правилу буравчика» или любому другому применимому.

Ось ординат j₁ находится как произведение k₁ x i₁ = a₁ x i₁ / R.

Теперь мы знаем проекции искомого направления на местные оси — (i₁c) и (j₁c), и угол найти довольно просто. В статье он уже был найден — но можно поступить еще проще и использовать стандартную библиотечную функцию «арктангенс двух аргументов» — в разных языках программирования она называется по разному, но суть у нее одна: вычисление угла поворота заданного вектора.

И так, вот что получилось:

a = ( cos ourLattitude * cos ourLongtitude, cos ourLattitude * sin ourLongtitude, sin ourLattitude )
b = ( cos mekkaLattitude * cos mekkaLongtitude, cos mekkaLattitude * sin mekkaLongtitude, sin mekkaLattitude )
c = - (b x a) x a
i1 = (0, 0, 1) x a
j1 = a x i1
angle = arctan (i1 * c, j1 * c)

Теперь осталось модифицировать angle с учетом heading. Разумеется, все эти формулы поддаются дальнейшему упрощению — только вот проводить его нет смысла.

Разумеется, код я не проверял, да и ночь на дворе, поэтому он может содержать ошибки. Извиняюсь заранее.

Полученного результата достаточно, чтобы примерно прикинуть направление. Ну а если нужна большая точность — то пора открывать специализированную литературу и не заниматься велосипедостроением.

О, а интересно, если Мекка ровно на противоположной стороне шарика, куда в таком случае нужно лицом садиться, и как?

Если хочется серьезно работать с гироскопом в телефоне, то без базового понимания кватернионов делать там нечего. С другой стороны, не так уж это и сложно, разобраться вполне реально.

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

Например, прямое (на сфере) направление на Мекку из Нью-Йорка выглядит так: yadi.sk/i/pKdLZT9IeMj49 (красная кривая). А синяя прямая — это, что получается по Вашему алгоритму. Сами видите, какая разница.

Радуют переменные с названиями: angle1, angle2, angle3