English ver. (original)

Программный инструментарий НАСА для расчета траекторий небесных тел Эфемериды DEххх JPL + интерфейс NAIF SPICE ... (The draft is in work! Черновик в работе!) <<< в начало сайта >>>

+38(066) 77-414-78

    Аннотация: Примеры использования инструментария SPICE подразделения NASA NAIF, или MICE в среде MatLab для расчетов траекторий искуственных и естественных космических тел солнечной системы в прошлое и будущее. Проводится сравнительный анализ точности доступного инструментария. В упрощенном виде представлены инструкции по установке и запуску MICE. Приводится множество примеров по практическому использованию.

Полтавский Саша (vasnas). Ставрополь-Курсавка-Крым. 2012-13 гг., e-mail: vasnas@ya.ru; cell: +7 (918) 868-22-68; +38 (066) 77-414-78

Введение.

    По 3-м условиям (высота луны и солнца + угол между ними), описанным в (CambridgePress.СС3) нужно было вычислить интервалы видимости новолуний , но астропроцессор ZET не давал доступа  ко всем 3-м условиям сразу и, даже с помощью его разработчика, я не смог приспособить программу к расчету. Решил посчитать напрямую из 'швейцарских эфемерид', от которых и отталкивается ZET, но фирменный интерфейс SWISS оказался платным (тестовая страничка).
     Эфемериды это 'опорные точки' траекторий небесных тел в прошлом и будущем, которые с высокой точностью заранее рассчитал суперкомпьютер, а координаты на конкретные моменты между опорными точками досчитывает Ваш персональный компьютер. Популярные 'швейцарские эфемериды' - основаны на эфемеридах НАСА.
    Обратился к сайту НАСА, где нашел on-line интерфейс HORIZONS System для расчета одного небесного тела, а на домене другого подразделения НАСА фирменный бесплатный интерфейс.  Выпускает эфемериды к небесным телам в NASA лаборатория JPL (Jet Propulsion Laboratory).  'Калифорнийский институт технологий' ( caltech.edu ), по контракту с НАСА, выпускает интерфейс SPICE (MatLab) - набор программ для доступа к эфемеридам (то же есть для сред Си и Фортран). Проект ведет подразделение NAIF (NASA's Navigation and Ancillary Information Facility - навигационные и вспомогательные возможности НАСА)

     Имея скудные астрономические познания, я потратил 5 дней на то, чтобы разобраться в англоязычной документации НАСА, написать программу и сделать свой расчет новолуний. Мне будет приятно, если я сэкономлю кому-то эти 5 дней. Получить помощь непосредственно у разработчиков и опытных пользователей можно послав свой вопрос на английском языке на адрес mailman-owner@naif.jpl.nasa.gov.
     В системе SPICE, для нахождения интервалов времени, удовлетворяющим нескольким условиям (расстояние, угол, перекрытие), массив временных интервалов последовательно сужается по условиям, от одной функции поиска событий к другой.

    В далекое прошлое на 32-разрядных процессорах система пока не считает. 4 апреля я обратился в NAIF и этот недостаток к началу мая 2012 в новом выпуске №65 обещали устранить. Потом выпуск был отложен к концу мая, а затем к концу лета, а теперь на 2013 год.

Сравнение точности инструментария

   С 2008 года я привык к астропроцессору ZET, это очень удобная программа - планетарий. Но программа лишь досчитывает траекторию между опорными точками эфемерид, расчитанных учеными. ZET вычисляет на основе 'швейцарских эфемерид'. Практически сравниваются не программы, а теории и эфемериды. Программа лишь дополнительный источник ошибок.

Переработанные насовские эфемериды продает швейцарская фирма, коммерческое название 'швейцарские эфемериды', видимо, поддерживает бренд дорогих и высокоточных швейцарских часов.

Самая подвижная звезда неба - 61 Cygin (61 Лебедя)

     В ZET '61 Cygin' за 1000 лет переместилась на 'пол неба', а в реальности лишь на доли градуса :

Снимок 61-я Лебедя (61 Cygin) Телескоп Хаббл 2002-10-16 20:35:37

PPM+SPICE+Matlab

Расчет движения звезд по небу

Функция get_star к программе.

Архив файла ppm_stars_vm3.def - звезды до 3 величины из каталога PPM

     Счастливые обладатели программы ZET лет 12 не замечали ошибок в расчете звезд, потому что ошибки расчета малоподвижных звезд не бросаются в глаза, а принципиально программу никто проверить не удосужился. Оказалось, что о какой-либо точности говорить не приходится, звезда 61 Cygin пролетела по небесной сфере 'пол неба', ~90°, всего за ~1000 лет. Разработчик признал ошибку : "22.04.2012 21:36 ... Эта ошибка возникает при использовании для расчёта звёзд Швейцарских эфемерид (...). При использовании формул Меёса такой ошибки не возникает.", но 'свалил' всё на 'швейцарские эфемериды'. У НАСА я не нашел эфемерид для звезд, поскольку движение звезд примитивно и его параметры просто хранятся в звездных каталогах (см. Подключение звезд), движение расчитывается по линейным формулам. Параметры звезды 61 Cygin взяты мною из каталога PPM на эпоху b1950 и рассчитаны по приведенным в starsb.ug формулам на время снимка телескопа Хаббл (см. фото).

Касания и покрытия небесных тел

   Касания и покрытия небесных тел, лучший способ сверки времени, они длятся минуты и секунды.
Данные на сайте НАСА http://eclipse.gsfc.nasa.gov/transit/transit.html .
Международная ассоциация времени покрытий http://www.lunar-occultations.com/iota/iotandx.htm

Сравнение SPICE, ZET и наблюдений.

Наблюдение касания Луны и Сатурна 03.11.2001 г. 22:05:24 UT Время ZET касания Луны и Сатурна 03.11.2001 г. 21:08 GMT Первый контакт (к источнику) Через 43 сек.	ZET: как показатель невероятной точности своей программы, разработчик приводит пример покрытия Луной Сатурна (сохраненный скриншот страницы). В реальности Сатурн появился из-за темного лимба Луны только в 03.11.2001 г. в 22:05:24 UT, через час от времени ZET (наблюдение с Lat 54.3045278 (54°18'16.30"N) Lon -0.5112778 (0°30'40.6"W))      Разработчик, избегая показать скриншот программы, на котором Сатурн за центром Луны, подстроив время программы на ~1 час раньше, сравнивает расстояние Сатурна от центра Луны по высотам и склонениям выданным ZET и системой HORIZONS NASA, пишет: "Отсюда можно определить, что разница в горизонтальных координатах между центрами дисков Луны и Сатурна составляет 18.25', радиус диска Луны 15.72', т.е. центр Сатурна находится на 2.5' дальше от края диска Луны. Радиус диска Сатурна 10". Это соответствует фотографии." (фото слева) - т.е. разработчик ZET уравнял угловые минуты ' с секундами '' ( ' = " , но в минуте 60 секунд, 1' = 60") . На самом деле Сатурн в этот момент в ~7.5 своих диаметров от края Луны - в центре Луны, т.к. 2.5' содержат 150", а диаметр самого Сатурна только 10". Несколько лет никто не замечал, что 1=60 и мне едва удалось убедить самого разработчика в этой умышленной или неумышленной многоходовой подтасовке. Касание в ZET происходит, но за ~час до реального, а это погрешность неприемлемая. SPICE, результаты : Касание MOON с SATURN, наблюдатель 54°18'16.3001" N 0°30'40.6001" W 03.11.2001 21:33.83 - 03.11 22:28.17 (UTC), 21:34.90 - 22:29.24 (TDT) Покрытие MOON с SATURN, наблюдатель 54°18'16.3001" N 0°30'40.6001" W 03.11.2001 21:35.04 - 03.11 22:26.96 (UTC), 21:36.11 - 22:28.03 (TDT) Выход из-за Луны (изоб. повернуто) Барицентр Сатурна (SPICE DE421): Положение планеты : 6 Место наблюдения: 54.30N 0.51W (54°18'16.3001"N °30'40.6001"W)В момент : 03.11.2001 22:05:24 A.D. юлианский день № JD 2452217.42 ЭКЛИПТИКАЛЬНАЯ сетка координат : Долгота, ° : 0.00000049 (0°0'0.0017735") Широта, ° : -1.88182970 (-1°52'54.5869") Широта (геодез.), ° : -1.88182977 (-1°52'54.5872") ЭКВАТОРИАЛЬНАЯ сетка координат : Прямое восх.,HR:MN:SC: 23:58:55.82 Прямое восх., : 0.00000049 (0°0'0.0017462") Склонение, ° : 20.57593351 (20°34'33.3606") ГОРИЗОНТАЛЬНАЯ сетка координат : Азимут, ° : 283.47272126 (283°28'21.7965") Высота, ° : 34.77169074 (34°46'18.0867") Высота (геодез.), ° : 34.77169167 (34°46'18.09")
	ZET: Критерий поиска: 0N 42E, Луна.0.Венера % непонятно по каким критериям поиск 13.12.2001 5:24(GMT+3)-15.12.2001 0:0(GMT+3) % взята ~половина углового размера луны ZET: Критерий поиска: Луна.[0-0.15].Венера 14.12.2001 4:31 (GMT) - 14.12.2001 5:02 (GMT) SPICE (эфемериды DE421): Наблюдатель 54°18'16.3001" N 0°30'40.6001" W Касание  : 14.12.2001 04:27.26 - 05:16.74 (UTC) Покрытие : 14.12.2001 04:27.83  - 05:16.13  (UTC)
2004 г. Венера прошла по диску Солнца Подробнее (Википедия)	Наблюдения прохождений Венеры по диску Солнца Дата место наблюдения наблюдения SPICE (UTC) примечание касание вошла. появл. отрыв касание вошла. появл. отрыв 04.12.1639 53°41′12″.9N, 2°48′54″.3W   15:15     14:57.30 15:15.68 21:36.70 21:55.39 Джереми Хоррокс 06.071761 51°30′09″.9N, 3°37′01″E                 Jan de Munck in Middelburg 1769                     09.12.1874                     06.12.1882                     08.06.2004 Лондон 05:19.9 05:39.6 11:04.1 11:23.5 05:19.99 05:39.72 11:04.17 11:23.54   06.06.2012                     ??.12.2117                     ??.12.2125                     Сравнение точности (прохождение Венеры по диску Солнца 8 июня 2008 г., UTC) Место наблюдения SPICE ZET касание вошла. появл. отрыв касание вошла. появл. отрыв касание вошла. появл. отрыв Геоцентр 05:13.5 05:32.8 11:06.6 11:25.9 05:13.57 05:32.85 11:06.64 11:25.92 - - - - Лондон 05:19.9 05:39.6 11:04.1 11:23.5 05:19.99 05:39.72 11:04.17 11:23.54 рис. рис. рис. рис. Москва 05:18.8 05:38.3 11:02.1 11:21.4 - - - - - - - - Каир 05:19.7 05:39.0 11:04.6 11:23.6 05:19.82 05:39.09 11:04.68 11:23.66 - - - - Время наблюдений взято: John Westfall A.L.P.O. Mercury/Venus Transit Coordinator В астропроцессоре ZET нет доступа к точным расчетам, но в программе, подстроив время, можно графически посмотреть на результаты полученные в SPICE.
Покрытие Венерой Меркурия 29.5.1737 г Гринвич, как это показывает астропроцессор ZET	Покрытие Меркурия Венерой 28.05.1737 г., UTC Место наблюдения SPICE ZET касание н.покр.. к.покр отрыв касание н.покрыт.. к.покрыт отрыв касание н.покр.. к.покр отрыв Гринвич - 21:52.06 - - 21:44.10 21:46.75 21:51.98 21:54.63 рис. рис. рис. рис. Обсерватория Гринвич, 51°29' (51.48333...)N 0°0'0" E: В 21:44:00 Джон Бевис отметил, что Меркурий "не ближе чем 1/10 части диаметра Венеры", набежали облака, в 21:52:06 "Венера сияет ярче, Меркурий фактически полностью скрылся за Венерой", опять облака. Примечание: Часы Бевиса спешили на 3 мин 57 сек. (Sky and Telescope, 09.1986, p220-222)
	Покрытие Луной теты Весов (PPM230911). Галилей. место: Падуя 45°25'0"N 11°53'0"E SPICE (DE406) ZET (SWISS) центры ~ касание касание Время,19.1.1610 UT 6:05:40 6:35 5:44:40 ☽ RA 238°26'41" 238°36'33" 233°10'54" Θ Lib RA 238°26'41" 238°26'41" 232°55'43" ☽ DEC -16°48'56" -16°52'21" -15°32'52" Θ Lib DEC -16°44'35" -16°44'35" -15°31'1"

     Таблица известных соединений

дата	время	место	описание	Примечание о наблюдателе и наблюдении
05.02.1570	7:47		Венера Юпитер
03.10.1590	4:55		Венера покрыла Марс,	Мёстлин
19.01.1610			покрытие Луной теты Весов (PPM230911)	Галилей
04.01.1613	2:08		Юпитер Нептун
15.08.1623	17:03		Юпитер Уран
19.09.1702			Юпитер покрыл Нептун
14.07.1708			Меркурий покрыл Уран
04.10.1708			Меркурий прошёл перед Юпитером
28.05.1737			Венера покрыла Меркурий
29.08.1771			Венера прошла перед Сатурном
21.07.1793			Меркурий покрыл Уран
09.12.1808			Меркурий прошёл перед Сатурном
03.01.1818			Венера прошла перед Юпитером

     Марс '01.01.2000 г. 12:00:00 TDT', место: с. Курсавка, '44.46N 42.51E'.

Параметр	SPICE						ZET 9
	Эфемериды DE406			Эфемериды DE421			швейц. эф.
	UTC	TDT	TDB	UTC	TDT	TDB	UTC
Долгота экл., °	327°58'1.075"	327°57'59.0156"	327°57'59.0156"	327°58'1.0743"	327°57'59.0149"	327°57'59.0149"	327°57'47.88"
Широта экл., °	-1°4'7.9653"	-1°4'8.0034"	-1°4'8.0034"	-1°4'7.9655"	-1°4'8.0036"	-1°4'8.0036"	-1°4'4.02"
Широта экл (геодез.), °	-1°4'7.9659"	-1°4'8.0039"	-1°4'8.0039"	-1°4'7.9661"	-1°4'8.0042"	-1°4'8.0042"	-
Азимут, °	351°16'34.5553"	350°58'16.4865"	350°58'16.4877"	351°16'34.5561"	350°58'16.4872"	350°58'16.4885"	351°20'13.37"
Высота, °	31°56'39.4527"	31°54'52.6059"	31°54'52.606"	31°56'39.4523"	31°54'52.6055"	31°54'52.6056"	31°56'38.21"
Высота (геодез.), °	31°56'39.467"	31°54'52.6202"	31°54'52.6203"	31°56'39.4666"	31°54'52.6198"	31°54'52.6199"	-
Удаление, А.Е.	1.8496652165239409	1.8496612070034069	1.8496612070079503	1.8496652160563021	1.8496612065357709	1.8496612065403140	1.849669

Вывод: Разница до 4 угловых секунд при расчете на ~12 лет назад.

Разница несколько дней по Марсу,при расчете минус ~32 года назад.

     Moon , location: Astrakhan 46°21'N (46.35) 48°03'E (48.05), time in TDT, J2000.0 r.frame

Параметр	JPL SPICE			SWISS ephemeris
	эфемериды DE406		эф. DE422	astro-processor ZET (SWISS)		SWISS test p.
	1.1.2000 12:00 JD 2451545	1.1.1000 12:00 JD 2086307	1.1.2000 12:00 JD 2451545	1.1.2000 12:00 JD 245154-5	2.1.1000 12:00 JD 2086309	1.1.1000 12:00 JD 2086308
Долгота экл., ° (LON)	222°55'50"	210°57'58"		222°55'29.65"	209°33'31.69"	210°51'22.2"
Широта экл., ° (LAT)	4°24'53"	0°2'1"		4°24'58.01"	-0°16'8.23"	-0° 8'15.99"
Широта (geodetic.), °	4°24'55"			-	-
Удаление, А.Е.	0.0027013587129823		0.0027013587215472	0.0026899914975339		0.002644007

Сравнение по Солнцу в 0-1 градусах Овна с 1.1.2000 по 1.1.2010 г. (формула для ZET Солнце.ARI[0-1], с шагом 10 мин)

Сравнение точности эфемерид НАСА между собой

     Вот как НАСА оценивает точность своих эфемерид DE421 и уточненных по данным межпланетных миссий DE423 (источник)

DE423 vs. DE421

Разница в километрах в прямом восхождении (RA) Венеры относительно Земли между эфемеридами DE421 и DE423	Разница в километрах в склонении (DEC) Венеры относительно Земли между эфемеридами DE421 и DE423
Разница в километрах в расстоянии (DIST) Венеры относительно Земли между эфемеридами DE421 и DE423	Разница в метрах в расстоянии (DIST) Венеры относительно Земли между эфемеридами DE421 и DE423
Ковариантность в километрах в прямом восхождении (RA) Меркурия относительно Земли, между эфемеридами DE421 и DE423	Ковариантность в километрах в склонении (DEC) Меркурия относительно Земли, между эфемеридами DE421 и DE423

Установка SPICE

     Здесь рассматривается версия SPICE для MatLab. После подключения библиотек SPICE написанных на С и подключения файлов эфемерид и др. , SPICE используется в Matlab вызовом функций. MatLab можно найти в "торрентах" (несколько Гб)

Что скачать?

     Интерфейс SPICE для MatLab: http://naif.jpl.nasa.gov/naif/toolkit_MATLAB.html Набор распаковывается, пути поиска к нему добавляются в переменные окружения MatLab, в командной строке >> addpath('c:\путь\mice\src\mice\'), >> addpath('c:\путь\mice\lib\' ), setenv MATLABPATH путь/mice/src/mice/ или в меню File->Set path. Проверим: ' >> which mice' должен вывести путь к файлу mice.mex*, если что-то не так ' 'mice' not found' или >> cspice_tkvrsn('toolkit'), выдаст версию продукта.
     В файлах SpiceGF.h, gfdist.c, gfposc.c, gfoclt.c, gfrfov.c, gfrr.c, gfsep.c, gfsntc.c, gfsubc, gftfov.c, zzgflong.c и др. исходного пакета константа погрешности расчетов SPICE_GF_CNVTOL установлена 1.e-6, что ограничивает расчеты в диапазоне десятилетий. В пакете для Windows x86 + Matlab я изменил значение на 10.e-4 и, с помощью инженера НАСА Эдварда Райта, перекомпилировал пакет. Измененными файлами необходимо перезаписать оригинальные (их копии сохраните).

     * - этот файл был испорчен на сайте НАСА, по моей просьбе 12.05.2012 г. его заменили.
     Светлые строки - файлы для основных расчетов.
     Если ошибка ссылки на закачку, значит файлы перемещены, ищите отсюда : ftp://naif.jpl.nasa.gov/pub/naif/generic_kernels/      (старые версии эфемерид)

     Барицентры планет, например, 3 - барицентр Земли, это точки и не подходят для операций с поверхностями планет (соединения, покрытия, светимость) и их ориентацией в пространстве. Поэтому для планет нужно скачать эфемериды с этими параметрами, например, 399 - геоид Земли. Программа brief.exe, входящая в набор для MatLab показывает в командной строке параметры бинарных (не текстовых) файлов данных (эфемерид и др.) системы SPICE. Например, '> brief de421.bsp' , выдает:

Загружается в программе cspice_furnsh( 'путь\standart.tm' );      

Типы текстовых файлов данных системы SPICE:
 • Spacecraft Clock kernels (SCLK) - параметры часов космического корабля
 • Leapseconds kernel (LSK) - файл временных поправок
 • Physical Constants kernels (PCK) - Ядро физических констант (обычно текстовый файл)
 • Instrument parameter kernels (IK) - Ппараметры инструмента (космического корабля, наземной станции и т.п.)
 • Frame definition kernels (FK) - Ядро определения систем координат
 • Some E-kernels (EK, although EKs are now rarely used)
 • Meta-kernels (MK, also called "furnsh kernels")
Типы бинарных файлов:
 • SP-kernels (SPK)
 • Binary-style (just a few) PC-kernels (PCK)
 • C-kernels (CK)
 • Some spacecraft Events kernels (EK/ESP) - Параметры некоторых событий коспических кораблей

Системы координат.

Системы координат

С. к. : J2000, ECLIPJ2000, CALACTIC, FK4, TOPO

Изменение координат Солнца в FK, 01.01.2000 - 01.01.2001 На графике видно, что если бы Земля не вращалась вокруг своей оси, Солнце бы делало вокруг наблюдателя на поверхности Земли, например находящегося на северном полюсе, 1 оборот в год, т.к. Земля наклонена на ~23.5° к плоскости своей орбиты, то угол между Солнцем и осью Земли колеблется в этих пределах.

Изменение долготы и широты Солнца в с.к. FK4 в моменты весенних равноденствий 1970-2012 гг., проявляются нутация и прецессия.

    'FK4' (фундаментальный каталог 4) - система координат неподвижна относительно звезд, но прецессия не учитывается, привязана к центру или топоточке Земли на момент B1950.0 ~совмещена с ECLIPJ2000, но потом не вращается с геоидом, как ECLIPJ2000, J2000 и TOPO.

     'ECLIPJ2000' - эклиптикальная С.К., широта из плоскости орбиты Земли (эклиптики), долгота 0 в точке весеннего равноденствия, т.к. С.К. вращается относительно звезд (С.К. FK4+прецессия) по долготе с периодом ~ 25 765 лет в плоскости орбиты Земли. Начало в центре Земли.

     'J2000' - экваториальная с. к. - , RA и DEC от экватора Земли. Начало в центре Земли. Можно представить эту систему координат горизонатальной системой координат 'TOPO', которая смещена в сфере на широту и долготу наблюдателя.

     'TOPO'- горизонтальная с.к., высота над горизонтом и азимут. Начало в точке наблюдения.

Относительно других объектов, например, Марса от экватора и IAU вектора от J2000 с.к -'MARSIAU', и т.д.

     Фиксированные к телам IAU модели вращения поставляются в PCK файлах, для Земли низкой точности 'IAU_EARTH' (всегда доступна) и высокой 'ITRF93', только на сравнительно близкие промежутки времени (нужны PCK файлы и подключение), для Марса 'IAU_MARS' и т.д.

    'IAU_EARTH' и 'ITRF93' переключаются через алиас 'EARTH_FIXED' в файле my_frames.tk :

Галактическая Система II координат 'GALACTIC' в плоскости галактики (центр Солнце, начало координат галактический центр). На B1950 RA = 12ч 49м (192°,25) DEC = +27,4°.

В файле my_sites.tpc, мною введены точки наблюдения : 'DYN', 'KURSAVKA' , 'ASTRAKHAN' , 'CAIRO', 'ISTANBUL', 'LONDON'.
'DYN' - зарезервировано мною для точки и С.К. с динамически именяющимися координатами.

Названия систем координат получаются добавлением '_TOPO', они описаны в my_frames.tk

Классы систем координат :

1. Внутренние, инерционные с.к.. Не вращаются относительно звездного окружения. К ним можно приложить законы ньютоновского движения.

2. Прикрепленные к телу) с.к. Необходимо сначала загрузить PCK файл в котором находятся параметры самого тела.

3. CK - относятся те с.к. ориентация которых соотносится к другим с.к. поддерживаемыми SPICE C-kernel. ... C-kernels использует таймер космического аппарата 'ticks'. ...

4. Фиксированные с.к. к другой с.к. с определенным смещением. Текстовые файлы - Text Kernel (TK) с.к.

5. Динамические с.к. Смещение относительно другой с.к. описывается формулами.

      Главный источник земной нутации приливные силы Луны, которые постепенно изменяют направление и таким образом
расшатывают земную ось. Этот компонент нутации имеет 18.6 летний цикл (6798 дней), за который узлы лунной орбиты делают оборот вокруг Земли. Нутация достигает ∓ 17″ по долготе и 9″ уклона. 183 дневный период (0.5 года), имеет амплитуду ∓ 1.3″ и 0.6″ соответственно. (Wikipedia) Пример расчета нутации - открыть (Ерана О. Офек).

Прецессия и динамический фрейм её учитывающий.

Динамический фрейм всегда привязан к базовому и определяется вращением вокруг осей бозового фрейма с течением времени.

Терминология:

Прямое восхождение — астрономический эквивалент земной долготы во второй экваториальной системе. И прямое восхождение, и долгота измеряют угол «восток-запад» вдоль экватора; обе меры берут отсчёт от нулевого пункта на экваторе. Начало отсчёта долготы на Земле — нулевой меридиан; начало отсчёта прямого восхождения на небе — точка весеннего равноденствия.
Склонение (δ) угол от плоскости небесного экватора, + к северу, - к югу.

Как добавить точку наблюдения?

    Необходимо приготовить текстовый файл с параметрами, где :
SITES = ('<имя1>','<имя2>', ... '<имяN>') - имена точек;
<имя>_CENTER - NAIF ID тела, к которому соотносится;
<имя>_FRAME - прикреплена к системе координат, для Земли должно быть 'IAU_EARTH', но в примере 'EARTH_FIXED', потому что переопределена строками :
TKFRAME_EARTH_FIXED_RELATIVE = 'IAU_EARTH'
TKFRAME_EARTH_FIXED_SPEC = 'MATRIX'
TKFRAME_EARTH_FIXED_MATRIX = ( 1 0 0
0 1 0
0 0 1 )
в текстовом файле my_frames.tf, для того, чтобы при необходимости переопределить на более точную модель :
TKFRAME_EARTH_FIXED_RELATIVE = 'ITRF93'. (в файле my_frames.tf определены эти же имена для SPICE, но уже в качестве начал систем их топоцентрических координат.)
<имя>_IDCODE - свой уникальный код, начинающийся с кода тела прикрепления;
<имя>_XYZ    - позиция относительно центра Земли или если подключить таблицу геоидов Земли pck00009.tpc то в виде :
<имя>_LATLON = (широта долгота высота), градусы десятичные, высота в метрах;
<имя>_DXYZ    - ускорение метры в год, для неподвижных точек не указывается;
<имя>_EPOCH   - эпоха фиксации в TDB , для неподвижных точек не требуется;
<имя>_BOUNDS - границы времени действия;

     Программой pinpoint.exe преобразуем в бинарный файл, в командной строке (я ввел в pinpoint.bat и запускал уже его):
>> pinpoint -def my_sites.def -pck pck00008.tpc -spk my_sites.bsp
Полученный my_sites.bsp подключаем cspice_furnsh( 'путь\my_sites.bsp' );

my_frames.tf

set_observer.m

get_observer.m

     Для того, чтобы динамически изменять место и привязанную к нему топоцентрическую горизонтальную систему координат можно воспользоваться скриптами :

Топоцентрическая система координат

     Дает возможность вовлечь в расчеты местные направления: 'верх', 'север', 'запад', 'горизонт' и т.д.
Координаты пересчитываем: lon=-lon, lat=lat-90. Например, Астрахань: 48°03'00"E -> -48, 46°21'00"N -> -43.24; (высота не используется поскольку система координат это ориентирование осей, а высота над поверхностью Земли вводится в точку наблюдателя (см. выше)
Создадим файл фреймов my_frames.tf и запишем в него:

Доступ к переменным пула

В обычный текстовый файл поместим в текстовом виде переменные, которые зарузим cspice_furnsh( 'путь\my_stars.def' );

[values, found] = cspice_gcpool( name, start, room ); % доступ к srt
[values, found] = cspice_gdpool( name, start, room ); % доступ к DP
[values, found] = cspice_gipool( name, start, room );  % доступ к integer
[string, found] = cspice_stpool( item, nth, contin );        % доступ к string
где start = номер в массиве, room - кол-во возвращаемых данных. found=1, если найдена.

cspice_clpool;  % очистить пул данных
cspice_lmpool; % данные хранятся в массиве вместо внешнего файла.
cspice_pcpool; % вставки STR прямо в пул данных
cspice_pdpool; % вставки double precision ...
cspice_pipool;  % вставки integer ...

[kvars, found] = cspice_gnpool( name, start, room );  % загрузка в массив даным из пула ПО ИМЕНИ ...

count = cspice_ktotal( kind ) % количество загруженных файлов, kind: 'ALL', 'SPK', 'CK', 'PCK', 'EK', 'TEXT', 'META'

[ filtyp, source, handle, found] = cspice_kinfo( file) % получение информации о файле по его имени
[ file, filtyp, source, handle, found ] = cspice_kdata( which, 'spk'); % загрузка данных о SPK в пуле, загруженным № which

     cspice_pdpool('TKFRAME_DYN_TOPO_ANGLES',[-lon; -90 + lat; 180]); % вставка координат в топоцентрическую сист.коорд.

values1 = cspice_bodvcd( 399, 'RADII', 3); % получить радиусы Земли - 3 значения
values1 = cspice_bodvrd( 'EARTH', 'RADII', 3)

Получение широты и долготы наблюдателя

Внутреннее
представление координат

   1 вариант:  Из разницы углов местной и геоцентрической систем координат:
observer= 'N30_E40_TOPO'; % топоцентрическая система координат
et=0; % et=момент времени, cspice_dpr - градусов в радиане
[lat,z,lon]=cspice_m2eul( cspice_pxform('IAU_EARTH', observer, et ), 2,1,3) ;
lat = lat * cspice_dpr; % вн. представление широты наблюдателя
lon = lon * cspice_dpr; % вн/ представление долготы наблюдателя
if lat < - 90; sslat='S'; else sslat='N'; end % строка знака долготы
if lon > 0; sslon='W'; else sslon='E'; end % строка знака широты
pnlat = 90 - abs(lat); % абс/ значение широты наблюдателя
pnlon = 180 - abs(lon); % абс. значение долготы наблюдателя
fprintf( 'Место наблюдения : %s, %2.2f%s %2.2f%s %s%s%s %s%s%s\n', ...
   observer, pnlat, sslat, pnlon, sslon, '(', angl2minsec(pnlat), sslat, angl2minsec(pnlon), sslon,')');

Выдает :

Место наблюдения : N30_E40, 30.00N 40.00E (30°0'2.558e-011"N 39°59'60"E)

       2 вариант: Из параметров заданной системы координат.
[val, found] = cspice_gdpool( ''TKFRAME_',obsrvr, '_TOPO_ANGLES'', 2, 1 ); lat=90+val;
[val, found] = cspice_gdpool( ''TKFRAME_',obsrvr, '_TOPO_ANGLES'', 1, 1 ); lon= - val;

Изменение широты и долготы наблюдателя

     Для того, чтобы получать данные из разных широт или долгот наблюдателя с каким-то шагом пропишем в my_frames.tf топоцентрическую систему координат с именем, например, DYN и в my_sites.def точку на геоиде DYN, как для обычных. Затем из программы будем изменять широту и долготу:

observer = 'DYN'; frame = strcat(observer,'_TOPO'); % frame = 'DYN_TOPO'
lat=30; lon=40; alt=0.0; % водим параметры строкой
cspice_lmpool(strcat('TKFRAME_' , obsrvr , '_TOPO_ANGLES = (', num2str(-lon), ' , ', num2str(-90 + lat), ', 180.0000000000000)') );
cspice_lmpool(strcat('DYN_LATLON = (', num2str(lat), ' , ', num2str(lon), ' , ', num2str(alt), ' )' )); % загрузили их в соотв. объекты

cspice_pdpool(strcat( 'TKFRAME_' , obsrvr , '_TOPO_ANGLES') , [-lon; -90 + lat; 180]); % или числами

'Подключение' звезд.

   Нужную звезду можно идентифицировать по различным каталогам на портале Asteroseismic Modeling Portal. Уточнить положение звезды можно по, привязанным к координатам, снимкам телескопа Хаббл.

ppm_stars_vm3.def звезды до 3 величины из каталога PPM для загрузки в SPICE

ppm_stars_vm5.def звезды до 5 величины из каталога PPM для выборки

     Прямого доступа к оригинальным звездным каталогам в SPICE нет. Современные звездные каталоги ( ftp://naif.jpl.nasa.gov/pub/naif/generic_kernels/stars/ - 'Гиппарх', 'Тихо Браге', PPM) НАСА 'прицепило' функциями к SPICE еще в 1996 году, но только в версии ФОРТРАН. Я сконвертировал PPM каталог в текстовый kernel - загрузить (<= 3 зв.в.) для обычной загрузки в SPICE, весь каталог загрузить не получится из-за перегрузки имен в среде SPICE. В этот kernel можно добавить необходимые звезды из файла со звездами до 5 величины.

    Mожно получить данные любых других звезд из этих каталогов и рассчитать их положение в RA и DEC на нужный момент. Для начала нужно скачать звездный каталог. Я использовал Positions and Proper Motions Star Catalogue (скачать с NAIF) - каталог позиций и собственных движений звезд, потомок SAO Star Catalog эпохи J1950. Распаковал и утилитой >> tobin.exe ppm.xdb ppm.bdb , сконвертировал в совместимый формат.

     Программой starsb.exe по инструкции из звездного каталога генерирует набор звезд упрощенного типа для работы в SPICE: prompt> starsb 'имя_командного_файла'

     Определим квадрат в 6 ° RA-DEC, расположенный в центре небесного экватора, с самым западным RA в 357 °:
STARSB_WESTRA = 357.0
STARSB_EASTRA = 3.0
STARSB_STHDEC = -3.0
STARSB_NTHDEC = 3.0

     Определим 5 ° RA-DEC квадрат с самым западным RA 10 °, самым восточным RA 15 °, самым южным склонением в 45 °: STARSB_WESTRA = 10.0
STARSB_EASTRA = 15.0
STARSB_STHDEC = 45.0
STARSB_NTHDEC = 50.0

     Все RA [0, 360] и DEC [-90,90] в градусах относительно инерциальной системы J2000.

В интерактивной утилите Inspekt.exe наберем команду вызова пакетного файла 'Inspekt> start my_ppm.bat'.
Содержимое my_ppm.bat :

Получим текстовый каталог звезд на эпоху J2000.0, который можно импортировать в файл Excel (разделитель - tab). В программе можно теперь использовать параметры звезд введенные вручную.

Чтобы получить данные расширенного типа и досчитать положение движущихся звезд в программе, используем данные мастер-каталога ppm.bdb, экспортируя его (файл Excel - ppm_vm5.xlsb),

     'Неподвижные' звезды по степени уменьшения неподвижности.

PPM		DEC	RA	V	Примечание (координаты на J2000.0 взяты см.)
86047	61 Cygin			4.84	Самая подвижная звезда на небе J2000 RA 21 06 54.6 DEC +38 44 45
400076	Alnitak	-1,9425	85,18958	1,8	Левая звезда пояса Ориона.
348537		-56,7335	306,4117	1,9
60323		45,28032	310,3579	1,3
53742		49,31354	206,8878	1,9
327928		-46,9585	332,0558	1,7
148916		6,349893	81,28287	1,6
91373		23,46529	31,79035	2
12597		60,7168	14,17649	1,6
336856		-59,5098	125,6292	1,9
251718		-26,3933	107,0979	1,9
33769		55,95996	193,5047	1,8
357200		-69,719	138,3059	1,7
46127		49,86173	51,08014	1,8
17705		61,7522	165,9356	1,8
335149	Canopus	-52,6961	95,98746	-0,7
94361		28,61068	81,57259	1,6
312662		-47,3367	122,3833	1,8
296488		-37,1033	263,4022	1,6
269078		-26,2957	283,8161	2
122774		16,40003	99,42736	1,9
47925	Capella	46,00723	79,17065	0,1
297655		-34,3825	276,0436	1,8
341305		-59,6885	191,9313	1,3
362330		-69,0271	252,1655	1,9
130442	alf Boo	19,22192	213,9332	0
172188		9,874996	326,046	0,7
48617		44,94744	89,88346	1,9
81558	Vega	38,77732	279,2305	0	DEC +38 47 01.291, RA 18 36 56.3364, alf Lyr, 3 Lyr.
127140		11,96707	152,0971	1,3
341058		-57,1092	187,7909	1,6
97924		28,02715	116,3407	1,1
168779		8,860283	297,6868	0,8
72938		31,89072	113,6533	1,9
120061		16,51333	68,97905	0,8
192393		-8,65933	141,8971	2
323186		-42,9978	264,3296	1,9
251347		-28,9721	104,6564	1,5
217099		-17,9559	95,67497	2
149643	alf Ori	7,406855	88,79251	0,4	Semi-regular pulsating Star: alf Ori
337198		-54,7072	131,1751	2
331199	alf Eri	-57,236	24,42549	0,5	ACHERNAR
153068	Procyon	5,247656	114,838	0,4
209214		-17,9873	10,89306	2
274426		-29,6191	344,406	1,2
227262		-11,1607	201,299	1
175945		-1,20191	84,05334	1,7
360515		-60,3726	210,957	0,6
359410		-63,0988	186,6509	1,3
184482		-2,97547	34,83662	2
187839	Rigel	-8,20164	78,63445	0,1
217626	Sirius	-16,6867	101,2978	-1,5	DEC -16 42 58.017, RA 06 45 08.9173, Dog Star
431		89,26444	37,89261	2
360911	alf Cen A	-60,8545	220,0556	0
265579		-26,4315	247,3521	0,9
89441	Sirrah	29,09347	2,094481	2	FK5 1, α Andromedae
25054	Caph	59,15267722	2,281568343	2.3	FK5 1, Variable Star of delta Sct type, * bet Cas

Время.

Аналемма: фото солнца в одно и то же время дня в течении года - видимый полдень не 12:00 UTC. (взято)

    Время это 4-тое, более важное чем 'XYZ', измерение, потому что оно оказывает влияние на все три пространственных измерения. Какой смысл добиваться правильности расчета координат, при неправильно определенном времени?! Все расчеты SPICE ведутся в равномерном 'эфемеридном времени' (без скачковых секунд) - Ephemeris Time (ET), которое имеет внутреннее представление в секундах от эпохи 'J2000.0'. В это же представление перед вычислениями переводятся UTC, JD и другие системы представления времени.

    Аналемма солнца
    Изначально время дня определяли по солнцу, полднем (12:00) считался момент наивысшего положение солнца над горизонтом, когда азимут = 0°. Этот момент наступает неравномерно в равномерном времени UT. Суточное вращение Земли ~ равномерно, но из-за эксентриситета своей орбиты, Земля летит то быстрее к перигелию ~ 3 января, то медленее к афелию ~ 4 июля. За счет доворота или недоворота вокруг своей оси движением по орбите , видимый солнечный полдень, день ото дня, то чуть запаздывает к сидерическому звездному времени, то опережает его. День ото дня может быть короче до ~20 сек или длиннее до ~30 сек, а полдень наступает раньше на ~ 14 минут в '06.02.???? 12:00 UTC' и через 16 мин после '03.11.???? 12:00 UTC'. Если cфотографировать солнце в одно и то же время, например, в '12:00 UTC', то получится аналемма (греч. ανάλημμα, «основа, фундамент») от рисунка - поправки на солнечных часах (анимационные графики).

     LST (Local Solar Time) - местное солнечное время, передает видимую высоту солнца на небосводе. Местное солнечное время отличается от планетоцентрического и не привязано к UTC, на экваторе восход ~6:00:00, полдень = 12:00, закат ~18:00. час М.С.В., например, на Марсе час будет длится ~62 мин UTC.

[hr,min,sec,time,ampm] = cspice_et2lst (et, body, lon, type);

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

     UTC - универсальное координированное время, Coordinated Universal Time синхронизировано по частоте с атомным временем TAI (TDT), но уводится (координируется) 'скачковыми' секундами за 'звездной' полночью 00:00:00 UT. Из-за неравномерности вращения Земли равномерное UTC уходит от астрономического времени UT (обновленное GMT) секунды которого переменны, потому что они определенная часть суток и года, определенных по прохождению звезд.

    Чтобы разница между их 00:00:00 не накапливалась более 0.9 сек UTC изредка согласовывают скачковыми секундами, см. файл поправок времени naif0010.tls. Когда разница накапливается до 0.7 с, IERS объявляет в IERS Bulletin C (Leap Second Announcements) о поправке, что будет внесена после последнего 'нормального' UTC времени 31.12 или 30.06 (см. табл.)

    UTC       GPS               TAI                    TT (TDT)
 ----+---------+-----------------+------------------------+-----
     |         |                 |                     ET 1984.0
     |<--34s Leap Seconds------->|<----32.184S fixed----->|
     |         |                 |                        |
     |<--15s-->|<----19s-------->|                        |
     |              fixed                                 |
     |                                                    |
                                                          |
 -----+----- DeltaT = 32.184 s + (TAI-UTC) - (UT1-UTC) ---+-----
     UT1 (UT)  66.65750 s          34 s     -0.47350 s TT (TDT)

UTC удерживается ±0.9s относительно UT1 'скачковыми' секундами
Разница между GPS  и TAI и TDT - константа.

    UT - усредненное солнечное время, неравномерное, без скачковых секунд. Истинные солнечные полдни усредняются за год. Секунды, часть 24 часов неравномерных суток.
     UT0 (GMT) - устанавливалось в одной обсерватории (Greenwich), поэтому включало эффект 'гуляния' ~10...20 м в 400-500 дней земной оси. UT0 = UT1 + tan(lat) * (x * sin(long) + y * cos(long)), где x и y смещения полюса в " опубликованные в бюллетене А, IERS и должны быть конвертированы во время 1 сек = 15 ", а lat и long координаты наблюдателя.
    UT1 - устанавливается многими обсерваториями, поэтому эффект 'гуляния' полюсов исключен.
    UT2 - UT1 поправленное на сезонные вариации вращения Земли. Сегодня, говоря UT, обычно имеют ввиду UT2.
UT2 = UT1 + 0.022 * sin(2*π*t) - 0.012 * cos(2*π*t) - 0.006 * sin(4*π*t) + 0.007 * cos(4*π*t), где t = 2000.0 + (MJD - 51544.03) / 365.2422
     GMST - Средне звездное время по Гринвичу (Greenwich Mean Sidereal Time) sec GMST = 365.25/366.25 sec UT1 second. День начинается в точке В.Р.Д. на Г.М. Час GMST = углу от средней позиции точки ♈, В.Р.Д., пренебрегая нутацией.
Это может показаться странным, но GMST может быть связано с UT1 по формуле GMST (в сек на UT1=0) = 24110.54841 + 8640184.812866 * T + 0.093104 * T^2 - 0.0000062 * T^3, где T юлианские века на J2000.0 UT1, T = d / 36525, d = JD - 2451545.0
    LMST - Local Mean Sidereal Time. LMST = GMST + (долгота наблюдателя/(360°/24), в часах)
     LST - Local Sidereal Time. Hour Angle = LST - Right Ascension
    GAST - Greenwich Apparent Sidereal Time. GAST = GMST + (уравнение равноденствий (времени))

    При расчетах в UTC на будущее невозможно предсказать, когда будет внесена поправка, примерно до следующей поправки время может соответствовать UTC, а после предполагаемой поправки нет. Время UTC с 58.999 ... по 59.999 ... секунду в момент поправки UTC на секунду назад не существующее в природе. А с 59.999 ... по 60.999 ... секунду UTC в период поправки на секунду вперед существует, например, '2012-06-30 23:59:60 UTC'

     TAI (атомное время) и TDT синхронны, TDT привязано к TAI, но TDT - TAI = 32.184 сек. До 1985 года JPL's Orbit Determination Program (ODP) и другие программы, CRS записи, использовали другое значение 32.1843817 сек., его можно изменить в файле поправок времени naif0010.tls. Так же JPL'S Optical Navigation Program (ONP) не использует периодическое изменение (K sin E) между TDB-TDT, установка K =0 решает проблему.

     ET - 'Эфемеридное время' (Ephemeris Time) равномерная шкала, общее название TDB и TDT. 
  TDB - барицентрическое динамическое время (Barycentric Dynamical Time) используется для описания движения тел, относительно барицентра Солнечной системы. TDT - земное динамическое время (Terrestrial Dynamical Time) используется при описании движения объектов около Земли.

     TDB и TDT подобны и равны по внутренней частоте, но не синхронны, по разному в разное время года, из-за релявисткого эффекта, т.к. движутся относительно друг друга. TDB - TDT = K * sin (E), где К константа, а E аномалия экцентриситета земной орбиты. TDB - TDT = ~ 0.000030 сек (при исключении малопериодных флуктуаций), т.е. TDB - TDT = периодическая функция с величиной ~ 0.001658 сек и периодом = сидерическому году. E = M + EB sin (M), где EB эксцнтриситет, а M (mean) средняя аномалия гелиоцентрической орбиты барицентра Земля-Луна. M = M0 + M1*t, где t эпоха TDB, в сек от J2000. Константы хранятся в пуле: DELTET/K = 1.657D-3, DELTET/EB = 1.671D-2, DELTET/M = ( 6.239996D0 1.99096871D-7 )

     TAI, TDT, TBT - равномерная секундная ± от J2000.0 шкала, никакого отношения к годам, месяца, дням, часам не имеющая. Только, ради удобства восприятия и сравнения с планетным земным временем - UTC, может вводиться и выводиться, 15.10.1582 12:00:00.000, подобно ему. Условно считается, что 'дне ET' 86400 секунд, хотя начало его суток 00:00 ET за соответствующий огромный период сместится по суткам шкалы времени UTC , полдень, которой скачковыми секундами, вставляемыми через месяцы и годы, привязывается к истинному солнечному полдню места.

    Julian Date - система юлианских дней - равномерная шкала, подобно 'ефемеридному времени', что без скачковых секунд, не имеет скачковых дней, что позволяет вычислять количество дней между датами. Существует 2 типа JD - Julian Ephemeris Date (JED) и Julian Date UTC (JDUTC). JED = cspice_j2000 + et/cspice_spd; Дробная часть JDUTC равна 0 в '12:00:00 UTC'. Префикс JD зависит от контекста, лучше использовать JDUTC, JDTDB, или JDTDT.

Формат времени

     Время вводится строкой через функцию :

ephemeris_time=cspice_str2et('1582-10-04 12:00 ');% по умолч. UTC-0

et1 = cspice_str2et ( '01-30-0100 BC' ); % до нашей эры
et2 = cspice_str2et ( '01-30-0100 AD' ); % после нашей эры

et=cspice_str2et('1582-10-04 12:00 UTC+2');% с местным часовым поясом
et=cspice_str2et('1582-10-04 12:00 UTC+2:30');% с местным часовым

   01.01.2000 г. 12:00 ТDT
= 01.01.2000 г. 11:58:55.816 UTC
= 01.01.2000 г. 11:59:27.816 TAI   (International Atomic Time)
TDT = TAI + 32.184; % Начало в 01.01.1958 00:00:00 Гринвич
JD 2451545.0            % юлианский день

Вывод ET по маске маркеров формата:
calendar_date = cspice_timout( ephemeris_time , 'DD.MM.YYYY ERA HR:MN:SC.#### ::RND::MCAL::UTC' );
, где ERA нужно для вывода эры B.C. или A.D., когда с минусом в формат не умещаются годы, например, -2999 даст ****
::MCAL - юлианско-григореанский смешанный календарь, т.е. ephemeris_time + cspiсe_spd (секунд в день) выдаст уже '15.10.1582 12:00:00.000', а не юлианский день '05.10.1582 12:00'
1995-18T % 18 день 1995 года, время 00:00

Разница = TDB-UTC сек. в 1965-2020 гг Регулирование UTC началось в 1972 г. С момента сост. графика (2012 г.) в будущее UTC тоже не регулируемо.

Конвертация между различными временными шкалами [unitim] = cspice_unitim( epoch, insys, outsys ) (см. табл.) :

Разница между TDB и UTC
delta = cspice_deltet( epoch, eptype ), где eptype 'UTC' или 'ET' TDB. Практически на графике отображено содержимое файла 'скачковых' секунд naif0010.tls
До 01.07.1972г. UTC = TDB + 41.183942417938496 сек.

Конвертер дат (SSD JPL NASA)
Конвертер даты в юлианский день (французское косм. институт).
Конвертер даты <-> юлианский день USNO ->

     SCLK - Часы космического корабля - Spacecraft Clock (SCLK) служат для внутренней координации и проставки штампов времени на записях (фото, датчики и т.п.) Для станции Galileo p/rrrrrrrr:mm:t:e, где p - номер части временной шкалы, r - rim counts, m - minor frame, t - real time interrupt, e - mod eight count. Для перевода внутреннего времени в UTC или ET нужно загрузить функцией cspice_furnsh файл коэффициентов - 'spacecraft clock coefficients'.

Григорианская реформа календаря 1582 г.

     В 1582 году по указу папы Григория XIII календарь Западной Католической церкви совершил скачок, оторвавшись от исходного юлианского календаря (старого стиля) на -10 дней. Это была знаменитая григорианская реформа календаря 1582 года. Бало еще три скачка по -1 дню - в 1700, 1800 и 1900 годах, после чего разница между старым и новым стилем достигла своего современного значения 13 дней. В разных странах календарь изменяли в разное время, см. неполную таблицу.

 • Полезные ресурсы по теме 'время' ->

Работа с функциями окон времени.

    В SPICE массив расмерностью 2n Х 1 т.е. [1;2;   3;4 ...   9;8] это n окон времени, с n= начало n+1 = конец интервала.

st_win=(reshape(st_win,[],2)); % преобразование последовательных окон времени в 2 колонки 1) начало 2) конец
st_win= st_win (:,1); % отбор 1-ой колонки

Абберация и поправка на скорость света (разница между видимым и реальным положением).

     Поправки на абберацию (смещение видимого положения объекта из-за призмы планетарного воздуха) и скорость света необходимы для того, чтобы разделить понятия видимого положение объекта и его реальных координат между небесными телами в прошлом и настоящем. Атмосфера приподнимает объекты над горизонтом, свет доносит до нас прошлые координаты объекта. Чтобы выяснить где на самом деле физически, а не визуально, находится объект учитывают эти поправки. В SPICE учтены следующие виды поправок:

Ни специальные ни общие эфекты относительности ни учитываются.(?)

Расчет координат планет и других космических тел.

     Из скрипта MatLab:

программа расчета координат планет и др. тел солнечной системы

target = 'Mars'; % Расчитываем координаты Марса, или любого другого небесного тела, на которое есть эфемериды

et = cspice_str2et( '2000-01-01, 12:00 TDT' ); % момент времени. TDT - земное динамическое время
TFMT='DD.MM.YYYY HR:MN:SC ERA::RND::TDT'; TFMTH='HR:MN:SC.## ::RND::TDT' ; % вывод в том же TDT

abcorr = 'LT+S'; % видимые положения планет, поправка на скорость света и аберацию

frame='KURSAVKA_TOPO'; % ГОРИЗОНТАЛЬНАЯ сетка координат
observer='KURSAVKA'; % 44.46 N 42.51 E
[pos, lt] = cspice_spkpos ( target, et, frame, abcorr, observer); % поз. планеты в картезианских коорд.
[r, lon, lat] = cspice_reclat(pos); % картезианские координаты в высотные
r = cspice_convrt( r, 'km' , 'au'); % Удаление в астрономических единицах
lon =cspice_convrt(lon,'radians','degrees'); az=-(-180+lon); % Азимут из долготы
lat=cspice_convrt(lat,'radians','degrees'); % Высота над горизонтом

frame='J2000'; % ЭКВАТОРИАЛЬНАЯ сетка координат
[pos, lt] = cspice_spkpos ( target, et, frame, abcorr, observer); % поз. планеты в картезианских коорд.
[r, ra, lat] = cspice_recrad(pos); % картезианские координаты в радиальные небесной сферы
r = cspice_convrt( r, 'km' , 'au'); % картезианские координаты в высотные
ra =cspice_convrt(ra,'radians','degrees'); lat=cspice_convrt(lat,'radians','degrees');
rah=cspice_timout((ra/15) * (cspice_spd/24)+ (cspice_spd/2) , TFMTH); % прямое восхождение в часах

frame='ECLIPJ2000'; % ЭКЛИПТИКАЛЬНАЯ сетка координат
[pos, lt] = cspice_spkpos ( target, et, frame, abcorr, observer); % позиция планеты
[r, ra, lat] = cspice_recrad(pos); % картезианские координаты в радиальные небесной сферы
r = cspice_convrt( r, 'km' , 'au'); ra =cspice_convrt(ra,'radians','degrees'); lat=cspice_convrt(lat,'radians','degrees');
disp(strcat( sprintf( 'Долгота, ° : % 12.8f\n', ra ), ' (', angl2minsec(ra), ')')); % пример вывода
disp(strcat( sprintf( 'Широта, ° : % 12.8f\n', lat ), ' (', angl2minsec(lat), ')'));

cspice_kclear

Угол между планетами

   Программа:

%Расчет окон времени из периода, когда угол между Солнцем и Луной > 10.6 °
MAXWIN = 1000; TIMFMT = 'DD.MM.YYYY HR:MN:SC.###### (TDB) ::TDB ::RND';
cspice_furnsh( 'd:\Documents and Settings\user\Мои документы\MATLAB\astro\mice\data\standard.tm' );
et = cspice_str2et( { '2012-01-01', '2012-03-01'} ); cnfine = cspice_wninsd( et(1), et(2) );
step = 6.*cspice_spd; % поиск с шагом 6 дней (в секундах) ???
adjust = 0.; refval =10.6 * cspice_rpd;
moon_10_sun = cspice_gfsep( 'MOON', 'POINT', 'NULL', 'SUN', 'POINT', 'NULL', 'LT+S', '399201', '>', refval, adjust, step, MAXWIN, cnfine );
cspice_kclear

Вывод:

Восходы солнца

    Программа :

Выдает:

Равноденствия и солнцестояния, длина сезонов

   В момент истинного равноденствие угол между наклоном земной оси и Солнцем = 90° (солнце слева, Земля движется 'вперед') и день везде на Земле 'равен ночи' (+- полсуток/360 по поверхности, от момента астрономического равноденствия). В момент летнего солнцестояния угол между Солнцем и наклоном земной оси = 0° (Земля наклонена к Солнцу). Принято считать, что моменты креста Солнца = 0°, 90°, 180°, 270° долготы Солнца по зодиаку.

В (Альмагест) описано, что во времена Гиппарха, (БСЭ) около 180 - 190 до н. э., Никея, - 125 до н. э., Родос, от В.Р.Д. до Л.С.С 94.5 дней, от Л.С.С. до О.Р.Д. 92.5 дня и что во времена Птоломея (БСЭ) 2 в.н.э. эти значения практически не изменились, из чего Птоломей сделал вывод, что длины сезонов неизменны. Было бы интересно по длинам сезонов рассчитать время Гиппарха и Птоломея (в конце мая выйдет новый пакет SPICE с которым это станет возможно)

equinoxes_solstices_seasons.m

Сегодня это такие даты креста солнца:

Расчет астрономических моментов новолуний

     Программа:

cspice_furnsh( 'd:\Documents and Settings\user\Мои документы\MATLAB\astro\mice\data\standard.tm' ); % эфемериды DE406
MAXWIN = 20000; MAXIVL = MAXWIN / 2; TIMFMT = 'DD.MM.YYYY HR:MN:SC ::RND::MCAL::UTC'; TIMLEN = 35;
st_win = cspice_wninsd ( cspice_str2et ( '1970-01-01 UTC' ), cspice_str2et ( '1971-05-01 UTC' ) ); % окно времени поиска

Эксцентриситет лунной орбиты

step = cspice_spd*14 ; % 14 дней менее 1/2 цикла
adjust = 0.0; refval = 0; abcorr = 'LT'; % учитывать скорость света
% угловое расстояние между центрами двух тел, наблюдатель на геоиде
observer='DYN' %топографическая точка на геоиде
res_sep = cspice_gfsep('SUN','POINT','NULL','MOON','POINT','NULL',abcorr ,observer,'LOCMIN',refval, adjust,step,MAXWIN,st_win);
% расстояние между 2-мя телам. Наблюдатель на Луне
observer='MOON';
res_dist = cspice_gfdist( 'SUN', abcorr, observer , 'LOCMIN', refval, adjust, step, MAXWIN, st_win);
fprintf ( '\n\n\nМежду луной и солнцем минимальные: \n' );
    for i=1:numel(res_sep)/2;  fprintf ( 'Угол: %s, Расстояние: %s\n', cspice_timout( res_sep(i*2-1), TIMFMT ), ...
       cspice_timout( res_dist(i*2-1), TIMFMT ));  end; cspice_kclear;

Тогда, когда угловое расстояние между солнцем и луной = 0 (наблюдатель в центре Земли) наступает момент астрономического новолуния, расстояние между Солнцем и Луной минимальные, это должно произойти примерно в одно время, а здесь разница часы и сутки! Из-за эксцентриситета лунной орбиты.

     На сайте http://www.moshier.net/ опубликована таблица новолуний, рассчитанная по этим же эфемеридам. DE406. Сравним с расчетом по моему скрипту на период григореанской реформы (формат времени TDT)

     Разница, видимо, оттого, что в другом расчете новолуние топоцентрицеское с того бока земли с которого Луна и Солнце уже соеденились, до того как это будет видно с "центра Земли".

Расчет видимых новолуний и старолуний

     Программа (скачать):

% (I.Bromberg.AHC) приводит из (CambridgePress.СС3) условия нового месяца, адапрированные им к Иерусалиму:
% 1. Центр солнечного диска > 4.5° под горизонтом;
% 2. Высота луны над горизонтом > 4.1°;
% 3. Угол зрения между центрами солнца и луны > 10.6° (геоцентрический);

clear all % очищаем среду МатЛаб от предыдущих, возможных, остатков в памяти после прерываний
cspice_furnsh( 'd:\Documents and Settings\user\Мои документы\MATLAB\astro\mice\data\standard.tm' );

MAXWIN = 200000; MAXIVL = MAXWIN / 2; TIMFMT = 'DD.MM.YYYY HR:MN ::RND::UTC+3'; TIMLEN = 35;
result = cspice_wninsd ( cspice_str2et ( '2012-01-01' ), cspice_str2et ( '2012-05-01' ) );
abcorr = 'LT+S'; % видимые положения планет

% интервал поиска нужно "вручную" выставить так, чтобы старолуние попало в интервал, а заканчивался он новолунием
et1=cspice_str2et ( '1972-04-01' ); et2=cspice_str2et ( '1972-05-01' ); st_win = cspice_wninsd (et1 ,et2 );

observer = 'DYN'; frame = strcat(observer,'_TOPO'); crdsys = 'LATITUDINAL'; coord = 'LATITUDE';
lat=30; lon=40; alt=0.0; % вводим координаты наблюдателя
itt=2; % иттераций
result_sum_new=cell(1,itt); result_sum_old=cell(1,itt); % массивы для результатов

% для размера массива результатов "астрономические новолуния" минимумы угла между луной и солнцем, наблюдатель в горизонтальной местной системе координат
step = cspice_spd*14;
moon_o_n = cspice_gfsep('Moon','POINT','NULL','Sun','POINT','NULL','LT', observer, 'LOCMIN',0,0,step,MAXWIN,st_win);
a=zeros(cspice_wncard( moon_o_n ),1); moon_o_n=horzcat((1:numel(moon_o_n)/2)',a,a, moon_o_n(1:2:numel(moon_o_n)),a,a,a,a,a);
moon_o_n_sum = zeros(itt*a,11); moon_o_n_surf= zeros(0); % массивы для итогов

for f=1:itt; % иттерация широты
% загрузка координат наблюдателя в систему координат
cspice_pdpool('TKFRAME_DYN_TOPO_ANGLES',[-lon; -90 + lat; 180]);
cspice_lmpool(strcat('DYN_LATLON = (', num2str(lat), ',', num2str(lon), ',', num2str(alt), ' )' ));

result = st_win; step = cspice_spd/2 ; refval = 10.6 * cspice_rpd; % Луна с солнцем не ближе 10.6 градусов
result = cspice_gfsep('MOON','POINT','NULL','SUN','POINT','NULL','LT+S',observer,'>',refval,adjust,step,MAXWIN,result);

step = cspice_spd/2 ; relate = '<'; refval = -4.5 * cspice_rpd; % Солнце менее 4.5 градусов под горизонтом
result = cspice_gfposc ( 'SUN', frame, abcorr, observer, crdsys, coord, relate, refval, adjust, step, MAXIVL, result);

step = cspice_spd*14 ; relate = '>'; refval = 4.1 * cspice_rpd; % Луна > 4.1 градуса над горизонтом
result = cspice_gfposc ( 'MOON', frame, abcorr, observer, crdsys,coord,relate,refval,adjust,step,MAXIVL,result);

% обработка результата для этой широты, к а.н. привязать старолуния и новолуния и записать
jj=1; n = cspice_wncard( result )-1;
for ii=1:n;
[ st_1, fin_1 ] = cspice_wnfetd ( result, ii ); % output1 = cspice_timout( [st_1,fin_1], TIMFMT );
[ st_2, fin_2 ] = cspice_wnfetd ( result, ii+1 ); % output2 = cspice_timout( [st_2,fin_2], TIMFMT );
if (( moon_o_n(jj,4) - st_1) > 0) ... % is end of old moon before conj. (конец старолуния перед а. новолунием) ?
&& (( moon_o_n(jj,4) - st_2) < 0) % is start new moon after conj. (начало новолуния после а. новолуния) ?
moon_o_n(jj,2)= st_1; moon_o_n(jj,3)= fin_1; moon_o_n(jj,5)= st_2; moon_o_n(jj,6)= fin_2;
moon_o_n(jj,7)= abs((moon_o_n(jj,3)- moon_o_n(jj,2))/60);
moon_o_n(jj,8)= abs((moon_o_n(jj,6)- moon_o_n(jj,5))/60); % продолжительность видимого новолуния
moon_o_n(jj,9)= abs((moon_o_n(jj,6)- moon_o_n(jj,2))/3600); % не видно луны между с.л. и н.л.
moon_o_n(jj,10)=lat; moon_o_n(jj,11)=lon; % широта и долгота
jj=jj+1; % следующее а. новолуние

if jj>numel(moon_o_n(:,1)); moon_o_n_sum = vertcat(moon_o_n_sum, moon_o_n); break ;end % аккумулирование данных в общий массив
end; end; lat=lat+1; step = cspice_spd*14;
moon_o_n = cspice_gfsep('Moon','POINT','NULL','Sun','POINT','NULL','LT', observer, 'LOCMIN',0,0,step,MAXWIN,st_win);
a=zeros(cspice_wncard( moon_o_n ),1);
moon_o_n=horzcat((1:numel(moon_o_n)/2)',a,a, moon_o_n(1:2:numel(moon_o_n)),a,a,a,a,a);
end

cspice_kclear

Не видно луны между старым и новым месяцем с 0° по 57° С.Ш. (часов) 1970-1989 годы

  Часы невидимости луны
между старолунием и новолунием (воскресенья)
от 0° до 57° северной широты 1970-1989 гг. (в часах)

   Расчетные условия определения новолуния (аналогично я определял старолуния) взяты из (I.Bromberg.AHC) и (CambridgePress.СС3) :

1. Центр солнечного диска > 4.5° под горизонтом;
2. Высота луны над горизонтом > 4.1°;
3. Угол зрения между центрами солнца и луны > 10.6° (геоцентрический);
(адаптировано к Иерусалиму)

Видно, что выше 50° широты лунный календарь становится нестабилен и его невозможно использовать.
Расчетные условия из (I.Bromberg.AHC). (скачать программу и данные формата MatLab и Excel)

Не видно луны между соединением и новолунием

Видно 19-летний лунный цикл. Когда луна на закате солнца выше эклиптики, она становится видна скорее после коньюнкции.

     Интересны средние значения, с 47° широты резкий скачок минимальных значений больше 2 суток и чуть выше средние значения - более 4 дней ... еще выше, выше широты Москвы лунный календарь не имеет смысла ... не потому ли большинство древних столиц примерно на этой широте - "разделяй и властвуй" - две системы календарей лунная и солнечная ...

Яркость Луны (фаза), Солнца и освещенность неба

 Израильское Общество Новой Луны (Israeli New Moon Society), возглавляемое доктором Роем Хоффманом, публикует на своем сайте, таблицу видимости нового месяца. Я проанализировал эту таблицу и нашел, что есть много наблюдений, где солнце, вовремя обнаружения нового месяца, над горизонтом. По условиям CC3 это невозможно - солнце должно быть под горизонтом на 4.6°. Свой вопрос я направил доктору Рою.

Ограничение видимости звезд из-за освещенности неба - Limiting Magnitude Calculations.

Этот скрипт я содрал с k3pgp.org/star.htm и перевел его в Matlab + SPICE. В скрипте был свой расчет положений планет, я заменил его на SPICE, больше никаких изменений не делал (правильность не гарантирую, возможно, из-за ручной конвертации могла закрастся ошибка - сверяйте с оригиналом)

Выходные параметры :
Процент освещенной поверхности Луны, % - ill_frac
Коеффициент затухания [V mag per airmass], - K
Затухание [V magnitudes], - D
Магнитуда наблюдаемой луны, V Mag, - magmn
Магнитуда наблюдаемого солнца,V Mag, magsn
Общая яркость неба в точке, NanoLamberts - Bnl
Общая яркость неба в точке, V Mag per arcsec2 - Bmag
Яркость неба в точке от света луны, NanoLamberts - Bmnl
Яркость неба в точке от света луны, V Mag per arcsec2 - Bmmag
Ограничение на видимость звезд, V Mag - limmag
+/- ошибка ограничения на видимость звезд - magerr
Удаление луны, ° rhomn
Удаление солнца, ° - rhosn

Другие примеры:

   Поворот 3х3 матриц одной системы координат к другой.

Преобразование координат на поверхности Земли в относительно J2000 + ускорения
% Определяем геодезические координаты: долгота, широта, высота (Астрахань)
lon = 48 * cspice_rpd; lat = 46.33 * cspice_rpd; alt = -0.002;
utc = 'January 1, 1990'; et = cspice_str2et( utc ); % задаем время
abc = cspice_bodvrd( 'EARTH', 'RADII', 3 ); equatr = abc(1); polar = abc(3); % получаем ширину Земли в экваторе и в полюсах
f = ( equatr - polar ) / equatr; % фактор сплющивания Земли
estate = cspice_georec( lon, lat, alt, equatr, f); % прямогольная координата на Земле
estate = [ estate; [0.; 0.; 0.] ]; % расширяем массив для значений ускорения
xform = cspice_sxform( 'IAU_EARTH', 'J2000', et ); % получаем матрицу преобр. из "IAU_EARTH" в "J2000" в момент 'et'.
jstate = xform * estate;
jstate(1:3) ); % прямоугольные координаты в J2000
jstate(4:6) ); % прямоугольные ускорения в J2000

Получение координат Марса как его видно из Астрахани, относительно системы координат J2000
cspice_furnsh( 'd:\Documents and Settings\user\Мои документы\MATLAB\astro\mice\data\standard.tm' ); % загрузим эфемериды
target = 'Mars'; epoch = '2000 JAN 01, 12:00 UTC'; frame = 'J2000'; abcorr = 'LT+S'; observer = '399201'; % загрузим данные
et = cspice_str2et( epoch ); % преобразуем дату в эфемеридное время. Дата - начало J2000
[ state, ltime ] = cspice_spkezr( target, et, frame, abcorr, observer); state=state(1:3,1); % расчет состояния планеты Марс
[radius, ra, dec] = cspice_recrad( state ); % конвертирует прямоугольные координаты в RA и Dec
radius = cspice_convrt( radius, 'km' , 'au'); ra = cspice_convrt( ra, 'radians' , 'degrees'); dec = cspice_convrt( dec, 'radians' , 'degrees');
utc_epoch = cspice_et2utc( et, 'C', 3 );

disp('Координата из функции cspice_spkezr( target, et, frame, abcorr, observer);')
txt = sprintf( 'Положение планеты : %s', target); disp( txt )
disp(strcat(sprintf( 'Как он наблюдается с : %s', observer ), ' (Астрахань - 399201)'));
txt = sprintf( 'Относительног с. к. : %s', frame ); disp( txt )
txt = sprintf( 'В момент : %s', epoch ); disp( txt )
txt = sprintf( ' : i.e. %s', utc_epoch ); disp( txt )
txt = sprintf( ['Удаление,А.Е.(recrad): %12.16f'], radius ); disp( txt )
disp(strcat( sprintf( ['Прямое в.,° (recrad): %12.16f'], ra ), ' (', angl2minsec(ra), ')'));
disp(strcat( sprintf( ['Склонение,° (recrad): %12.16f'], dec ), ' (', angl2minsec(dec), ')'));

[radius, lon, lat] = cspice_reclat(state);
lon = cspice_convrt( lon, 'radians' , 'degrees'); lat = cspice_convrt( lat, 'radians' , 'degrees');
disp(strcat( sprintf( ['Долгота,° (reclat): % 12.8f'], lon ), ' (', angl2minsec(lon), ')'));
disp(strcat( sprintf( ['Широта, ° (reclat): % 12.8f'], lat ), ' (', angl2minsec(lat), ')'));

% [radius, lon, lat] = cspice_reclat(state);
% % Convert the cylindrical to rectangular.
% [rectan] = cspice_latrec(radius, lon, lat);
% state=(rectan-state) ./ state; % Calculate the relative error against the original position
% abc = cspice_bodvrd( 'EARTH', 'RADII', 3 ); equatr = abc(1); polar = abc(3); % получаем ширину Земли в экваторе и в полюсах
% f = ( equatr - polar ) / equatr; % фактор сплющивания Земли
% [ lon, lat, alt ] = cspice_recpgr( '399', rectan, equatr, f);
% alt = cspice_convrt( alt, 'km' , 'au');
% lon = cspice_convrt( lon, 'radians' , 'degrees'); lat = cspice_convrt( lat, 'radians' , 'degrees');
% disp(strcat( sprintf( ['Удаление, А.Е. : % 12.8f'], alt )));
% disp(strcat( sprintf( ['Долгота,° : % 12.8f'], lon ), ' (', angl2minsec(lon), ')'));
% disp(strcat( sprintf( ['Широта, ° : % 12.8f'], lat ), ' (', angl2minsec(lat), ')'));

[r, lat, lon] = cspice_recsph(state);
r = cspice_convrt( r, 'km' , 'au');
lon = cspice_convrt( lon, 'radians' , 'degrees'); lat = cspice_convrt( lat, 'radians' , 'degrees');
disp(sprintf( ['Удаление,А.Е.(recsph): %12.16f'], r));
disp(strcat( sprintf( ['Долгота,° (recsph): % 12.8f'], lon ), ' (', angl2minsec(lon), ')'));
disp(strcat( sprintf( ['Широта, ° (recsph): % 12.8f'], lat ), ' (', angl2minsec(lat), ')'));
cspice_kclear

Выдает:

Расшифровка дат древних зодиаков

   Тестовый пример Курсавского гороскопа 'на сегодня'.
    Пример программы для поиска даты зодиака, составленного в Курсавке 25.04.2012 2:30 UTC+3, при поиске с 1970 г. по 2020 г. нашел событие +- 1.5 дня, 23.04.2012 07:27 ... 26.04.2012 06:42 UTC. Такую 'неточность' можно объяснить тем, что при составлении гороскопа долготы планет огрублялись до созвездия. Система координат ECLIPJ2000 для расчета древних зодиаков не подходит из-за прецессии, достигающей больше градуса в столетие, нужно использовать FK4.

Круглый Дендерский древний зодиак (Египет)

Разбивка зон зодиака в созвездиях

№ Знак Границы зодиака,° Границы* в созвездиях,° Границы** в созвездиях,° 1 ♈ Овен ARI 0-30 31-56 26-51 2 ♉ Телец TAU 30-60 56-92 51-89 3 ♊ Близнецы GEM 60-90 92-118 89-118 4 ♋ Рак CNC 90-120 118-137 118-142 5 ♌ Лев LEO 120-150 137-172 142-174 6 ♍ Дева VIR 150-180 172-215 174-215 7 ♎ Весы LIB 180-210 215-236 215-236 8 ♏ Скорпион SCO 210-240 236-266 236-267 9 ♐ Стрелец SGR 240-270 266-296 267-301 10 ♑ Козерог CAP 270-300 296-326 301-329 11 ♒ Водолей AQR 300-330 326-349 329-346 12 ♓ Рыбы PSC 330-360 349-31 346-26 * ФиН исп. сейчас см. Египетские, русские итал. зодиаки. Приложение ** ФиН использовали с 2002 года по 17.11.2007 г.

Подробно о его расшифровке в книге "Новая хронология Египта" ->

      Имеем разбивку зон зодиака в созвездиях (см. табл.). Имеем уже подготовленную ФиН расшифровку 'Дендерского круглого' зодиака:

1. Луна в Весах
2. Солнце в Рыбах
3. Меркурий в Водолее или Рыбах
4. Сатурн в Деве или Весах
5. Марс в Козероге
6. Венера в Овне или Рыбах
7. Юпитер в Раке или Близнецах

Моя программа для SPICE в MatLab.

cspice_furnsh( 'd:\Documents and Settings\user\Мои документы\MATLAB\astro\mice\data\standard.tm' );
TIMFMT = 'DD.MM.YYYY HR:MN ::RND::MCAL::UTC'; TIMLEN = 35; % выводимый формат времени. Внимание, UTC+?!
et1=cspice_str2et ( '0000-01-01' ); et2=cspice_str2et ( '2008-01-01' ); res = cspice_wninsd (et1 ,et2 ); % исследуемый интервал
abcr = 'LT+S'; % учитывать оберацию и скорость света, хотя это излишне
obsrvr='CAIRO' ; % место записи/наблюдения Каир
crds = 'CYLINDRICAL' ; crd = 'LONGITUDE' ; adj = 0.0; st = cspice_spd; MAXWIN = 20000; MAXIVL = MAXWIN / 2;

по старой разбивке (в программе по новой) :

... текст программы ...

% решение/output для ECLIPJ2000 т.е. без учета прецессии:
% с 07.03.-963 08:01 по 09.03.-963 11:46 UTC
% с 07.03.-583 22:59 по 08.03.-583 12:12 UTC
% с 22.03. 331 22:39 по 22.03. 331 22:39 UTC
% с 31.03.1185 02:49 по 31.03.1185 02:49 UTC

% решение/output для FK4 по старой разбивке (более ранние даты здесь не представлены):
% с 07.03.-963 12:35 UTC (14:50:57 LST) по 09.03.-963 15:14 UTC (17:30:54 LST)
% с 07.03.-583 19:06 UTC (21:19:58 LST) по 08.03.-583 15:48 UTC (18:02:29 LST)
% с 18.03.-583 02:03 UTC (04:20:16 LST) по 18.03.-583 02:03 UTC (04:20:16 LST)
% с 13.03. 568 09:04 UTC (11:11:56 LST) по 13.03. 568 09:04 UTC (11:11:56 LST)

Результат по новой разбивке:
с 16.02.2969 B.C. 21:06 UTC (23:34:37 LST) по 16.02.2969 B.C. 21:06 UTC (23:34:37 LST)
с 24.02.2115 B.C. 19:37 UTC (21:59:34 LST) по 24.02.2115 B.C. 19:37 UTC (21:59:34 LST)
с 27.02.1438 B.C. 20:20 UTC (22:37:47 LST) по 02.03.1438 B.C. 08:37 UTC (10:55:02 LST)
с 10.03.1438 B.C. 11:33 UTC (13:52:45 LST) по 10.03.1438 B.C. 11:33 UTC (13:52:45 LST)
с 07.03. 964 B.C. 12:35 UTC (14:50:57 LST) по 09.03. 964 B.C. 15:14 UTC (17:30:54 LST)
с 07.03. 584 B.C. 19:06 UTC (21:19:58 LST) по 08.03. 584 B.C. 15:48 UTC (18:02:29 LST)
с 18.03. 584 B.C. 02:03 UTC (04:20:16 LST) по 18.03. 584 B.C. 02:03 UTC (04:20:16 LST)
с 17.03. 128 A.D. 00:10 UTC (02:22:25 LST) по 17.03. 128 A.D. 00:10 UTC (02:22:25 LST)
с 13.03. 568 A.D. 09:04 UTC (11:11:56 LST) по 13.03. 568 A.D. 09:04 UTC (11:11:56 LST)
с 29.03. 568 A.D. 04:32 UTC (06:45:23 LST) по 31.03. 568 A.D. 09:06 UTC (11:20:50 LST)

Решение ФиН : утро 20 марта 1185 г.н.э. см. http://www.chronologia.org/xpon3/17_04.html

Для ECLIPJ2000 у меня было близкое решение 31.03.1185 г.н.э., получается ФиН могли не учесть прецессию земной оси.

Зодиак Фальконетто (Мантуя)

     Разбивка зон зодиака в созвездиях

Расшифровка гороскопа Фальконетто (см. Г.В.Носовский, А.Т.Фоменко Египетские, русские и итальянские зодиаки. Исследования 2005-2008 годов. Глава 3. Итальянские средневековые зодиаки): Солнце (Гидра); Луна (женщина на облаке с раком в руке); Юпитер (Геркулес); Сатурн (человек в черном); Венера (красивая женщина на переднем плане) - в Раке Марс и Меркурий (двое мужчин на заднем плане, в отдалении) были не в Раке, а где-то в других созвездиях. Планета Lon,экл. созв. созв. Солнце 116 2.9 Рак/Близ. Луна 141.1 3.9 Лев/Рак Сатурн 132,9 3,6 Рак Юпитер 116,1 2,9 Рак/Близ. Марс 168.1 4.8 Лев Венера 124.6 3.2 Рак Меркурий 141.5 3.9 Лев/Рак Желательно - сильное сближение Юпитера с Солнцем.      Положения планет, рассчитанная ФиН - полночь по Гринвичу 3-4 июля старого стиля 1741 года (JD 2357142,5). Долгота планет на эклиптике J2000 в градусах, а также положения планет по шкале созвездий и названия созвездий Моя программа для SPICE в MatLab

Гороскоп Рождества Христа из книги Эбенезера Сибли (зодиак "AE") - расшифровка
(см. Г.В.Носовский, А.Т.Фоменко Египетские, русские и итальянские зодиаки. Исследования 2005-2008 годов. Глава 2. Русские и европейские зодиаки)

В древней книге астролога Сибли был напечатан гороскоп рождения Христа, на котором написана его дата 25 декабря 45 года (по юлианскому календарю), полночь, т.е. 23 декабря 45 г. по григореанскому календарю. В книге Фин рисунок гороскопа очень плохого качества, я привожу опубликованный в Википедии. Трактовка ФиН расположения планет в созвездиях долгота по эклиптике J2000 Планета в знаке границы созвезд. +- допуск 3/4 сосед. созвездия Lon Lon Lon Lon Юпитер Весы 215 239 182,75 259,25 Сатурн Близнец 92 118 65 132,25 Венера Козерог 296 326 273,5 343,25 Марс Рыбы 349 31 331,75 49,75 Меркурий Стрелец 266 296 245,75 318,5 Солнце Стрелец 266 296 245,75 318,5 Луна Телец 56 92 37,25 111,5 Решения в правой колонке.   моя трактовка расположения планет в созвездиях долгота по эклиптике J2000 Планета в знаке границы созвезд. +- допуск 3/4 сосед. созвездия Lon Lon Lon Lon Юпитер Весы 215 236     Сатурн Близнец 92 118     Венера Водолей 326 349     Марс Овен 31 56     Меркурий Козерог 296 326     Солнце Козерог 296 326     Луна Телец 56 92     Я задал программе свою расстановку планет, такую, какой увидел ее на гороскопе и получил единственное решение по Юлианскому календарю: с 25.12. 1 B.C. 01:21 UTC (03:39:50 LST) по   27.12. 1 B.C. 22:35 UTC (00:52:26 LST) расчет допусков к гороскопу Сибли AE программа расчета гороскопа Сибли AE

рис. взят с http://en.wikipedia.org/wiki/Ebenezer_Sibly Результат для строгих условий по ФиН, допуск +- 5°: с 12.12. 855 B.C. 13:53 UTC (16:27:00 LST) по 12.12. 855 B.C. 13:53 UTC (16:27:00 LST) с 25.12. 914 A.D. 07:32 UTC (09:41:51 LST) по 25.12. 914 A.D. 07:32 UTC (09:41:51 LST) с 29.12. 914 A.D. 16:42 UTC (18:48:59 LST) по 30.12. 914 A.D. 05:18 UTC (07:25:08 LST) Результат для нестрогих условий по ФиН, допуск +-5° + 3/4 соседнего созвездия: Результат: ... ... ... с 02.01. 237 B.C. 00:15 UTC (02:32:29 LST) по 07.01. 237 B.C. 16:37 UTC (18:51:47 LST) с 23.11. 1 B.C. 14:06 UTC (16:38:55 LST) по 23.11. 1 B.C. 14:06 UTC (16:38:55 LST) с 26.11. 1 B.C. 06:18 UTC (08:50:07 LST) по 01.12. 1 B.C. 14:11 UTC (16:41:30 LST) с 20.12. 1 B.C. 23:21 UTC (01:42:16 LST) по 20.12. 1 B.C. 23:21 UTC (01:42:16 LST) с 23.12. 1 B.C. 15:03 UTC (17:22:09 LST) по 27.12. 1 B.C. 03:08 UTC (05:25:25 LST) с 23.12. 60 A.D. 06:39 UTC (08:57:49 LST) по 25.12. 60 A.D. 15:29 UTC (17:46:52 LST) с 19.11. 617 A.D. 10:03 UTC (12:31:15 LST) по 21.11. 617 A.D. 08:09 UTC (10:36:53 LST) с 13.12. 617 A.D. 01:15 UTC (03:34:05 LST) по 18.12. 617 A.D. 15:23 UTC (17:38:32 LST) с 03.01. 678 A.D. 01:20 UTC (03:28:02 LST) по 03.01. 678 A.D. 01:20 UTC (03:28:02 LST) с 05.01. 678 A.D. 11:24 UTC (13:30:23 LST) по 12.01. 678 A.D. 06:06 UTC (08:10:06 LST) с 30.11. 914 A.D. 21:46 UTC (00:07:20 LST) по 07.12. 914 A.D. 10:43 UTC (13:01:20 LST) с 28.12. 914 A.D. 04:02 UTC (06:10:01 LST) по 03.01. 915 A.D. 16:15 UTC (18:20:34 LST) с 17.12.1151 A.D. 01:47 UTC (03:57:49 LST) по 17.12.1151 A.D. 01:47 UTC (03:57:49 LST) 19.12.1151 23:33 UTC (01:42:28 LST) по 25.12.1151 A.D. 11:03 UTC (13:09:49 LST) с 20.01.1295 A.D. 15:41 UTC (17:37:21 LST) по 20.01.1295 A.D. 15:41 UTC (17:37:21 LST) с 23.01.1295 A.D. 12:02 UTC (13:57:51 LST) по 24.01.1295 A.D. 09:43 UTC (11:39:06 LST) с 18.12.1531 A.D. 20:43 UTC (22:49:04 LST) по 24.12.1531 A.D. 10:50 UTC (12:54:05 LST) с 15.01.1532 A.D. 03:12 UTC (05:07:17 LST) по 20.01.1532 A.D. 23:28 UTC (01:22:35 LST) с 23.01.1592 A.D. 10:33 UTC (12:28:04 LST) по 27.01.1592 A.D. 08:46 UTC (10:40:36 LST)

Считаю, что гороскоп Сибли Рождества Христова был расчетным, а не зарисованным с неба.

Конвертации форматов.

  Меры измерений y = cspice_convrt( x, 'из', 'в');

Расстояние Углы Время 'METERS' 'KM' 'CM' 'MM' 'FEET' 'INCHES' 'YARDS' 'STATUTE_MILES' 'NAUTICAL_MILES' 'AU' 'PARSECS' 'LIGHTSECS' 'LIGHTYEARS', юл. 'RADIANS' 'DEGREES' 'ARCMINUTES' 'ARCSECONDS' 'HOURANGLE' 'MINUTEANGLE' 'SECONDANGLE' 'SECONDS' 'MINUTES' 'HOURS' 'DAYS' 'JULIAN_YEARS'('YEARS') 'TROPICAL_YEARS'

     Моя функция перевода углов из десятичной системы в десятично-шестидесятеричную.
% Целая часть десятичная, а остальное в шестидесятеричной системе.
function anstr = angl2minsec (an, prec, varargin)
if nargin == 1 ; sprec='0'; else; sprec=int2str(prec); end
if an<0 ; s_sign='-'; else s_sign=' '; end
an1=fix(an); an2=abs(rem(an,1)); an2=an2*0.6; an3=rem(an2*100,1)*100;
an2=fix(an2*100); an3=an3*0.6; anstr=sprintf('%s%s°%s''%s"', s_sign, int2str(an1), int2str(an2), num2str(an3, strcat('%2.',sprec,'f')));
return

 Перевод угла 0-360 в часы 0-24
     rah=cspice_timout((ra/15)* (cspice_spd/24) + (cspice_spd/2), ...
         'HR:MN:SC.## ::RND::TDT' );
Констранта cspice_spd содержит в себе ровно день в секундах TDT (не UT). В примере угол переводится в часы и + 12 часов.

Преобразование массивов

    st_win=(reshape(st_win,2,[])'); % преобразование окон времени в 2 колонки 1) начало события 2) конец события
       st_win= st_win (:,1); % выбор одной колонки, там где колонки равны, например, равноденствие

Некоторые функции Ерана Офека

   Оригинал: http://wise-obs.tau.ac.il/~eran/matlab.html (или выборочно здесь)

Константы

      Земля

	Наименование	Данные	или	или
	Перигелий	147 098 074 км	0,9832898912 а. е.	~ 3 января
	Афелий	152 097 701 км	1,0167103335 а. е.	~ 4 июля
	Сидерический период обращения	365,256366 дней	365 дн. 6 ч. 9 мин. 10 сек.
	Звёздные сутки	23:56:4.091
	Полярное сжатие	0,0033528
	Полярный радиус	6 356,8 км
	Средний радиус	6 371,0 км
	Окружность большого круга	40009,88 км

    Интересные снимки

Взрыв галактического ядра

Снимок вселенной телескопом космического аппарата "Планк". В центре, горизонтально, расположена наша галактика "Млечный путь" Топ-100 Хаббл и другие телескопы

http://www.biguniverse.ru/posts/news/plank-delaet-poputnyie-otkryitiya

Источники (литература) >>>