2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему
 
 Формулы сферической тригонометрии точные или приближенные ?
Сообщение02.09.2010, 07:59 
Аватара пользователя


22/05/06

358
Волгоград
Занимаясь вопросом эфемеридной поправки ET-UT для современных теорий и для теорий древних астрономов, я создал еще одну форму для своей программы Solsys7 (скриншот смотрите ниже), где у меня воспроизводится анимация солнечных и лунных затмений и транзита планет по диску Солнца. Проверяя время наступления затмений или транзита планет, я заметил, что в экваториальных координатах по анимационной картинке (нижняя картинка на скриншоте) получаются хорошие результаты, как для современных наблюдений, так и для древних, а в горизонтальных координатах наблюдается какая то странность. Большинство наблюдений тоже хорошо воспроизводилось и в горизонтальных координатах, но иногда наблюдения очень сильно расходились с расчетными данными, которые я получал, наблюдая анимационную картинку затмения или транзита.


Изображение



Скриншот формы 16 программы Solsys7, где воспроизводится транзит Венеры 8.6.2004 (на нижней картинке анимация в горизонтальных координатах), где точками 0…7 обозначены восемь точек на диске Солнца, координаты которых определены в экваториальной системе координат и затем пересчитаны в горизонтальные координаты. Кружок на диске Солнца это Венера.


В центре картинки я рисовал видимый диаметр Солнца (для солнечных затмений и транзитов), а затем по разности координат (экваториальных или горизонтальных) между центром диска Солнца и планеты или Луны находил положение на картинке центра планеты или Луны и рисовал их видимый диаметр. При этом, первоначально я получал для современных теорий Ньюкома или JPL (Лаборатория реактивного движения – подразделение НАСА) гелиоцентрические эклиптические координаты планет (или в эпохе даты или в стандартной эпохе J2000) для планет и геоцентрические эклиптические координаты Луны по теории Брауна (в эпохе даты), а затем приводил их к видимым координатам, т.е. к наблюдаемым с конкретной точки Земли. Аналогичные преобразования я делал и с расчетными данными, полученными по таблицам древних астрономов (Птолемей, Аль Хорезми, Кеплер и т.д.). На скриншоте в окошках справа отражена последовательность преобразований. Цепочка преобразований получается довольно таки длинной по этому в первую очередь я, конечно же, заподозрил какую то ошибку в коде программы, но проверка (прямое преобразование координат и обратное) показала, что ошибки нет.


Тогда я занялся проверкой стандартных функций преобразования, которые я взял из книги О. Монтенбрук и Т. Пфлегер – Астрономия с персональным компьютером и перевел их с языка программирования Pascal на язык Visual Basic 6.0, а в первую очередь я занялся функцией перевода экваториальных координат в горизонтальные и необходимой для этого функцией расчета Гринвичского звездного времени GMST (в коде функций после апострофа ` идет комментарий к коду).


Код:
Private Sub EQU_HOR(ByVal j As Integer)
TAU = ObsLMST - GeoL(j): If TAU < 0 Then TAU = TAU + 360  'часовой угол в градусах
Call EQUHOR(j) 'функция перевода координат из книги О. Монтенбрук и Т. Пфлегер
Call AlfaBetta(15, j) 'находим по декартовым координатам полярные координаты
Azimut(j) = AlfaGR(j): If Azimut(j) > 180 Then Azimut(j) = Azimut(j) - 360
Visota(j) = Betta1
End Sub



Код:
Public Sub EQUHOR(i As Integer) 'conversion of equatorial into horizontal coordinates
'аналог PROCEDURE EQUHOR (DEC,TAU,PHI: REAL; VAR H,AZ: REAL);
'(*   DEC = GeoB(I) : declination (-90 deg .. +90 deg)                 *)
'(*   TAU  : hour angle (0 deg .. 360 deg)                             *)
'(*   PHI = ObsB : geographical latitude (in deg)                      *)
Dim CS_PHI#, SN_PHI#, CS_DEC#, SN_DEC#, CS_TAU#
CS_PHI = Cos(ObsB / kUgol): SN_PHI = Sin(ObsB / kUgol)
CS_DEC = Cos(GeoB(i) / kUgol): SN_DEC = Sin(GeoB(i) / kUgol)
CS_TAU = Cos(TAU / kUgol)
X(i) = CS_DEC * SN_PHI * CS_TAU - SN_DEC * CS_PHI
Y(i) = CS_DEC * Sin(TAU / kUgol)
Z(i) = CS_DEC * CS_PHI * CS_TAU + SN_DEC * SN_PHI
'POLAR (X,Y,Z, DUMMY,H,AZ)
End Sub



Код:
Public Function LMST(MJD As Double) As Double ' рассчитывается GMST в часах
'FUNCTION LMST(MJD,LAMBDA:REAL):REAL; local mean sidereal time
'MJD = JD - 2400000.5
T2 = (MJD - Fix(MJD)) * 24: T1 = (Fix(MJD) - 51544.5) / 36525
LMST = 6.697374558 + 1.0027379093 * T2 + (8640184.812866 + (0.093104 - 0.0000062 * T1) * T1) * T1 / 3600# End Function


К сжалению, проверка показала, что функции дают те же результаты, что и примеры описанные в литературе, например, в книге Jean Meeus – Astronomical algorithmus. В книге даются расчетные топоцентрические экваториальные координаты Венеры на 19:21:00 UT для обсерватории USNO (ObsL= -77,0656 ObsB= 38,9214) - долгота 347,3193 градусов и широта -6,72 градуса при GMST=08:34:56,853. А у меня программа дает расчетные экваториальные координаты – долгота 347,3168 и широта -6,7216 при GMST=08:34:57 (секунды округлены и GMST посчитано без учета нутационной поправки). При этих исходных данных в книге получаются следующие горизонтальные координаты – азимут 68,0337 и высота 15,1249, а у меня программа дает – азимут 68,035 и высота 15,1211. Расхождения получились в пятом знаке, что вполне допустимо при небольших различиях в исходных экваториальных координатах и при отсутствие у меня нутационной поправки при расчете GMST.


Вообще-то проведенная мною проверка была просто шагом отчаяния, т.к. все эти функции я и так многократно тестировал на других формах программы Solsys и никаких замечаний они у меня не вызывали. Но так как полученные мною результаты в экваториальных и горизонтальных координатах сильно отличались (чего не должно быть в принципе) я решил воспользоваться еще и формулами перевода экваториальных координат в горизонтальные и в полярных координатах (приведенная выше процедура EQUHOR производит перевод координат в декартовой системе координат, а потом вычисляются полярные координаты). Я взял стандартные функции перевода (приводятся в любом учебнике) и написал процедуру EQU_HOR2. При этом, т.к. в языке программирования Visual Basic 6.0 нет многих тригонометрических функций (как и во многих других языках), например, нужных мне arcsin и arccos, а есть только arctan (в коде Atn), то мне пришлось написать для них код самому (используя известные тригонометрические преобразования).


Код:
Private Sub EQU_HOR2(ByVal j As Integer)
TAU = ObsLMST - GeoL(j): If TAU < 0 Then TAU = TAU + 360  ' TAU часовой угол в градусах
'H=sin(Visota(j)) – временная переменная
H = Sin(GeoB(j) / kUgol) * Sin(ObsB / kUgol) + Cos(GeoB(j) / kUgol) * Cos(ObsB / kUgol) * Cos(TAU / kUgol)
Visota(j) = Atn(H / (1 - H * H) ^ 0.5) * kUgol '               функция  arcsin -1...+1
H = Sin(GeoB(j) / kUgol) - Sin(ObsB / kUgol) * Sin(Visota(j) / kUgol)
H = H / Cos(ObsB / kUgol) / Cos(Visota(j) / kUgol) '       H=cos(Azimut(j)) – временная переменная
Azimut(j) = Atn((1 - H * H) ^ 0.5 / H) * kUgol '               функция arccos  0...+1 
AZ = Azimut(j)
If AZ >= 0 And TAU <= 180 Then Azimut(j) = 180 - AZ
If AZ >= 0 And TAU > 180 Then Azimut(j) = AZ - 180
If AZ < 0 And TAU <= 90 Then Azimut(j) = -AZ
If AZ < 0 And TAU > 90 Then Azimut(j) = AZ
End Sub



Результаты получились те же самые, что и при преобразованиях в декартовой системе координат, и я был просто в тупике не зная, что можно еще проверить. Тогда я решил просто тупо посмотреть, а как переводятся из экваториальных координат в горизонтальные не только координаты центра Солнца, но еще и восьми точек расположенных на диске Солнца. И тут оказалось, что в каких то положениях Солнца на горизонте точки расположенные по окружности Солнца в экваториальных координатах переводятся в горизонтальные правильно и мы видим ту же окружность, а в каких то положениях мы видим в горизонтальных координатах эллипс. Например, на скриншоте, показан вид Солнца при азимуте 0 градусов (направление на Юг), а при азимутах -90, +90 и +/-180 градусов приведен на рисунке ниже.


Изображение


Как видим при азимуте +/-180 градусов координаты точек переводятся так, что совпадают с диском Солнца нарисованном относительно центра Солнца, координаты которого тоже пересчитаны в горизонтальные координаты, а при азимутах -90 и +90 градусов у нас координаты восьми точек незначительно отстоят от окружности, т.е. переводятся почти без ошибки. Но в разное время года максимальная ошибка при переводе координат получается при разных азимутах и, например, в декабре она будет максимальной уже не при 0 градусов как летом, а при азимуте +/-180, а в марте и сентябре и при 0 градусов и при +/-180 градусов. Вот только при этом величина максимальной ошибки по видимому радиусу Солнца в марте и сентябре будет только 0,1 градуса, а в июне и декабре уже почти 0,4 градуса, что очень много и сопоставимо с видимым диаметром Солнца (0,5 градуса). И если мы будем обрабатывать данные наблюдений за Солнцем, которые выполнены не по центру диска, а по какому-то его краю, то мы получим очень не точные данные.



Таким образом, я понял откуда у меня получаются разные данные в разных системах координат для затмений и транзитов планет и возникает вопрос. А может быть, формулы сферической тригонометрии являются не точными, а приближенными и при определенных углах дают значительную погрешность, что я и получил. Вот ответ на этот вопрос я бы и хотел получить. Или может быть кто-то знает более точные формулы перевода экваториальных координат в горизонтальные. А если кто-то считает, что у меня в программе имеется какой то глюк, то у меня к нему будет просьба взять экваториальные координаты 8-и точек по окружности Солнца 22 июня в 12 часов местного времени для пункта наблюдения находящегося в северном полушарии, перевести их своими программами или обычным расчетом в горизонтальные и нанести их на рисунок.



Кстати, замеченный мною эффект эллипсности Солнца в горизонтальных координатах можно наблюдать и в культовой программе StarCalc (у меня версия 5.7) если отслеживать движение Солнца при очень большом масштабе (20000%…40000%), как это показано на скриншоте ниже. При этом в настройках для координат лучше выбрать для изображения неба не параллельную, а коническую сетку координат, т.к. при этом на анимационной картинке Солнце будет рисоваться в виде окружности, а его эллипсность можно будет определить по координатной сетке (при параллельной проекции будет и неравномерность масштаба координатной сетки и эллипсность Солнца, что смажет этот эффект). На скриншоте даже визуально видно, что по координатной сетке размер по азимуту примерно в два раза больше чем по высоте, а конкретно по высоте размер будет 67 градусов 17 минут – 66 градусов 45 минут = 32 минуты (реальный видимый диаметр Солнца), а по азимуту 0 градусов 40 минут – (-0 градусов 40 минут) = 80 минут. В результате получаем, что ошибка по видимому радиусу Солнца получается (80-32)/2=24 минуты, т.е. те же самые 0,4 градуса, что я получил по анимационной картинке программы Solsys7.


Изображение


Таким образом, мы еще раз убеждаемся, что при переводе экваториальных координат нескольких точек расположенных по окружности Солнца в горизонтальные мы получаем не окружность, а эллипс, чего не должно быть в принципе, т.к. при простом повороте системы координат (масштабы по осям остались неизменными) такой эффект не должен наблюдаться. Но, если у Вас при таком переводе координат получится окружность, то Вы обязательно сообщите мне какими формулами Вы пользовались. А, если какая то ошибка имеется в моих рассуждениях, то укажите где у меня ошибка в рассуждениях, по тому, что я, например, не понимаю почему программа StarCalc при использовании функции определения углового расстояния между двумя точками дает по азимуту результат 32 минуты, хотя по координатной сетке между этими же двумя точками мы получили 80 минут. Хотя, быть может, такой проблемы с формулами сферической тригонометрии и не существует и просто не хватает точности вычислений компьютера. Я, например, использую в программе 16 разрядов, а в программе StarCalc, как я догадываюсь, используется 32 разряда, но по грубому прикиду для устранения ошибки в десятых долях градуса хватило бы и 8 разрядов (или я не прав).


С наилучшими пожеланиями Сергей Юдин.

 Профиль  
                  
 
 Re: Формулы сферической тригонометрии точные или приближенные ?
Сообщение06.09.2010, 20:33 
Заблокирован
Аватара пользователя


03/05/09

288
Gomel BY
Формулы, может, и точные, но все зависит от метода вычислений.
Особенно внимательно надо смотреть, если, например, бесконечности вычитаются или матрицы обращать.
В лоб не получиться - придется теорией заняться, в чстности выделить предел, решить, а потом вокруг вертеться.
По этому поводу мне понравился Фихтенгольц.
В вашем случае, возможно, другое, но я не разбираюсь.
Есть даже наука методы вычислений на ЭВМ, если надо, могу немного литературы.

 Профиль  
                  
 
 Re: Формулы сферической тригонометрии точные или приближенные ?
Сообщение07.09.2010, 10:28 
Аватара пользователя


22/05/06

358
Волгоград
iig в сообщении #350152 писал(а):
Есть даже наука методы вычислений на ЭВМ, если надо, могу немного литературы.


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

Код:
cosa = Sin(DEC1 * qq) * Sin(DEC2 * qq) + Cos(DEC1 * qq) * Cos(DEC2 * qq) * Cos((RA1 - RA2) * qq)
Angle_Distance = ArcCos(cosa) / qq


Вот, например, получаются там координаты точки, которая имеет высоту примерно равную высоте Солнца (68,4327 и 68,4342), а азимут -0,7138 (у Солнца -0,0003) и по этим формулам получилась разность в координатах 0,2623. Но ведь при одинаковой высоте двух точек длина дуги между ними будет равна разности координат по азимуту, т.е. 0,7135. Если интересно можете посмотреть исходный код программы, который я дополнил и своим расчетом, и исполняемый файл http://ser.t-k.ru/Arhiv/TestSun.zip . Исполняемый файл это для тех, у кого нет языка Visual Basic 6.0, и они могут запускать его на исполнение, также как и все остальные программы, двойным кликом мышки по нему. Вот только, если у Вас старая версия Windows и в ней установлена виртуальная машина для 5-ой версии языка (файл msvbvm50.dll), то Вам надо будет поместить в папку WINDOWS/system32 файл msvbvm60.dll (скачать можно с моей домашней страницы http://ser.t-k.ru из раздела Программы).

С наилучшими пожеланиями Сергей Юдин.

 Профиль  
                  
 
 Re: Формулы сферической тригонометрии точные или приближенные ?
Сообщение07.09.2010, 15:55 
Заслуженный участник
Аватара пользователя


11/03/08
9490
Москва
Цитата:
Но ведь при одинаковой высоте двух точек длина дуги между ними будет равна разности координат по азимуту


Нет. Не будет. Давайте рассмотрим на Земле, взамен неба. Один путешественник на экваторе, и градус по долготе для него составит более 100 км, а второй на Полюсе, и для него изменение долготы на эту величину - сколь угодно малая величина.

 Профиль  
                  
 
 Re: Формулы сферической тригонометрии точные или приближенные ?
Сообщение08.09.2010, 14:39 
Аватара пользователя


22/05/06

358
Волгоград
Евгений Машеров в сообщении #350301 писал(а):
Цитата:
Но ведь при одинаковой высоте двух точек длина дуги между ними будет равна разности координат по азимуту


Нет. Не будет. Давайте рассмотрим на Земле, взамен неба. Один путешественник на экваторе, и градус по долготе для него составит более 100 км, а второй на Полюсе, и для него изменение долготы на эту величину - сколь угодно малая величина.


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

Код:
kMFequ = 1 / Cos(DecSun * qq)
kMFhor = 1 / Cos(HSun * qq)

Picture1.Circle ((RaSun - Ra2(i)) / kMFequ, Dec2(i) - DecSun), 0.005, RGB(0, 0, 250)
Picture2.Circle ((Az(i) - AzSun) / kMFhor, H(i) - HSun), 0.005, RGB(0, 0, 250)


С наилучшими пожеланиями Сергей Юдин.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 5 ] 

Модераторы: photon, whiterussian, Jnrty, Aer, Парджеттер, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group