Матан №2

RawonaM · ноября 5, 2012, 11:36

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

GaLL · ноября 5, 2012, 11:57

Цитата: RawonaM от ноября 5, 2012, 11:36
Все-таки мне кажется, что должен быть способ из данных двух точек в угловых координатах в градусах найти угловое расстояние в градусах между ними без перехода в радианы, потом в декартовы координаты.

Есть способ с координатами в радианах, благодаря сферической версии теоремы косинусов:
(wiki/en) Spherical_law_of_cosines
Если надо найти расстояние между A и B, в качестве C можно взять один из полюсов.

Hellerick · ноября 5, 2012, 12:00

Это фактически в (wiki/en) Great-circle_distance и было написано:

RawonaM · ноября 5, 2012, 12:07

Цитата: GaLL от ноября 5, 2012, 11:57
Если надо найти расстояние между A и B, в качестве C можно взять один из полюсов.

Что-то я не понял как пользоваться этой теоремой. Как в таком случае найти расстояние от С до А и В?

RawonaM · ноября 5, 2012, 12:08

Цитата: Hellerick от ноября 5, 2012, 12:00

Я пытался применить эту формулу, у меня выходят неправильные результаты (при переходе в декартову систему выходят другие), что-то я делаю не так.

GaLL · ноября 5, 2012, 12:13

Цитата: RawonaM от ноября 5, 2012, 12:07
Что-то я не понял как пользоваться этой теоремой. Как в таком случае найти расстояние от С до А и В?

Расстояние на сфере единичного радиуса от полюса до А равно π/2 минус широта точки А (или плюс, если взять противоположный полюс).

RawonaM · ноября 5, 2012, 12:15

Вот как я считываю данные:

Код Выделить

      phi   = ra * dpi/180.0;
      theta = (90.0-dec)*dpi/180;
      x[i] = sinf(theta)*cosf(phi);
      y[i] = sinf(theta)*sinf(phi);
      z[i] = sinf(cosf(theta));
      raa[i] = phi;
      deca[i] = theta;

А вот тут я печатаю для точек p и q:

Код Выделить

  theta = px*qx + py*qy + pz*qz;
  printf("theta: %f, mytheta: %f\n", theta, cosf(pra)*cosf(qra)*cosf(pdec-qdec)+sinf(pra)*sinf(qra));

Цитата: outputtheta: 0.999557, mytheta: 0.937411
theta: 0.999557, mytheta: 0.885613
theta: 0.999562, mytheta: 0.843993
theta: 0.999561, mytheta: 0.700190
theta: 0.999562, mytheta: 0.819516
theta: 0.999562, mytheta: 0.620748
theta: 0.999561, mytheta: 0.996170
theta: 0.999561, mytheta: 0.991670
theta: 0.999566, mytheta: 0.687178
theta: 0.999565, mytheta: 0.543625
theta: 0.999567, mytheta: 0.920863
theta: 0.999567, mytheta: 0.854585
theta: 0.999560, mytheta: 0.976778
theta: 0.999563, mytheta: 0.671846
theta: 0.999563, mytheta: 0.513703
theta: 0.999564, mytheta: 0.942698
theta: 0.999563, mytheta: 0.852605
theta: 0.999563, mytheta: 0.684127
theta: 0.999563, mytheta: 0.562309
theta: 0.999564, mytheta: 0.882153
theta: 0.999563, mytheta: 0.857082
theta: 0.999567, mytheta: 0.972607
theta: 0.999569, mytheta: 0.405220
theta: 0.999569, mytheta: 0.215598
theta: 0.999568, mytheta: 0.205926
theta: 0.999568, mytheta: 0.056741
theta: 0.999569, mytheta: 0.926480

Hellerick · ноября 5, 2012, 12:17

Цитата: RawonaM от ноября 5, 2012, 12:08
Цитата: Hellerick от ноября 5, 2012, 12:00

Я пытался применить эту формулу, у меня выходят неправильные результаты (при переходе в декартову систему выходят другие), что-то я делаю не так.

ЦитироватьThis arccosine formula above can have large rounding errors if the distance is small

Она действительно плохо подходит для малых расстояний. Поэтому, когда я рассчитывал картографические проекции, мне приходилось просчитывать именно через декартовы координаты.

RawonaM · ноября 5, 2012, 12:19

Понятно, значит похоже надо делать через декартову.

GaLL · ноября 5, 2012, 13:07

Цитата: RawonaM от ноября 5, 2012, 12:15
printf("theta: %f, mytheta: %f\n", theta, cosf(pra)*cosf(qra)*cosf(pdec-qdec)+sinf(pra)*sinf(qra));

Для вычислений лучше не юзать float'ы, так как плохо сохраняют точность. Впрочем, зависит от характера задачи.

RawonaM · ноября 5, 2012, 13:09

Цитата: GaLL от ноября 5, 2012, 13:07
Для вычислений лучше не юзать float'ы, так как плохо сохраняют точность. Впрочем, зависит от характера задачи.

Да я в курсе, но почему-то в этом коде уже было так и я не буду это менять.
Мне только надо распараллелить этот код для кластера, просто я подумал было оптимизировать вычисления, не мог понять зачем они в декартову переходят, однако похоже этим вопросом уже до меня задавались.

GaLL · ноября 5, 2012, 13:11

Цитата: RawonaM от ноября 5, 2012, 12:19
Понятно, значит похоже надо делать через декартову.

Тут явно не в этом дело.

Цитировать
theta = px*qx + py*qy + pz*qz;

А почему тета равна просто скалярному произведению векторов?

RawonaM · ноября 5, 2012, 13:17

Цитата: GaLL от ноября 5, 2012, 13:11
Цитата: RawonaM от ноября 5, 2012, 12:19Понятно, значит похоже надо делать через декартову.
Тут явно не в этом дело.

Да я тоже так думаю, но сейчас не могу выделить на это время.

Цитата: GaLL от ноября 5, 2012, 13:11
Цитироватьtheta = px*qx + py*qy + pz*qz;
А почему тета равна просто скалярному произведению векторов?

На самом деле theta это не тета, а cos(тета). т.е. настоящая тета там дальше берется арккосинусом.

GaLL · ноября 5, 2012, 13:31

Цитата: RawonaM от ноября 5, 2012, 12:15
printf("theta: %f, mytheta: %f\n", theta, cosf(pra)*cosf(qra)*cosf(pdec-qdec)+sinf(pra)*sinf(qra));

pra и qra отсчитываются от полюса или от экватора?

RawonaM · ноября 5, 2012, 13:35

Цитата: GaLL от ноября 5, 2012, 13:31
Цитата: RawonaM от ноября 5, 2012, 12:15printf("theta: %f, mytheta: %f\n", theta, cosf(pra)*cosf(qra)*cosf(pdec-qdec)+sinf(pra)*sinf(qra));
pra и qra отсчитываются от полюса или от экватора?

От полюса, они тут вычисляются:

Код Выделить

      phi   = ra * dpi/180.0;
      theta = (90.0-dec)*dpi/180;
      x[i] = sinf(theta)*cosf(phi);
      y[i] = sinf(theta)*sinf(phi);
      z[i] = sinf(cosf(theta));
      raa[i] = phi;
      deca[i] = theta;

Тут запутано с названиями переменных, по-хорошему надо было по-другому называть, прошу прощения

GaLL · ноября 5, 2012, 13:38

Раз от полюса, то должно быть:
sinf(pra)*sinf(qra)*cosf(pdec-qdec)+cosf(pra)*cosf(qra)

RawonaM · ноября 5, 2012, 13:47

Цитата: GaLL от ноября 5, 2012, 13:38
Раз от полюса, то должно быть:
sinf(pra)*sinf(qra)*cosf(pdec-qdec)+cosf(pra)*cosf(qra)

Поменял, вот результат:

Цитироватьtheta: 0.999560, mytheta: 0.989693
theta: 0.999563, mytheta: 0.633222
theta: 0.999563, mytheta: 0.428785
theta: 0.999564, mytheta: 0.945083
theta: 0.999563, mytheta: 0.884490
theta: 0.999563, mytheta: 0.586955
theta: 0.999563, mytheta: 0.530499
theta: 0.999564, mytheta: 0.919175
theta: 0.999563, mytheta: 0.857468
theta: 0.999567, mytheta: 0.439090
theta: 0.999569, mytheta: 0.408172
theta: 0.999569, mytheta: 0.270013
theta: 0.999568, mytheta: 0.163122
theta: 0.999568, mytheta: 0.054385
theta: 0.999569, mytheta: 0.987898

GaLL · ноября 5, 2012, 14:08

Цитата: RawonaM от ноября 5, 2012, 12:15
phi = ra * dpi/180.0;
theta = (90.0-dec)*dpi/180;
x = sinf(theta)*cosf(phi);
y = sinf(theta)*sinf(phi);
z = sinf(cosf(theta));

Надо z = cosf(theta);

GaLL · ноября 5, 2012, 14:34

Цитата: Квас от июля 22, 2011, 13:45
Цитата: RawonaM от июля 22, 2011, 13:35
Цитата: Квас от июля 22, 2011, 13:29В обосновании надо быть аккуратным, потому что ряды с общим членом, эквивалентным 1/n, расходятся. Как указано в вики, железобетонно будет через предел.
Так я через предел и нашел. А как еще?

Да мало ли? Например,
$[tex] \sum_{k=1}^{\infty}\frac{1}{(k+1)(k+2)}=\sum_{k=1}^{\infty}\left(\frac{1}{k+1}-\frac{1}{k+2}\right)=\\=\left(\frac{1}{2}-\frac{1}{3}\right)+\left(\frac{1}{3}-\frac{1}{4}\right)+\left(\frac{1}{4}-\frac{1}{5}\right)+\ldots=\\ =\frac{1}{2}+\left(-\frac{1}{3}+\frac{1}{3}\right)+\left(-\frac{1}{4}+\frac{1}{4}\right)+\ldots=\frac{1}{2} [/tex]$

Кстати, разложение на элементарные дроби работает и вообще для $[tex]\sum_{k=1}^{\infty}\frac{1}{P(x)}, [/tex]$ где P(x) - многочлен k-ой степени с k различными действительными корнями. Это нетрудно показать при помощи алгоритма разложения на элементарные дроби.

Цитировать
Это решение не годится без серьёзных пояснений, потому что производится перегруппировка членов условно сходящегося ряда.

Достаточно записать в таком виде i-ую частичную сумму, там получится
$[tex] \frac{1}{2}-\frac{1}{i+2} [/tex]$
стремящееся к 1/2.

RawonaM · ноября 5, 2012, 15:59

Цитата: GaLL от ноября 5, 2012, 14:08
Надо z = cosf(theta);

Действительно, спасибо. Я не перепроверял их алгоритм, а стоило.

Но результат теперь выходит странный:

Цитироватьtheta: 1.000000, mytheta: 0.989693
theta: 0.999999, mytheta: 0.633222
theta: 0.999999, mytheta: 0.428785
theta: 0.999999, mytheta: 0.945083
theta: 0.999998, mytheta: 0.884490
theta: 0.999999, mytheta: 0.586955
theta: 0.999999, mytheta: 0.530499
theta: 0.999999, mytheta: 0.919175
theta: 0.999998, mytheta: 0.857468
theta: 1.000000, mytheta: 0.439090
theta: 1.000000, mytheta: 0.408172
theta: 1.000000, mytheta: 0.270013
theta: 1.000000, mytheta: 0.163122
theta: 1.000000, mytheta: 0.054385
theta: 1.000000, mytheta: 0.987898

GaLL · ноября 5, 2012, 16:27

Это выводятся только те случаи, когда два варианта вычислений сильно расходятся?
Можно попробовать заменить float'ы на double'ы. Или добавить
if (fabs(pra - qra) < eps && fabs(pdec - qdec) < eps)
{ вычислить расстояние через декартовы координаты }
что не отразится значительно на времени работы, если точки лишь изредка расположены близко.

GaLL · ноября 5, 2012, 17:18

А вообще, вычисление синуса и косинуса раньше было гораздо медленнее, чем сложение и умножение вещественных типов данных. Так что здесь ещё надо разбираться, насколько замедляет вычисление расстояния переход к декартовым координатам.
Если требуется вычисление расстояния между всеми парами точек, то выгоднее сначала перевести координаты в декартовы, а потом уже считать все пары расстояний. Если каждая точка встречается в запросах на расстояние в среднем раз 10 и более, то, очевидно, тоже (если конечно они уже все известны заранее).

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

RawonaM · ноября 5, 2012, 18:11

Цитата: GaLL от ноября 5, 2012, 16:27
Это выводятся только те случаи, когда два варианта вычислений сильно расходятся?

Да они везде сильно расходятся, это только отрывок вывода.

Цитата: GaLL от ноября 5, 2012, 17:18
Если требуется вычисление расстояния между всеми парами точек, то выгоднее сначала перевести координаты в декартовы, а потом уже считать все пары расстояний. Если каждая точка встречается в запросах на расстояние в среднем раз 10 и более, то, очевидно, тоже (если конечно они уже все известны заранее).

Ну видимо поэтому они так и сделали, наверное.

Лингвофорум

Матан №2

RawonaM

GaLL

Hellerick

RawonaM

RawonaM

GaLL

RawonaM

Hellerick

RawonaM

GaLL

RawonaM

GaLL

RawonaM

GaLL

RawonaM

GaLL

RawonaM

GaLL

GaLL

RawonaM

GaLL

GaLL

RawonaM

Быстрый ответ