Расчет посадки космического аппарата ‘Луна-Глоб’ на Луну

Расчет посадки космического аппарата ‘Луна-Глоб’ на Луну

Введение

Данная научная работа проводилась в рамках проекта «Луна-Глоб», реализуемого в НПО имени С.А.Лавочкина. Целью этого проекта является исследование лунной поверхности с помощью посадочного аппарата проводящего мониторинг автоматического зонда, движущегося по замкнутой периодической орбите с миссией дистанционного исследования Луны. Запуск космического аппарата «Луна-Глоб» запланирован на 2015 год на ракетном носителе «Зенит».

Так как для решения поставленной задачи оценки пригодности выбранного места посадки на поверхности Луны не было найдено готовое программное обеспечение (ПО), было принято решение реализовать свое ПО, удовлетворяющее необходимым критериям, которое для любой заданной точки наблюдения:

·Вычисляет угол места над горизонтом для космического аппарата, движущегося по заданной траектории в течение определенного промежутка времени;

·Вычисляет максимальный угол возвышения высот неровностей лунного рельефа в окрестности выбранного пункта наблюдения;

·Оценивает возможность виденья движущегося космического аппарат с заданной точки наблюдения на каждый момент времени;

·Вычисляет максимальный период времени, в течение которого виден рассматриваемый объект на орбите.

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

Глава 1.Выбор оптимального места для посадки аппарата

При выборе места для посадки аппарата на Луну необходимо выбрать оптимальное положение на ее поверхности. Критерии выбора данного положения зависят от цели миссии, рельефа лунной поверхности, видимости Солнца и других факторов. Их совокупность определяет требования, по которым отбираются районы, пригодные для посадки аппарата.

Одним из таких факторов, играющих важную роль в выборе оптимального места для посадки космического аппарата (КА), является состав почвы на выбранном участке поверхности. Лунные запасы воды в форме льда, которые могут быть использованы для получения ракетного топлива, составляют 1,6 млрд. тонн и сосредоточены в основном на полюсах. Также на Луне в изобилии присутствуют редкоземельные элементы, залегающие под поверхностью. Таким образом, Луна является территорией, богатой ценными ресурсами.

Наиболее интересной для исследования является поверхность в северной части Луны, так как в этой области почва богата водородом.

Актуальной целью посадки за последнее время является наблюдение за космическим аппаратом, двигающимся по гало-орбите в точке либрации .

.1 Ограниченная задача трех тел

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

Два обращающихся тела называются основными телами или точками конечной массы, а третье тело — точкой «нулевой массы» или «пассивно-гравитирующей точкой». Масса третьего тела намного меньше масс и . Это следует из принятого предположения, что масса не влияет на движение масс и . Круговое движение основных тел и показано на рис. 1.

Рис. 1. Взаимно притягивающееся движение масс , и в ограниченной задаче трех тел.

В 1772 г. Лагранж доказал, что существует определенное количество частных случаев в задаче о трех телах, в которых может быть найдено точное решение. Если заданы массы двух тел и их положение на плоскости (рис. 2), то рассматриваемые частные случаи движения в этой плоскости получаются при расположении третьего тела в одной из пяти точек, называемых точками либрации Лагранжа. Первые три точки либрации располагаются в определенных точках прямой, соединяющей обе заданные массы, причем одна между ними, а две другие — вне их. Их называют коллинеарными точками либрации. Четвертая и пятая точки являются вершинами двух равносторонних треугольников, в которых остальные вершины заняты заданными массами. Они называются треугольными точками либрации.

Рис. 2. Положение двух тел на плоскости и пять точек либрации Лагранжа.

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

Таким образом:

)если три тела расположены на одной прямой, то они обращаются, оставаясь на ней, вокруг общего центра масс;

)если три тела расположены в вершинах равностороннего треугольника, то они обращаются вокруг общего центра масс так, что треугольник все время остается равносторонним.

.2 Точки либрации

В системе Земля-Луна первые три точки либрации находятся на вращающейся линии, соединяющей Землю и Луну: точка лежит между планетами, вторая точка находится за Луной, а третья коллинеарная точка расположена с обратной стороны Земли по отношению к Луне. Остальные две точки либрации и находятся с двух сторон вне вращающейся линии.

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

Рис. 3.Пять точек либрации в системе Земля-Луна.

Точки , и являются неустойчивыми. Так как, если объект, помещенный в коллинеарную точку Лагранжа, слегка смещается вдоль прямой, соединяющей Землю и Луну, то сила, притягивающая объект к тому телу, к которому он приближается, увеличивается, а сила притяжения со стороны другого тела, наоборот, уменьшается. В результате объект будет все больше удалятся от положения равновесия.

Однако существуют стабильные замкнутые квазипериодические и периодические орбиты, такие как Лиссажу и гало-орбиты, которые колеблются около этих точек. То есть космический аппарата, совершающий движение по гало-орбите, будет оставаться на ней в течение длительного времени (рис.4).

Рис. 4. Гало-орбиты в системе Земля-Луна.

Объект, такой как космический аппарат, который смещен от точки либрации, будет колебаться вокруг точки с периодом определенным тем, насколько далеко он смещается в Y и Z направлениях (рис.5). Параметр τ является углом, определяющим положение космического аппарата на заданной гало-орбите и аналогичен истинной аномалии при полете по эллиптической орбите. Он измеряется в положительном направлении от оси +Z около оси +X от 0 ̊ до 360 ̊.

Рис 5. Движение космического аппарата по гало-орбите.

Среди пяти точек либрации системы Земля-Луна более актуальными для исследования человеком космического пространства являются две, находящиеся ближе всего к Луне — и . Они расположены около ближней и дальней сторон Луны соответственно, если смотреть с Земли. Однако лучше исследовать обратную сторону Луны, которая является одним из приоритетных мест для исследования космогонии и истории Солнечной системы. Луна защищает поверхность на её обратной стороне от наземных радиошумов, что облегчает изучение низкочастотных сигналов (ниже 100 МГц).

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

Таким образом, в данной работе рассматривается космический аппарат на лунной поверхности, находящейся ближе к северу (диапазон по широте равен от 60̊ до 90̊ ), с которого ведется наблюдение за объектом, совершающим движение по гало-орбите в точке либрации (рис. 6).

Рис. 6. Космический аппарат в точке либрации .

Глава 2. Определение видимости КА без учета рельефа лунной поверхности

Видимость космического аппарата, двигающегося по заданной траектории, в данный момент времени будет зависеть от нескольких основных факторов:

·его положения на гало-орбите;

·координат точки наблюдения с поверхности Луны, то есть выбора места посадки аппарата, с которого ведется наблюдение;

·высоты рельефа лунной поверхности.

Однако для начала рассмотрим упрощенную задачу, где неровности горизонта на поверхности Луны не будут преградой для отслеживания КА. То есть будем считать, что заданная точка наблюдения находится на гладкой сфере с постоянным радиусом равным радиусу Луны.

Тогда видимость КА в данный момент времени будет зависеть только от его положения на гало-орбите и координат выбранной точки наблюдения с поверхности Луны. Информация о положении КА на гало-орбите в данный момент времени может быть получена из файла «ОРБИТА.txt», в котором хранятся данные по координатам и скорости аппарата в системе координат J200 (она будет рассмотрена ниже) и дате и времени в формате UTC (рис. 7).

Рис. 7. Данные о КА, двигающемся по гало-орбите, в файле «ОРБИТА.txt».

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

.1 Преобразование данных в единую систему координат

Перед выполнением задачи определения видимости КА с заданной точки на Луне необходимо представить все данные о положении в одной системе координат — СК J200. Для этого сначала переведем координаты точки на поверхности Луны в декартовы, а далее в систему координат J200.

.1.1 Сферическая система координат

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

Рис. 8. Сферические координаты и .

Зенит — это направление вертикального подъема над произвольно выбранной точкой (точкой наблюдения), принадлежащей так называемой фундаментальной плоскости. Азимут — угол между произвольно выбранным лучом фундаментальной плоскости с началом в точке наблюдения и другим лучом этой плоскости, имеющим общее начало с первым. В сферической системе координат, фундаментальная плоскость — это плоскость XY. Зенит — некая удаленная точка, лежащая на оси Z и видимая из начала координат. Азимут отсчитывается от оси X до проекции радиус-вектора на плоскость XY (рис.9).

Рис.9. Сферическая система координат.

Три координаты определены следующим образом:

· — расстояние от начала координат до заданной точки .

· ̊ — угол между осью и отрезком, соединяющим начало координат и точку .

· ̊ — угол между осью и проекцией отрезка, соединяющего начало координат с точкой , на плоскость .

.1.2 Связь сферических координат с декартовыми

Если заданы сферические координаты точки, то переход к декартовым осуществляется по нижеприведенным формулам:

(2.1)

где — декартовы координаты точки;

— сферические координаты.

Чтобы перейти обратно от декартовых координат к сферическим, нужно применить формулы:

(2.2)

Таким образом, пользуясь формулами (2.1) переводим сферические координаты заданной точки на поверхности Луны в декартовые .

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

2.1.3 Селенографическая вращающаяся система координат (СВСК)

Положение и скорость объекта относительно Луны удобно определяется в подвижной системе координат (СК), оси которой вращаются вместе с Луной. Такая СК называется селеноэкваториальной луноцентрической системой координат. Основной координатной плоскостью является плоскость истинной экватора Луны. За основную точку отсчета принята точка пересечения первого радиуса с лунным экватором. Первый радиус определяется пересечением плоскости лунного меридиана, проведенной через центр масс Земли, с плоскостью лунного экватора в момент времени, когда средняя долгота Луны равна средней долготе ее восходящего узла, и направлен в сторону Земли (рис.10).

Рис.10. Селеноэкваториальная луноцентрическая система координат.

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

Селенографическая долгота отсчитывается по лунному экватору от основной точки к востоку. Селенографическая широта — это острый угол между луноцентрическим радиусом-вектором и плоскостью лунного экватора, отсчитывается от экватора Луны по лунным меридианам.

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

Рис. 11. Селенографические долгота и широта .

2.1.4 Инерциальная система координат J2000 (СК J2000)

В середине 80-х годов Международный астрономический союз (IAU) рекомендовал использовать для отсчета координат и угловых параметров небесных тел инерциальную систему координат J2000, базирующуюся на фундаментальном каталоге FK5. Ранее основой астрономической системой координат была В1950.0 (каталог звезд FK4). СК J2000 относится к эпохе J2000, соответствующей юлианской дате 2451545.0 или 11 час. 58 мин. 55.816 сек. 1 января 2000 года времени UTC. Плоскость XY СК J2000 совпадает с плоскостью среднего экватора J2000 эпохи. Ось Х направлена в точку среднего весеннего равноденствия (пересечения экватора и эклиптики). Ось Z направлена ортогонально плоскости XY в сторону северного полюса. Ось Y дополняет систему координат до правой тройки (рис.12).

Рис. 12. Направления осей в системе координат J2000.

Переход в систему координат J2000

Чтобы перевести координаты из селенографический вращающейся системы координат в СК J2000, нужно умножить их на матрицу перехода , которая вычисляется по следующему алгоритму:

Пусть − интервал времени в сутках, начиная с 1 января 12 часов 2000 года (эпоха J2000), который вычисляется по формуле:

, (2.3)

где — юлианская дата в момент времени t;

— юлианская дата, соответствующая эпохе J2000.

Тогда стандартные координаты и среднего лунного экватора рассчитываются по соотношениям:

(2.4)

(2.5)

А угол поворота нулевого меридиана по соотношению:

(2.6)

В формулах (2.5) и (2.6) − угловая скорость вращения Луны (º/сутки).

А величина T определяется по формуле:

, (2.7)

Величины находятся по нижеприведенным формулам:


Тогда матрица перехода из инерциальной СК J2000 к СВСК вычисляется следующим способом:

(2.8)

Матрица перехода от СВСК к инерциальной СК J2000 вычисляется по соотношению:

, (2.9)

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

В формуле (2.3) можно увидеть, что для нахождения интервала времени необходимо знать юлианскую дату для данного момента t. Таким образом, прежде чем переводить координаты в СК J2000, надо перевести время из формата UTC в формат юлианской даты.

.1.6 Юлианская дата

Юлиа́нская дата (JD) — астрономический способ измерения времени, при котором считается число дней, прошедших начиная с полудня 1 января 4713 до н. э. юлианского календаря или, что то же самое, 24 ноября 4714 г. до н. э. григорианского календаря. Первый день имел номер 0. С тех пор по настоящее время прошло немногим менее 2,5 миллионов дней. Даты сменяются в полдень UT или TT. Для точного обозначения времени применяют дробную часть, например, JD = 2451545,25 соответствует 18 часам 1 января 2000 г., 3 часа дня 2 августа 1942 г. — JD 2430574,125, 13,5 июня 1944 г. — JD 2431255,0.

Вычисление юлианской даты по дате календаря:

Для дальнейших вычислений введем новые обозначения:

· — год, если использовать его для дат до нашей эры, необходимо перевести год до н. э. в отрицательный год (например, 10 г. до н. э. = −9);

· — номер месяца;

· — день месяца;

· — часы от 0 до 23;

· — минуты от 0 до 59;

· — секунды от 0 до 59, могут содержать дробную часть;

· — это номер юлианского дня (англ. Julian Day Number), который начинается в полдень числа, для которого производятся вычисления;

· — юлианская дата, содержащая дробную часть.

Сначала вычислим промежуточные коэффициенты:

После этого можно вычислить номер юлианского дня () по формуле:

(2.10)

Все деления целочисленные, то есть остатки деления отбрасываются.

Формула (2.10) справедлива для дат после 23 ноября — 4713 г. (то есть 23 ноября 4714 г. до н.э.)

Для перехода к юлианской дате, содержащей дробную часть, можно воспользоваться формулой:

(2.11)

При делении в формуле (2.11) дробная часть не отбрасывается. Сутки не должны содержать високосной секунды.

Таким образом, вычислив значение юлианской даты для заданного момента времени t, переводим координаты точки на поверхности Луны в инерциальную систему координат J2000 и получаем точку наблюдения . Убедившись, что все данные находятся в одной системе координат, можно приступать к следующему самому важному этапу определения видимости КА.

2.2 Расчет угла места КА над горизонтом Луны

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

.2.1 Угол места над горизонтом

Угол места (элевация) — угловая высота наблюдаемого объекта (земного предмета, летательного аппарата, небесного светила и др.) над истинным горизонтом. То есть это угол возвышения или угол подъема космического аппарата над линией горизонта (рис. 13). Очевидно, что если угол места меньше нуля, то аппарата над горизонтом не видно.

Рис. 13. Угол места КА над горизонтом.

Угол места для космического аппарата, совершающего движение по гало-орбите, рассчитывается по нижеприведенному алгоритму.

2.2.2 Последовательность действий для определения угла места

На рис. 14 изображен угол возвышения КА над лунным горизонтом . — точка наблюдения на поверхности Луны, — положение КА на гало-орбите.

Рис. 14. Угол места КА над горизонтом Луны.

Находим радиус планеты , в данной задаче он равен радиусу Луны , и расстояние от центра Луны до космического аппарата по координатам положения КА на орбите.

, .

Учитывая, что вектора и выходят из центра Луны, следовательно, имеют координаты:

,

Получаем:

(2.12)

(2.13)

.Находим угол между векторами и :

,

где — скалярное произведение векторов и , рассчитываемое по формуле:

.

и — длины векторов и , рассчитываемые по формулам (2.12) и (2.13).

.Находим вектор как разность векторов и :

.

Получаем координаты .

Вектор равен расстоянию от точки наблюдения на поверхности Луны до положения КА на гало-орбите.

Находим длину вектора :

.Аналогично пункту 2 находим угол между векторами и :

.Из треугольника находим угол места :

2.2.3 Частный случай положения КА относительно точки наблюдения

Отдельно стоит рассмотреть случай, когда КА находится прямо над точкой наблюдения . Тогда вектора и — коллинеарные и угол места . Однако, если рассматриваемый объект находится прямо над наблюдателем, то очевидно, что КА будет видно (рис. 15). Поэтому в случае если вектора и — коллинеарные, то есть выполняется условие:

, (2.14)

то КА видно.

В формуле (2.14) — векторное произведение векторов и , которое рассчитывается по формуле:

Рис. 15. Случай расположения КА непосредственно над наблюдателем.

Еще одним важным этапом в определении видимости космического аппарата (КА) является нахождение границ области на поверхности Луны, из которой будет виден КА («области интересов»). Это позволит сократить количество рассматриваемых точек наблюдения на поверхности Луны и сократит объем вычислительной работы. Таким образом, если фиксированная точка не входит в область интересов, то космический аппарат, находящийся на орбите, не будет виден из нее.

.3 Определение границ области на поверхности Луны, из которой виден КА

Область, состоящая из множества точек на поверхности Луны, из которых виден наблюдаемый объект, двигающийся по заданной траектории, называется зоной видимости объекта.

В рассматриваемой задаче подвижным объектом является КА, двигающийся по заданной гало-орбите. Зона видимости КА определяется несколькими факторами:

·Касательными, проведенными из точки положения КА на данный момент времени к лунной сфере;

·Проекцией точки положения КА на поверхность Луны;

Рис. 16. Зона видимости КА с поверхности Луны.

На рис. 16 красным цветом заштрихована область лунной поверхности, откуда виден КА на данный момент времени. O — центр Луны, точки B и C — точки касания лунной сферы, точка — проекция положения КА на сферу.

Нетрудно заметить, что зоной видимости КА является параболический гиперболоид, образующийся при пересечении Лунной сферы с центром в точке O и радиуса и конуса с вершиной в точке КА, чьими образующими являются касательные к сфере.

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

Как известно, через любые три точки в пространстве можно провести плоскость. Тогда проведем плоскость через точку наблюдения на поверхности , центр Луны O и проекцию положения КА на лунную поверхность K ̍. Эта плоскость пересекает Лунную сферу, и сечением является окружность с центром в точке O и радиусом (рис. 17).

Рис. 17. Зона видимости КА в двумерной задаче.

— угол между радиус-вектором к точке касания С и вектором, соединяющим центр Луны с точкой положения КА.

— угол между радиус-вектором к точке наблюдения P и вектором, соединяющим центр Луны с точкой положения КА.

То есть — это угловое отклонение точки наблюдения от проекции положения КА на орбите . По рис. 17 нетрудно заметить, что, если , то точка наблюдения принадлежит зоне видимости КА, если , то — не принадлежит.

Нижеприведенный алгоритм позволяет определить, принадлежит ли заданная точка на поверхности Луны зоне видимости КА или нет.

.Находим координаты точки проекции :

Так как точка лежит на том же векторе, что и точка положения КА на орбите, то первые две координаты у них одинаковые, а третья координата для точки равна . То есть, учитывая, что точка КА имеет координаты , то координаты проекции равны:

.Из прямоугольного треугольника видно, что:

(2.15)

.Далее по формуле для угла между векторами находим угол :

(2.16)

Так как арккосинус убывающая функция, то для того, чтобы , нужно чтобы выполнялось условие: .

Учитывая формулы (2.15) и (2.16), получаем:

(2.17)

Так как и , то формула (2.17) принимает вид:

(2.18)

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

.4 Реализация разработанного алгоритма на языке С++

По разработанному алгоритму определения видимости КА с заданной точки на поверхности Луны была написана программа на языке С++, определяющая степень обзора на движущийся по гало-орбите объект (Приложение1).

Программа считывает из файла «ОРБИТА.txt» данные о положении КА на гало-орбите для каждого момента времени t и для заданной пользователем точки на поверхности Луны рассчитывает угол места над горизонтом.

Ниже приведены основные фрагменты кода программы.

Перевод координат точки из сферических в декартовы:

//декартовы (x,y,z): (dec[1],dec[0],dec[2])

//сферические (r, teta,fi): (sph[2],sph[1],sph[0])

void S_to_SG(double sph[3], double dec[3])

{//широта долгота радиус[0] = sph[2] * sin(DTR*sph[0])*sin(DTR*sph[1]);[1] = sph[2] * cos(DTR*sph[0])*sin(DTR*sph[1]);[2] = sph[2] * cos(DTR*sph[1]);}

Расчет матрицы и переход из СВСК в СК L2000:

void martix_J_to_SG(double Rj[3], double M[3][3], int time[6]){dt, T, alph, bet, w, e[13];

int p = 10;gaml = 13.17635815; // угл. скорость вращения Луны в градусах.

dt = DateTimeToJulDate(time[0], time[1], time[2], time[3], time[4], time[5]) — 2415020.0;[0] = (125.045 — 0.0529921*dt)*DTR;[1] = (250.089 — 0.1059842*dt)*DTR;[2] = (260.008 + 13.0120009*dt)*DTR;[3] = (176.625 + 13.3407154*dt)*DTR;[4] = (357.529 + 0.9856003*dt)*DTR;[5] = (311.589 + 26.4057084*dt)*DTR;[6] = (134.963 + 13.0649930*dt)*DTR;[7] = (276.617 + 0.3287146*dt)*DTR;[8] = (34.226 + 1.7484877*dt)*DTR;[9] = (15.134 — 0.1589763*dt)*DTR;[10] = (119.743 + 0.0036096*dt)*DTR;[11] = (239.961 + 0.1643573*dt)*DTR;[12] = (25.053 + 12.9590088*dt)*DTR;= dt / 36525;= (269.9949 + 0.0031*T — 3.8787*sin(e[0]) — 0.1204*sin(e[1]) + 0.07*sin(e[2]) — 0.0172*sin(e[3]) + 0.0072*sin(e[5]) — 0.0052*sin(e[9]) — 0.0043*sin(e[12]))*DTR;= (66.5392 + 0.013*T + 1.5419*cos(e[0]) + 0.0239*cos(e[1]) — 0.0278*cos(e[2]) + 0.0068*cos(e[3]) — 0.0029*cos(e[5]) + 0.0009*cos(e[6]) + 0.0008*cos(e[9]) — 0.0009*cos(e[12]))*DTR;= (38.3213 + gaml*dt + gaml + dt — 1.4*(p, (-12))*pow(dt, 2) + 3.561*sin(e[0]) + 0.1208*sin(e[1]) — 0.0642*sin(e[2]) + 0.0158*sin(e[3]) + 0.0252*sin(e[4]) — 0.0066*sin(e[5]) — 0.0047*sin(e[6]) — 0.0046*sin(e[7]) + 0.0028*sin(e[8]) + 0.0052*sin(e[9]) + 0.004*sin(e[10]) + 0.0019*sin(e[11]) — 0.0044*sin(e[12]))*DTR;[0][0] = -sin(alph)*cos(w) — cos(alph)*sin(bet)*sin(w);[0][1] = cos(alph)*cos(w) — sin(alph)*sin(bet)*sin(w);[0][2] = sin(w)*cos(bet);[1][0] = sin(alph)*sin(w) — cos(alph)*sin(bet)*cos(w);[1][1] = -cos(alph)*sin(w) — sin(alph)*sin(bet)*cos(w);[1][2] = cos(w)*cos(bet);[2][0] = cos(alph)*cos(bet);[2][1] = sin(alph)*cos(bet);[2][2] = sin(bet);SG_to_J(double Rsg[3], double Rj[3], int time[6])

{double M[3][3];_J_to_SG(Rj, M, time);_matrix(M);_on_matrix(Rsg, M, Rj);

}

Перевод времени из формата UTC в юлианскую дату:

double DateTimeToJulDate(int Year, int Month, int Day, int Hour, int Min, double Sec){ long MjdMidnight;b;(Month <= 2){ Month += 12; —Year; }((10000L * Year + 100L * Month + Day) <= 15821004L)= -2 + ((Year + 4716) / 4) — 1179;= (Year / 400) — (Year / 100) + (Year / 4);= 365L * Year — 679004 + b + int(30.6001*(Month + 1)) + Day;MjdMidnight + ((Hour + Min / 60.0 + Sec / 3600.0) / 24.0) + 2400000.5;

Расчет угла места над горизонтом:

// OP — вектор в J2000 из центра Луны в точку на поверхности Луны;

// OK — вектор в J2000 из центра Луны в координаты КА на гало-орбите.

double elevation_angle(double * OP, double * OK)

{int k;lenOP, lenOK, lenPK, betta, gamma, alpha;

double PK[3];// вектор от точки P на Луне в положение КА= length(OP);//длина вектора OP= length(OK);//длина вектора OK

//Вектор PK(k = 0; k < 3; k++) PK[k] = OK[k] — OP[k];

lenPK = length(PK);//длина вектора PK= angle(OP, OK);//угол между векторами OP и OK= angle(PK, OK);//угол между векторами PK и OK

//Перевод из радиан в градусы.

alpha = 90.e0 — gamma — betta;alpha;}

В результате выполнения программы, получаем файл «угол_места.txt» с данными (рис. 18):

·координаты точки наблюдения Point в сферических координатах: point x, point y, point z;

·координаты КА на гало-орбите в СК J2000: KA J2000 x, KA J2000 y, KA J2000 z;

·Время и дата, в момент которых наблюдается КА в формате UTC: year.month.day hours:minutes:seconds;

·Угол места над горизонтом в градусах;

·Виден КА в данный момент времени с точки Point или нет.

Рис. 18. Данные в файле «угол_места.txt» после выполнения программы.

На рис. 18 представлены результаты программы для точки наблюдения на поверхности Луны с координатами по широте и долготе равными 74 ̊ и 65 ̊ .

Если построить график зависимости угла места для КА от времени (рис. 19), то можно увидеть, что угол места КА над горизонтом периодически меняет свой знак. Это объясняется тем, что наблюдаемый КА двигается по периодической гало-орбите, поэтому он время от времени то появляется в зоне видимости, то исчезает из нее.

Рис. 19. График зависимости угла места от времени для точки с координатами по широте 74 ̊ и долготе 65 ̊ .

Возьмем другую точку на поверхности Луны. На рис. 20 представлен график зависимости угла места КА от времени для точки с координатами по широте 74 ̊ и долготе 175 ̊ . Если сравнить рис. 19 и рис. 20, то можно заметить, что видимость в точке (74, 175) хуже, так как угол места КА в ней меньше. Это означает, что видимость КА с заданной точки на поверхности Луны зависит от координат точки наблюдения за объектом.

Рис. 20. График зависимости угла места от времени для точки с координатами по широте 74 ̊ и долготе 175 ̊ .

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

Глава 3. Определение видимости КА с учетом рельефа лунной поверхности

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

Информация о неровностях рельефа считываются из файла «moon.txt», где находятся данные высот, как отклонения от среднего радиуса Луны. Размеры файла составляют , так как в нем представлены значения высот для точек, полученных путем разбиения лунной поверхности на 15520 частей по широте (от 60 ̊ до 90 ̊ ) и на 15520 частей по долготе (от 0 до 360 ̊ ).

Для того чтобы исследовать видимость КА с точки на поверхности Луны с учетом рельефа, нужно понять, есть ли хоть одна высота, закрывающая обзор на КА в данный момент времени. Параметром, характеризующим степень возвышения высоты над поверхностью, является угол места. То есть необходимо вычислить угол места высоты над горизонтом и сравнить его с углом места для КА. Если угол места над горизонтом для высоты больше или равен угла места для КА, то эта высота закрывает обзор на объект (рис.21).

Рис. 21. Угол возвышения высоты над горизонтом.

На рис. 21 — угол места космического аппарата, находящегося в положении КА на гало-орбите, — угол места для неровности лунной поверхности высотой , точка O — центр Луны.

Таким образом, если выполняется условие:

, (3.1)

то высота не закрывает обзор из выбранной точки на КА.

Если условие (3.1) не выполняется, то есть , то КА невидно.

.1 Расчет угла места для высоты

При расчете угла места для высоты необходимо учитывать тот факт, что точка наблюдения P может находиться на высоте над уровнем горизонта Луны (рис. 22). Поэтому значения высот относительно точки наблюдения будут меньше на величину . То есть в ниже приведенном алгоритме расчета угла места, высоты неровностей будут равны: .

Рис. 22. Положение точки наблюдения P на высоте над горизонтом.

3.1.1 Алгоритм нахождения угла места для высоты

1.Находим радиус-вектор от центра Луны до точки наблюдения на поверхности и радиус-вектор до точки , где возвышается неровность рельефа высотой :

,

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

.Находим угол между векторами и :

(3.2)

Так как точка наблюдения находится на высоте над горизонтом, то , где .

Таким образом, от задачи с точкой наблюдения на высоте над горизонтом сферы радиуса переходим к задаче с точкой наблюдения на поверхности сферы радиуса (рис. 22).

.В треугольнике по теореме косинусов находим :

(3.3)

Учитывая, что , формула (3.3) принимает вид:


.В треугольнике по теореме синусов находим угол :


Учитывая, что , получаем:

(3.4)

.Вычислив углы и из формул (3.2) и (3.4), находим угол :

3.1.2 Частный случай значения высоты неровности рельефа

В случае, когда высота неровности отрицательна или равна нулю (кратер, море), а высота , то есть точка находится на возвышенности, нет смысла рассчитывать для нее угол места над горизонтом. Так как при неровность не будет преградой для обзора КА из точки .

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

.2 Нахождение области высот на Луне, которые являются возможными преградами для обзора КА

3.2.1 Определение области высот на поверхности Луны

Для нахождения интересующей нас области на лунной поверхности, содержащей неровности Луны, высоты которых могут стать преградой для наблюдения КА из данной точки , делаем проекцию точки положения КА на лунную сферу — .

Из рис. 23 можно увидеть что, так называемая «область высот» представляет собой дугу на лунной сфере, начало которой лежит в точке наблюдения , а конец — в точке проекции .

Рис. 23. Дуга, представляющая собой область высот на лунной поверхности.

На рис. 23 красным цветом обозначена дуга , представляющая собой «область высот».

Таким образом, для определения видимости КА из точки на поверхности Луны с учетом рельефа необходимо вычислить угол места для каждой высоты , находящейся на дуге . Для этого разбиваем дугу на равных частей (маленьких дуг) с заданным шагом разбиения . На каждой дуге , где , , находим угол места для i-ой высоты (рис. 24).

Рис. 24. Разбиение дуги на интервалы.

.2.2 Разбиение дуги, представляющей область высот, на интервалы

Ниже представлен алгоритм разбиения дуги на интервалы.

.Для начала нужно найти координаты проекции точки положения на лунную сферу. Так как вектора и коллинеарные, то координаты точек и отличаются только третьей координатой. Учитывая, что точка находится на высоте над горизонтом Луны, то координаты проекции точки равны: , где

.Зададим единичный вектор :

, (3.5)

где — векторное произведение.

Для векторов и вычисляется по формуле:

; ; .

В формуле (3.5) модуль векторного произведения равен:

.

Вектор задает плоскость, в которой лежат точка наблюдения и точка положения на гало-орбите, и перпендикулярен ей (рис. 25). Он будет использоваться в дальнейшем для равномерного перемещения по дуге с шагом s.

Рис. 25. Единичный вектор .

3.Вычисляем угол между векторами и (рис. 26):

,

где — скалярное произведение;

а длины векторов равны: .

Рис. 26. Угол между векторами и .

4.Далее задаем шаг разбиения . От выбора его величины зависит точность решения задачи определения видимости КА: чем шаг меньше, тем точность больше.

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

(3.6)

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

.Далее производить разбиения дуги будем, сдвигая вектор на угол , зависящий от шага разбиения :

Для этого вычислим матрицу поворота :

(3.7)

В формуле (3.7) написана матрица поворота относительно оси вращения, заданной единичным вектором , на угол .

Таким образом, умножая вектор на матрицу , получаем новый вектор , конец которого лежит на дуге (рис. 27).

Рис. 27. Вектора на дуге , полученные поворотом на угол .

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

3.3. Реализация разработанного алгоритма на языке С++.

По разработанному алгоритму определения видимости КА с заданной точки на поверхности Луны с учетом ее рельефа была написана программа на языке С++, определяющая степень обзора на движущийся по гало-орбите объект (Приложение 2).

Программа считывает из файла «ОРБИТА.txt» данные о положении КА на гало-орбите для каждого момента времени t и из файла «moon.txt» информацию о высотах лунной поверхности. Затем для заданной пользователем точки на поверхности Луны рассчитывает угол места КА над горизонтом и максимальный угол возвышения высот , находящихся на дуге . И в зависимости от результата сравнения углов и получает характеристику видимости космического аппарата из заданной точки.

Ниже приведены основные фрагменты кода программы.

Считывание из файла «moon.txt» значения высоты для заданной точки P:

//считывание значения высоты для точки (lat, lon)

double moon_height_get(moon_height*Height, double lat, double lon){

int iLat, iLon;

// переводим lat в целочисленное значение

iLat = (int)(Height->lat_N * (lat — Height->lat_min) / (Height->lat_max — Height->lat_min));(iLat < 0) iLat = 0;(iLat >= Height->lat_N) iLat = Height->lat_N — 1;

// переводим lon в целочисленное значение= (int)(Height->lon_N * (lon — Height->lon_min) / (Height->lon_max — Height->lon_min));(iLon < 0) iLon = 0;(iLon >= Height->lon_N) iLon = Height->lon_N — 1;Height->height[iLat][iLon];}

moon_height_read(FILE_MOON, 60, 90, 0, 360, -15520, 15520, &MoonHeight);

Так как в файле «moon.txt» высоты заданы дискретно, то для произвольной точки (lat, lon) считаем, что значение высоты равно ее значению в ближайшей точке из файла.

Разбиение дуги, представляющей область высот, на интервалы:

AngleNStep = (int)(angle(vecOP, vecOK) / AngleStep + 1);(i = 1; i <= AngleNStep; i++){_rotation(vecN, i*AngleStep, MatrixRotation);_on_matrix(vecOP, MatrixRotation, vecOS);(k = 0; k < 3; k++) pointS_J[k] = MoonCenterJ[k] + vecOS[k];_to_S(pointS_J, pointS_S,KA[iO].time);_S[2] = MoonRadius + moon_height_get(MoonHeight, pointS_S[0], pointS_S[1]);_to_J(pointS_S, pointS_J, KA[iO].time);(k = 0; k < 3; k++) vecOS[k] = pointS_J[k] — MoonCenterJ[k];= elevation_angle(vecOP, vecOS);(i == 1){_Moon[iO] = a;

}else{(angle_Moon[iO] < a) angle_Moon[iO] = a;}

//Получить матрицу поворота вокруг вектора v на угол angle

void matrix_rotation(double v[3], double angle, double M[3][3]){cosa, sina;x, y, z;= cos(angle*DTR);= sin(angle*DTR);= v[0];= v[1];= v[2];[0][0] = cosa + (1 — cosa)*x*x;[0][1] = (1 — cosa)*x*y — sina*z;[0][2] = (1 — cosa)*x*z + sina*y;[1][0] = (1-cosa)*y*x + sina*z;[1][1] = cosa + (1 -cosa)*y*y;[1][2] = (1 — cosa)*y*z — sina*x;[2][0] = (1-cosa)*z*x — sina*y;[2][1] = (1 — cosa)*z*y + sina*x;[2][2] = cosa+(1 — cosa)*z*z;}

В данной программе при разбиении дуги на интервалы значение шага разбиения принималось равным 0,5 ̊. Для увеличения точности вычислений следует уменьшать его значение.

Вычисление максимального периода времени (в секундах), в течении которого виден КА:

double TimeMax,Time,TimeStart,TimeFinish;= 0.e0;;= -1.e0;(i = 0; i < NKA; i++){= DateTimeToJulDate(KA[i].time[0], KA[i].time[1], KA[i].time[2], KA[i].time[3], KA[i].time[4], KA[i].time[5]);(angle_KA[i] > angle_Moon[i]){(TimeStart < 0){= Time;= Time;

}else{= Time;}

}else{(TimeStart>0){(TimeFinish — TimeStart > TimeMax)= TimeFinish — TimeStart;

}else{//do nothing

}}}(TimeStart > 0){(TimeFinish — TimeStart > TimeMax)= TimeFinish — TimeStart;}("Time max = %g days = %g secn", TimeMax,TimeMax*24*60*60);

В результате выполнения программы получаем файл «углы.csv» со следующими данными:

·Количество проделанных итераций (равно количеству точек положения КА в файле «ОРБИТА.txt»);

·Дата и время в формате UTC: year-month-day hours:minutes:seconds;

·Угол места КА над горизонтом Луны;

·Угол возвышения рельефа.

Вся информация сохраняется в текстовом формате csv, предназначенном для представления табличных данных, для того, чтобы полученные данные было легче обработать, например, в программе Excel (рис. 28).

Рис. 28. Данные в файле «углы.csv», полученные после выполнения программы

Также программа выводит на экран значение максимального периода времени в днях и секундах, в течение которого виден наблюдаемый КА с заданной точки на поверхности Луны (рис. 29). Эта информация является очень полезной для выбора пункта наблюдения.

Рис. 29. Вывод на экран максимального периода времени, в течение которого виден КА.

На рис. 28 и рис. 29 представлены результаты выполнения программы для заданной точки на поверхности Луны с координатами по широте 71 ̊ и долготе 8 ̊ .

Если построить график зависимости угла места КА и угла возвышения лунного рельефа от времени (рис. 30), то можно увидеть, что начиная с определенного момента времени все время меньше . То есть КА постоянно находится вне зоны видимости, так как обзор на него перекрывается возвышенностями рельефа Луны. Это объясняется тем, что выбранная точка с координатами по широте 71 ̊ и долготе 8 ̊ находится внутри древнего кратера Барроу, расположенного в северной приполярной части Луны (рис.31). Координаты центра кратера равны 71,28 ̊ с.ш. и 7,59 ̊ в.д., диаметр — 94 км, а максимальная глубина составляет 2,38 км.

Рис. 30. График зависимости КА и высот от времени для точки с координатами 71 ̊ с.ш. и 8 ̊ в.д.

Рис. 30. Кратер Барроу в северной части видимой стороны Луны.

Рассмотрим на поверхности Луны другую точку наблюдения, находящуюся на возвышенности на севере в районе кратера Белькович (рис.31). То есть точку с координатами по широте 65 ̊ и долготе 87 ̊ .

Рис. 31. Точка с координатами 65 ̊ с.ш. и 87 ̊ в.д. вблизи кратера Белькович.

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

Рис. 32. График зависимости КА и высот от времени для точки с координатами 65 ̊ с.ш. и 87 ̊ в.д.

Отсюда следует, что из заданной точки на поверхности Луны КА будет виден хорошо. Так же это подтверждается значением максимального периода времени, в течение которого можно наблюдать КА. Для точки с координатами 65 ̊ с.ш. и 87 ̊ в.д. он равен 177045 секунд, или же примерно 2 дня (рис.33). В то время как, для точки внутри кратера Барроу максимальный период был равен чуть больше 1 дня (рис. 29).

Рис. 33. Вывод максимального периода времени, в течение которого виден КА, для точки с координатами 65 ̊ с.ш. и 87 ̊ в.д

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

Заключение

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

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

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

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

Библиографический список

1.P. K. Seidelmann (Chair), V. K. Abalakin, M. Bursa, M. E. Davies, C. de Bergh, J. H. Lieske, J. Oberst, J. L. Simon, E. M. Standish, P. Stooke, P. C. Thomas : Report of the IAU/IAG Working Group on Cartographic Coordinates and Rotational Elements of the Planets and Satellites: 2000.

.Joshua B. Hopkins, Lockheed Martin Corporation, United States: William Pratt, Caley Buxton, Selena Hall, Andrew Scott, Lockheed Martin Corporation, United States: Robert Farquhar, David Dunham, KinetX, Inc, United States: Proposed orbits and trajectories for human missions to the Earth-Moon L2 region.

.Dunham, D.W., Farquhar, R.W., Eismont, N., Chumachenko, E.N., Aksenov, S.A., Genova A., Horsewood J., Furfaro R., Kidd J: Using lunar swingbys and libration-point orbits to extend human exploration to interplanetary destinations, in: Proceedings of the International Astronautical Congress, IAC, 2013, 64th International Astronautical Congress 2013, IAC 2013; Beijing; China; 23 September 2013 through 27 September 2013 Vol. 2. International Astronautical Federation, 2013. P. 1932-1941.

.Подбельский В.В. Практикум по программированию на языке Си: Учебное пособие. М.: Финансы и Статистика, 2004.

.Бакулин П.И., Кононович Э.В., Мороз В.И. Курс общей астрономии: Учебник для студентов высших учебных заведений специальности «Астрономия». 4-е изд. М.: Наука, 1977.-544 с.

.Абалакин В.К., Аксенов Е.П., Гребеников Е.А., Демин В.Г., Рябов Ю.А., Справочное руководство по небесной механике и астродинамике. М.: Наука, 1976.

.Выгодский М.Я., Справочник по высшей математике. М.: Астрель, 2006.

.Себехей В., Теория орбит: Ограниченная задача трех тел. М.: Наука, 1982.

.Левантовский В.И., Механика космического полета. М.: Наука, 1980.

Приложение

/*Нахождение угла места для точки КА для данного момента времени. Считыванм из файла "ОРБИТА.txt" координаты, записываем все результаты в "угол_места.txt"*/

#define _CRT_SECURE_NO_DEPRECATE

#include <time.h>

#include <stdio.h>

#include <cstdio>

#include <conio.h>

#include <locale.h>

#include <math.h>

#include <stdlib.h>

#include <string.h>struct{coord[3];JDtime;

}vector;struct{teta, fi, r;

}SPoint;

//Получение Юлианской даты по Дате и Времени

double DateTimeToJulDate(int Year, int Month, int Day, int Hour, int Min, double Sec){MjdMidnight;

//double FracOfDay;b;(Month <= 2){ Month += 12; —Year; }((10000L * Year + 100L * Month + Day) <= 15821004L)= -2 + ((Year + 4716) / 4) — 1179;= (Year / 400) — (Year / 100) + (Year / 4);= 365L * Year — 679004 + b + int(30.6001*(Month + 1)) + Day;MjdMidnight + ((Hour + Min / 60.0 + Sec / 3600.0) / 24.0) + 2400000.5;}

//Нахождение длины вектора по координатам.

double length(double * vect)

{double l;= sqrt(pow(vect[0], 2) + pow(vect[1], 2) + pow(vect[2], 2));

return l;}

//Нахождение угла между векторами

double angle(double *v1, double *v2)

{double ang;= acos((v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]) / (length(v1)*length(v2))); ang;}

//Расчет угла места КА над горизонтом относительно точки point

//point[3] — массив координат x y z точки на поверхности Луны

//KA[3] — массив координат x y z КА

double elevation_angle(double * point, double * KA)

{double Rl,Ra,r, betta, b, gamma, g, alpha;R[3];rtd = 57.295779513082320876798161804285;= length(point);//длина вектора Rl

Ra = length(KA);//длина вектора Ra

//Вектор R[0] = KA[0] — point[0];[1] = KA[1] — point[1];[2] = KA[2] — point[2]; = length(R);//длина вектора R

betta = angle(point, KA);//угол между векторами Rl и Ra= angle(R, KA);//угол между векторами R и Ra

//Перевод из радиан в градусы.

b = betta*rtd;= gamma*rtd;= 90 — g — b;alpha;}

//умножение строки на матрицу* str_on_matrix(double m1[3], double m2[3][3])

{double * result = (double *)malloc(3 * sizeof(double));[0] = m1[0] * m2[0][0] + m1[1] * m2[1][0] + m1[2] * m2[2][0];[1] = m1[0] * m2[0][1] + m1[1] * m2[1][1] + m1[2] * m2[2][1];[2] = m1[0] * m2[0][2] + m1[1] * m2[1][2] + m1[2] * m2[2][2];

return result;}

// Перевод из СК SG в СК J2000.

double * SG_to_J(double Rsg[], double dt)

{double * Rj = (double *)malloc(3 * sizeof(double));i, j;//i-строка, j-столбецT, alph, bet, w, e[13], M1[3][3], M2[3][3];

int p = 10;dtr = 0.017453292519943295769236905555556; // коэффициент перевода градусов в радианыgaml = 13.17635815; // угл. скорость вращения Луны в градусах.

dt = dt — 2415020.0;[0] = (125.045 — 0.0529921*dt)*dtr;[1] = (250.089 — 0.1059842*dt)*dtr;[2] = (260.008 + 13.0120009*dt)*dtr;[3] = (176.625 + 13.3407154*dt)*dtr;[4] = (357.529 + 0.9856003*dt)*dtr;[5] = (311.589 + 26.4057084*dt)*dtr;[6] = (134.963 + 13.0649930*dt)*dtr;[7] = (276.617 + 0.3287146*dt)*dtr;[8] = (34.226 + 1.7484877*dt)*dtr;[9] = (15.134 — 0.1589763*dt)*dtr;[10] = (119.743 + 0.0036096*dt)*dtr;[11] = (239.961 + 0.1643573*dt)*dtr;[12] = (25.053 + 12.9590088*dt)*dtr;= dt / 36525;= (269.9949 + 0.0031*T — 3.8787*sin(e[0]) — 0.1204*sin(e[1]) + 0.07*sin(e[2]) — 0.0172*sin(e[3]) + 0.0072*sin(e[5]) — 0.0052*sin(e[9]) — 0.0043*sin(e[12]))*dtr;= (66.5392 + 0.013*T + 1.5419*cos(e[0]) + 0.0239*cos(e[1]) — 0.0278*cos(e[2]) + 0.0068*cos(e[3]) — 0.0029*cos(e[5]) + 0.0009*cos(e[6]) + 0.0008*cos(e[9]) — 0.0009*cos(e[12]))*dtr;= (38.3213 + gaml*dt + gaml + dt — 1.4*(p, (-12))*pow(dt, 2) + 3.561*sin(e[0]) + 0.1208*sin(e[1]) — 0.0642*sin(e[2]) + 0.0158*sin(e[3]) + 0.0252*sin(e[4]) — 0.0066*sin(e[5]) — 0.0047*sin(e[6]) — 0.0046*sin(e[7]) + 0.0028*sin(e[8]) + 0.0052*sin(e[9]) + 0.004*sin(e[10]) + 0.0019*sin(e[11]) — 0.0044*sin(e[12]))*dtr;[0][0] = -sin(alph)*cos(w) — cos(alph)*sin(bet)*sin(w);[0][1] = cos(alph)*cos(w) — sin(alph)*sin(bet)*sin(w);[0][2] = sin(w)*cos(bet);[1][0] = sin(alph)*sin(w) — cos(alph)*sin(bet)*cos(w);[1][1] = -cos(alph)*sin(w) — sin(alph)*sin(bet)*cos(w);[1][2] = cos(w)*cos(bet);[2][0] = cos(alph)*cos(bet);[2][1] = sin(alph)*cos(bet);[2][2] = sin(bet);

//Траспонирование матрицы(i = 0; i < 3; i++)

{for (j = 0; j < 3; j++)

{M2[i][j] = M1[j][i];}

}= str_on_matrix(Rsg, M2);Rj;}read_line(FILE * STREAM, char *line){i;(fgets(line, 256, STREAM) == NULL)(-1);= strlen(line);(i>0 && (line[i — 1] == ‘r’ || line[i — 1] == ‘n’)){-; line[i] = ‘’;}i;}

/*Перевод из сферических в декартовы*/

void S_to_D(double sph[], double dec[])

{//широта долгота радиус[0] = sph[2] * sin(sph[0])*sin(sph[1]);[1] = sph[2] * cos(sph[0])*sin(sph[1]);[2] = sph[2] * cos(sph[1]);}main()

{(LC_ALL, "RUSSIAN");sch = 0, Nstr = 0;// sch-счетчик для case; Nstr-номер строкиstring[300], *pch;

//Объявление переменныхR_Moon = 1738.57;// радиус Луны

double point[3] = { 175, 74, R_Moon };// точка на поверхности Луны

double decpoint[3] = { 0 };*p;x_sc, y_sc, z_sc;//координаты КА в СК J200year, month, day, hour, min, sec;KA[3];line[1024];JD = 0;//Юлианская датаugol = 0;//угол местаmoon[3];

int npl = 10;Rl, Ra;//радиус Луны и расстояние до КА* shir, * dolg;//диапозоны по широте и долготеflag[10];//флаг. Если КА виден, пишет "да", если нет — "нет"* f;//f-файл, куда будем записывать

if ((f = fopen("D:\image\Катя\Documents\Visual Studio 2013\Projects\vidimost_bes_reliefa\угол_места.txt", "wt")) == 0)

{printf("Ошибка!");0;}(f, " tttt point tttttttttt KA J200 tttttt time tttt угол места tt виден? n");(f, " point x t point y tt point z tt KA J2000 x t KA J2000 y t KA J200 z n");

//Переводим из сферических в декартовы_to_D(point, decpoint);* o;//o-файл, откуда читаемл

if ((o = fopen("D:\image\Катя\Documents\Visual Studio 2013\Projects\Проект 8(с файлами)\ОРБИТА.txt", "rt")) == 0)

{printf("Ошибка!");0;}(string, 300, o);(!feof(o))

{fgets(string, 300, o);(string != NULL)

{pch = strtok(string, " ");[0] = atof(pch);(pch != NULL){(sch<5) pch = strtok(NULL, " ");((sch>5) && (sch<9)) pch = strtok(NULL, " .");(sch >= 9) pch = strtok(NULL, ":");(sch) {0: KA[1] = atof(pch); break;1: KA[2] = atof(pch); break;6: day = atof(pch); break;7: month = atof(pch); break;8: year = atof(pch); break;9: hour = atof(pch); break;10: min = atof(pch); break;11: sec = atof(pch); break;}

sch++;}= 0;

//Переводим из UTC в Юлианскую дату

JD = DateTimeToJulDate(year, month, day, hour, min, sec);

//Переводим из SG в СК J2000

p = SG_to_J(decpoint, JD);

//Находим угол места= elevation_angle(p, KA);

//Определяем виден КА или нет.

if (ugol > 0)

{flag[0] = ‘д’;[1] = ‘а’;[2] = ‘ ‘;}

{flag[0] = ‘н’;[1] = ‘е’;[2] = ‘т’;}(f, "%lft %lft %lft %15.5lf %15.5lf %15.5lftt %d.%d.%d %d:%d:%dt %15.5lf t %5.5c%c%c n", point[0], point[1], point[2], KA[0], KA[1], KA[2], year, month, day, hour, min, sec, ugol, flag[0], flag[1], flag[2]);++;}}(f);(o); 0;}

Файл «main.cpp».

/*нахождение угла места КА и угла места высот и максимального времени виденья для заданной точки*/

#define _CRT_SECURE_NO_DEPRECATE

#include <time.h>

#include <stdio.h>

#include <cstdio>

#include <conio.h>

#include <locale.h>

#include <math.h>

#include <stdlib.h>

#include <string.h>

#include "orbita.h"

#include "moon.h"

#include "dop.h"

#define FILE_ORBIT "D:\image\Катя\Documents\Visual Studio 2013\Projects\vidimost_s_reliefom\ОРБИТА.txt"

#define FILE_MOON "D:\image\Катя\Desktop\ДИПЛОМ\Moon\moon.txt"

#define FILE_ANGLES "D:\image\Катя\Documents\Visual Studio 2013\Projects\vidimost_s_reliefom\углы.csv"

//"D:\image\Катя\Documents\Visual Studio 2013\Projects\vidimost_s_reliefom\видимость.txt"

#define MoonRadius 1738.57

#define AngleStep 0.5e0main()

{int i, j, k;_point*KA;NKA;iO; // index of Orbit point_height*MoonHeight;*angle_KA;*angle_Moon;

//Объявление переменныхpoint_S[3] = { 65, 87, MoonRadius };// точка на поверхности Луны P'(lat,lon,R_Moon)point_J[3];MoonCenterSG[3];MoonCenterJ[3];vecOP[3],vecOK[3],vecN[3],vecOS[3];pointS_J[3],pointS_S[3];MatrixRotation[3][3];AngleNStep;a;_read(FILE_ORBIT, &KA, &NKA);_height_read(FILE_MOON, 60, 90, 0, 360, -15520, 15520, &MoonHeight);[0] = 0.e0;[1] = 0.e0;[2] = 0.e0;_KA = (double*)calloc(NKA, sizeof(double));_Moon = (double*)calloc(NKA, sizeof(double));(iO = 0; iO < NKA; iO++){_S[2] = MoonRadius + moon_height_get(MoonHeight, point_S[0], point_S[1]);_to_J(point_S, point_J, KA[iO].time);_to_J(MoonCenterSG, MoonCenterJ, KA[iO].time);(k = 0; k < 3; k++) vecOP[k] = point_J[k] — MoonCenterJ[k];(k = 0; k < 3; k++) vecOK[k] = KA[iO].coords[k] — MoonCenterJ[k];(fabs(angle(vecOP, vecOK)) < 1.e-10){

//TODO видно, угол 90_KA[iO] = 90.e0;_Moon[iO] = 0.e0;;}_KA[iO] = elevation_angle(vecOP, vecOK);

(vecOP, vecOK, vecN);_norm(vecN);= (int)(angle(vecOP, vecOK) / AngleStep + 1);(i = 1; i <= AngleNStep; i++){_rotation(vecN, i*AngleStep, MatrixRotation);_on_matrix(vecOP, MatrixRotation, vecOS);(k = 0; k < 3; k++) pointS_J[k] = MoonCenterJ[k] + vecOS[k];_to_S(pointS_J, pointS_S,KA[iO].time);_S[2] = MoonRadius + moon_height_get(MoonHeight, pointS_S[0], pointS_S[1]);_to_J(pointS_S, pointS_J, KA[iO].time);(k = 0; k < 3; k++) vecOS[k] = pointS_J[k] — MoonCenterJ[k];= elevation_angle(vecOP, vecOS);(i == 1){_Moon[iO] = a;

}else{(angle_Moon[iO] < a) angle_Moon[iO] = a;}}("i=%dtangle_KA=%gtangle_Moon=%gn", iO, angle_KA[iO], angle_Moon[iO]);}(LC_ALL, "RUSSIAN");* fd;((fd = fopen(FILE_ANGLES, "wt")) == 0){("Can’t open file Angels to write ‘%s’n",FILE_ANGLES);();0;}(fd, "Итерация;Дата;Время;Угол КА;Угол рельефаn");(i = 0; i < NKA; i++){(fd, "%d;%d-%d-%d;%d:%d:%d;%g;%gn", i,[i].time[0], KA[i].time[1], KA[i].time[2],[i].time[3], KA[i].time[4], KA[i].time[5],_KA[i], angle_Moon[i]);}(fd);TimeMax,Time,TimeStart,TimeFinish;= 0.e0;;= -1.e0;(i = 0; i < NKA; i++){= DateTimeToJulDate(KA[i].time[0], KA[i].time[1], KA[i].time[2], KA[i].time[3], KA[i].time[4], KA[i].time[5]);(angle_KA[i] > angle_Moon[i]){(TimeStart < 0){= Time;= Time;

}else{= Time;}

}else{(TimeStart>0){(TimeFinish — TimeStart > TimeMax)= TimeFinish — TimeStart;

}else{//do nothing}

}}(TimeStart > 0){(TimeFinish — TimeStart > TimeMax)= TimeFinish — TimeStart;}("Time max = %g days = %g secn", TimeMax,TimeMax*24*60*60);

getch();0;}

Файл с вспомогательными программами «dop.cpp».

#define _CRT_SECURE_NO_DEPRECATE

#include <time.h>

#include <stdio.h>

#include <cstdio>

#include <conio.h>

#include <locale.h>

#include <math.h>

#include <stdlib.h>

#include <string.h>

#include "dop.h"struct{coord[3];JDtime;

}vector;struct{teta, fi, r;

//Получение Юлианской даты по Дате и Времени

double DateTimeToJulDate(int Year, int Month, int Day, int Hour, int Min, double Sec){ long MjdMidnight;b;(Month <= 2){ Month += 12; —Year; }((10000L * Year + 100L * Month + Day) <= 15821004L)= -2 + ((Year + 4716) / 4) — 1179;= (Year / 400) — (Year / 100) + (Year / 4);= 365L * Year — 679004 + b + int(30.6001*(Month + 1)) + Day;MjdMidnight + ((Hour + Min / 60.0 + Sec / 3600.0) / 24.0) + 2400000.5;}

//Нахождение длины вектора по координатам.

double length(double * vect)

{double l;= sqrt(pow(vect[0], 2) + pow(vect[1], 2) + pow(vect[2], 2));l;}

//Скалярное произведение векторовskalyr(double *v1, double *v2){v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];}

//Нахождение угла между векторамиangle(double *v1, double *v2){RTD*acos(skalyr(v1, v2) / (length(v1)*length(v2)));}

//Векторное произведение векторовmvector(double v1[3], double v2[3],double v3[3]){[0] = v1[1] * v2[2] — v1[2] * v2[1];[1] = v1[2] * v2[0] — v1[0] * v2[2];

v3[2] = v1[0] * v2[1] — v1[1] * v2[0];}

//Нормировать вектор

void vec_norm(double v[3]){len;= length(v);[0] = v[0] / len;[1] = v[1] / len;

v[2] = v[2] / len;}

//Получить матрицу поворота вокруг вектора v на угол angle

void matrix_rotation(double v[3], double angle, double M[3][3]){cosa, sina;x, y, z;= cos(angle*DTR);= sin(angle*DTR);= v[0];= v[1];= v[2];[0][0] = cosa + (1 — cosa)*x*x;[0][1] = (1 — cosa)*x*y — sina*z;[0][2] = (1 — cosa)*x*z + sina*y;[1][0] = (1-cosa)*y*x + sina*z;[1][1] = cosa + (1 -cosa)*y*y;[1][2] = (1 — cosa)*y*z — sina*x;[2][0] = (1-cosa)*z*x — sina*y;[2][1] = (1 — cosa)*z*y + sina*x;[2][2] = cosa+(1 — cosa)*z*z;}

//Расчет угла места

// OP — векртор в J2000 из центра Луны в точку на поверхности Луны;

// OK — вектор в J2000 из центра Луны в координаты КА на гало-орбите.

double elevation_angle(double * OP, double * OK)

{int k;lenOP, lenOK, lenPK, betta, gamma, alpha;

double PK[3];// вектор от точки на поверхности луны в положение КА= length(OP);//длина вектора OP= length(OK);//длина вектора OK

for (k = 0; k < 3; k++) PK[k] = OK[k] — OP[k];

lenPK = length(PK);//длина вектора PK= angle(OP, OK);//угол между векторами OP и OK= angle(PK, OK);//угол между векторами PK и OK

alpha = 90.e0 — gamma — betta;alpha;}

//умножение вектора на матрицуvec_on_matrix(double m1[3], double m2[3][3], double m3[3]){[0] = m1[0] * m2[0][0] + m1[1] * m2[0][1] + m1[2] * m2[0][2];[1] = m1[0] * m2[1][0] + m1[1] * m2[1][1] + m1[2] * m2[1][2];

m3[2] = m1[0] * m2[2][0] + m1[1] * m2[2][1] + m1[2] * m2[2][2];

}

//Траспонирование матрицы

void trans_matrix(double M[3][3]){M2[3][3];i, j;(i = 0; i < 3; i++)(j = 0; j < 3; j++)[i][j] = M[j][i];(i = 0; i < 3; i++)(j = 0; j < 3; j++)[i][j] = M2[i][j];}martix_J_to_SG(double Rj[3], double M[3][3], int time[6]){dt, T, alph, bet, w, e[13];p = 10;gaml = 13.17635815;= DateTimeToJulDate(time[0], time[1], time[2], time[3], time[4], time[5]) — 2415020.0;[0] = (125.045 — 0.0529921*dt)*DTR;[1] = (250.089 — 0.1059842*dt)*DTR;[2] = (260.008 + 13.0120009*dt)*DTR;[3] = (176.625 + 13.3407154*dt)*DTR;[4] = (357.529 + 0.9856003*dt)*DTR;[5] = (311.589 + 26.4057084*dt)*DTR;[6] = (134.963 + 13.0649930*dt)*DTR;[7] = (276.617 + 0.3287146*dt)*DTR;[8] = (34.226 + 1.7484877*dt)*DTR;[9] = (15.134 — 0.1589763*dt)*DTR;[10] = (119.743 + 0.0036096*dt)*DTR;[11] = (239.961 + 0.1643573*dt)*DTR;[12] = (25.053 + 12.9590088*dt)*DTR;= dt / 36525;= (269.9949 + 0.0031*T — 3.8787*sin(e[0]) — 0.1204*sin(e[1]) + 0.07*sin(e[2]) — 0.0172*sin(e[3]) + 0.0072*sin(e[5]) — 0.0052*sin(e[9]) — 0.0043*sin(e[12]))*DTR;= (66.5392 + 0.013*T + 1.5419*cos(e[0]) + 0.0239*cos(e[1]) — 0.0278*cos(e[2]) + 0.0068*cos(e[3]) — 0.0029*cos(e[5]) + 0.0009*cos(e[6]) + 0.0008*cos(e[9]) — 0.0009*cos(e[12]))*DTR;= (38.3213 + gaml*dt + gaml + dt — 1.4*(p, (-12))*pow(dt, 2) + 3.561*sin(e[0]) + 0.1208*sin(e[1]) — 0.0642*sin(e[2]) + 0.0158*sin(e[3]) + 0.0252*sin(e[4]) — 0.0066*sin(e[5]) — 0.0047*sin(e[6]) — 0.0046*sin(e[7]) + 0.0028*sin(e[8]) + 0.0052*sin(e[9]) + 0.004*sin(e[10]) + 0.0019*sin(e[11]) — 0.0044*sin(e[12]))*DTR;[0][0] = -sin(alph)*cos(w) — cos(alph)*sin(bet)*sin(w);[0][1] = cos(alph)*cos(w) — sin(alph)*sin(bet)*sin(w);[0][2] = sin(w)*cos(bet);[1][0] = sin(alph)*sin(w) — cos(alph)*sin(bet)*cos(w);[1][1] = -cos(alph)*sin(w) — sin(alph)*sin(bet)*cos(w);[1][2] = cos(w)*cos(bet);[2][0] = cos(alph)*cos(bet);[2][1] = sin(alph)*cos(bet);

M[2][2] = sin(bet);}

/*Перевод из СК SG(декартовых) в сферические*/

void SG_to_S(double dec[3], double sph[3])

{double x, y, z;= dec[0];= dec[1];= dec[2];[0] = RTD*atan2(x, y); // lat[1] = RTD*atan2(sqrt(x*x + y*y),z); // lon[2] = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); }J_to_S(double Rj[3], double Rs[3], int time[6]){Rsg[3];_to_SG(Rj, Rsg, time);_to_S(Rsg, Rs);}S_to_J(double Rs[3], double Rj[3], int time[6]){Rsg[3];_to_SG(Rs, Rsg);_to_J(Rsg, Rj, time);}read_line(FILE * STREAM, char *line){i;(fgets(line, 256, STREAM) == NULL)(-1);= strlen(line);(i>0 && (line[i — 1] == ‘r’ || line[i — 1] == ‘n’)){-; line[i] = ‘’;}i;}

Файл «moon.cpp».

#define _CRT_SECURE_NO_DEPRECATE

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include "moon.h"moon_height_read(char*file,lat_min, double lat_max,lon_min, double lon_max,lat_N, int lon_N, moon_height**_Height){

FILE * fd;//o-файл, откуда читаем

int i, j, i2,j2;h;_height*Height;= (moon_height*)calloc(1, sizeof(moon_height));>lat_min = lat_min;>lat_max = lat_max;>lon_min = lon_min;>lon_max = lon_max;>lat_N = abs(lat_N);>lon_N = abs(lon_N);>height = (double**)calloc(Height->lat_N , sizeof(double*));(i = 0; i < Height->lat_N;i++)>height[i] = (double*)calloc(Height->lon_N, sizeof(double));((fd = fopen(file, "rt")) == 0){("Can’t read moon file ‘%s’n", file);0;}(i = 0; i < Height->lat_N; i++){(j = 0; j < Height->lon_N; j++){= (lat_N < 0) ? (Height->lat_N — 1 — i) : i;= (lon_N < 0) ? (Height->lon_N — 1 — j) : j;(fd, "%lf", &h);>height[i2][j2] = h;}}

*_Height = Height;("Readed file moon. Read %dx%d from lat[%lf, %lf] x lon[%lf, %lf].n", Height->lat_N,Height->lon_N,>lat_min,Height->lat_max,Height->lon_min,Height->lon_max);(fd);0;}moon_height_get(moon_height*Height, double lat, double lon){iLat, iLon;= (int)(Height->lat_N * (lat — Height->lat_min) / (Height->lat_max — Height->lat_min));(iLat < 0) iLat = 0;(iLat >= Height->lat_N) iLat = Height->lat_N — 1;= (int)(Height->lon_N * (lon — Height->lon_min) / (Height->lon_max — Height->lon_min));(iLon < 0) iLon = 0;(iLon >= Height->lon_N) iLon = Height->lon_N — 1;Height->height[iLat][iLon];}

Файл «orbita.cpp».

#define _CRT_SECURE_NO_DEPRECATE

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include "orbita.h"

// Считать орбиту спутникаorbit_read(char*file, orbita_point**_KA,int*_NKA){ * o;//o-файл, откуда читаем

orbita_point*KA;NKA;NKA_allocated;sch = 0, Nstr = 0;// sch-счетчик для case; Nstr-номер строкиstring[301], *pch;((o = fopen(file, "rt")) == 0){("Can’t read orbit file ‘%s’n", file);0;}= 0;= 0;_allocated = 0;

//считываем координаты KA в СК J2000 и время из "ОРБИТА.txt"

// пропускаем заголовки

fgets(string, 300, o);

while (!feof(o)){(string, 300, o);(string != NULL){(NKA + 1 >= NKA_allocated){(NKA_allocated == 0){_allocated = 16;= (orbita_point*)calloc(NKA_allocated, sizeof(orbita_point));}{(NKA + 1 >= NKA_allocated) NKA_allocated = NKA_allocated * 2;= (orbita_point*)realloc(KA, NKA_allocated * sizeof(orbita_point));}}= strtok(string, " ");[NKA].coords[0] = atof(pch); // X(pch != NULL){(sch<5) pch = strtok(NULL, " ");((sch>5) && (sch < 9)) pch = strtok(NULL, " .");(sch >= 9) pch = strtok(NULL, ":");(sch) {0: KA[NKA].coords[1] = atof(pch); break; // Y1: KA[NKA].coords[2] = atof(pch); break; // Z6: KA[NKA].time[2] = atoi(pch); break; // Day7: KA[NKA].time[1] = atoi(pch); break; // Month8: KA[NKA].time[0] = atoi(pch); break; // Year9: KA[NKA].time[3] = atoi(pch); break; // Hour10: KA[NKA].time[4] = atoi(pch); break; // Min11: KA[NKA].time[5] = atoi(pch); break; // Sec

}++;}= 0;++;++;}}("Readed file orbit. Count of points %d.n", NKA);

*_KA = KA;

*_NKA = NKA;0;}