Занимаясь вопросом эфемеридной поправки 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 разрядов (или я не прав).
С наилучшими пожеланиями Сергей Юдин.
