Тангаж, рысканье в норме. Выводим на орбиту ИСЗ.

 

Мои поздравления по случаю состоявшейся вчера 54-летней годовщины начала космической эры. В честь этого знаменательного события, которое отмечало все прогрессивное человечество, давайте тоже что-нибудь запустим. Начнем с момента, когда изделие вышло на орбиту, сориентировалось и приступило к штатной задаче мониторинга земной поверхности на предмет несанкционированных очагов возгорания, раннего обнаружения нефтяных пятен в районах, к тому не предназначенных, и прочей экологии. При помощи появившейся в SQL Server Denali поддержки FULLGLOBE проложим на глобусе его наземную трассу, чтобы понять, где конкретно оно в данный момент летит.

Для начала понадобится привязаться к какой-нибудь стандартной геодезической системе координат. Для операций планетного масштаба распространена геоцентрическая эллипсоидальная система WGS-84 (World Geodetic System), разработанная Военно-картографическим ведомством Министерства обороны США, последняя, 3-я редакция которой опубликована в 1997 г., она же используется в GPS, которая в SQL Server значится под SRID = 4326. Эллипсоид - тело, полученное вращением эллипса вокруг его в данном случае малой оси (Земля имеет неправильную сплюснутую с полюсов форму), размеры осей которого подобраны так, чтобы среднеквадратичное отклонение от поверхности Земли (геоида) было минимально либо по всей поверхности в целом, либо для заданной территории, например, России.

 

Топографической службой Вооруженных сил РФ разработана система ПЗ-90 на основе эллипсоида, построенного группой под руководством директора ЦНИИГАиК Ф.Н.Красовского и проф.А.А.Изотова, которая используется, в частности в ГЛОНАСС. Постановлением Правительства от 28 июля 2000 г. N 568 "Об установлении единых государственных систем координат" в России приняты в качестве стандарта геоцентрическая система координат "Параметры Земли 1990 года" (ПЗ-90) для использования в целях геодезического обеспечения орбитальных полетов и решения навигационных задач и система геодезических координат 1995 года (СК-95) для использования при осуществлении геодезических и картографических работ, начиная с 1 июля 2002 г. В SQL Server она значится под SRID = 4200.

 

Из-за того, что большинство карт, которые мне удалось найти в бесплатном доступе, привязаны к WGS-84, использовать в данном посте отечественную разработку не удастся :(

 

Рис.1

 

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

 

use tempdb

 

if object_id('Константы_WGS84', 'U') is not null drop table Константы_WGS84

create table Константы_WGS84 (name varchar(5) primary key, value float, description nvarchar(300))

insert Константы_WGS84 (name, value, description) values ('Ra', 6378137, 'Большая полуось эллипсоида Земли, м'),

                                                                                    ('f', 298.257223563, 'Сплюснутость'),

                                                                                    ('Wz', 7.2921158553e-5, 'Угловая скорость вращения Земли, рад/с'),

                                                                                    ('Gz', 398600.5e9, 'Геоцентрическая гравитационная постоянная, м3/с2')

go

if object_id('Константа_WGS84', 'FN') is not null drop function Константа_WGS84

go

create function Константа_WGS84(@name varchar(5)) returns float as

begin

 return (select value from Константы_WGS84 where name = @name)

end

go

insert Константы_WGS84 (name, value, description) values ('Rb', dbo.Константа_WGS84('Ra') * (1 - 1 / dbo.Константа_WGS84('f')), 'Малая полуось эллипсоида Земли, м'),

                                                                                    ('eZ', sqrt(dbo.Константа_WGS84('Ra') * dbo.Константа_WGS84('Ra') - dbo.Константа_WGS84('Rb') * dbo.Константа_WGS84('Rb')) / dbo.Константа_WGS84('Ra'), 'Эксцентриситет')

Рис.2

 

Частично их можно найти тут:

 

select * from sys.spatial_reference_systems where spatial_reference_id = 4326

Рис.3

 

остальные - поиском в Интернете.

 

Сплюснутость по науке называется коэффициентом сжатия и считается как большая полуось / (большая полуось - малая полуось). Геоцентрическая гравитационная постоянная считается как глобальная гравитационная постоянная (это которую в 6 классе на крутильных весах определял Кавендиш) умножить на массу Земли с учетом атмосферы.

 

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

 

Для описания невозмущенного движения достаточно трех законов Кеплера, которые все тоже проходили в школе то ли по физике, то ли по астрономии ближе к 10 кл. Первый из них гласит, что в поле тяготения тела движутся по коническим сечениям, в частности, по эллипсу, в одном из фокусов которого расположен притягивающий центр, если, конечно, у них не хватает сил оторваться по параболе. Не говоря про гиперболу, по которой вообще можно оторваться так, что мама не горюй. Как известно из курса военной подготовки, эллипс - это круг, вписанный в прямоугольник. Это отточенная, но не вполне корректная формулировка. Если плющить круг, получится овал. Форму эллипса можно задать, например, при помощи большой полуоси a и эксцентриситета e.

Теперь этот эллипс надо сориентировать вокруг Земли. Ставим его сперва торчком перпендикулярно экватору, чтобы эллипс располагался в плоскости сечения нулевого меридиана и линия апсид (длинная ось, соединяющая перигей с апогеем) совпадала с осью вращения Земли. Повернем плоскость эллипса вокруг линии апсид на угол Omega. Он называется долготой восходящего узла. Точка восходящего узла задает также направление движения по орбите. В ней спутник переходит из полупространства южного полушария Земли в северное. Проложим новую ось как результат пересечения плоскости эллипса с плоскос��ью экватора. Она называется линией узлов. Отпускаем вертикаль. Эллипс начнет свободно покачиваться вокруг линии узлов. Угол i между плоскостью экватора и плоскостью эллипса называется наклонением орбиты - см., напр., рис. Кстати, это неточный рис. Согласно военного устава центр Земли должен находиться в фокусе эллипса, а не в его центре. Будем считать, что на нем изображена круговая орбита. Зафиксируем плоскость эллипса и повернем его в ней вокруг фокуса. Угол w между линией узлов и направлением из фокуса на перигей, - см., напр., рис.4 - называется аргументом перигея.

 

Осталось задать положение спутника как точки на эллипсе и вычислить его от времени. Я помню, что это делается на основе остальных двух законов Кеплера, но как это делается, я не помню. Открываем любой популярный учебник по астрономии, параграф 3.10.2 - "Параметры и аномалии кеплеровской орбиты". Для задания положения точки на эллипсе будем откладывать угол v между перигеем и радиус-вектором, фокуса и упирающимся из фокуса в спутник - см. рис.3.14. Этот угол называется истинной аномалией. И еще осталось задать, где спутник находился в начальный момент времени. Это делается при помощи T0 - времени прохождения перигея. Замечательно. Я предлагаю держать параметры орбиты тоже в табличке, где ж еще? В качестве примера возьмем характеристики первого искусственного спутника Земли, выведенного на орбиту 4 октября 1957 г. Долготы восходящего узла не нашел, поставил по аналогии с последующими запусками, время прохождения перигея - наобум в 0.

 

if object_id('Параметры_орбиты', 'U') is not null drop table Параметры_орбиты

create table Параметры_орбиты (спутник tinyint identity, a float, e float, Omega float, i float, w float, T0 float)

insert Параметры_орбиты (a, e, Omega, i, w, T0) values (6955200, 0.05201, 28 * pi() / 180, 65.1 * pi() / 180, 58 * pi() / 180, 0)

go

if object_id('Параметр_орбиты', 'FN') is not null drop function Параметр_орбиты

go

create function Параметр_орбиты(@id tinyint, @name varchar(5)) returns float as

begin

 return (select value from Параметры_орбиты unpivot (Value for Field in (a, e, Omega, i, w, T0)) unpvt where спутник = @id and field = @name)

end

Рис.4

 

Из третьего закона Кеплера - квадраты периодов обращений спутников соотносятся, как кубы их больших полуосей, получается период обращения спутника по орбите

 

T = 2 * pi * sqrt( a * a * a / Gz)

 

Обратная периоду величина называется средним движением (3.51).

 

n = 2 * pi / T = sqrt (Gz / a / a / a)

 

По смыслу это действительно средняя угловая скорость.

Средняя скорость умножить на время называется средней аномалией:

 

M = n * (t - T0)

 

Из второго закона Кеплера - радиус-вектор из фокуса к спутнику за равные промежутки времени заметает равные площади - после интегрирования и преобразований получается равенство (3.63):

 

E = M + e * sin(E)

 

где Е - эксцентрическая аномалия (см. рис.3.3). Аналитически ее вычислить не удастся, зато по этой формуле она легко вычисляется итеративно с нужной точностью. За нулевое приближение можно принять E0 = M.

Эксцентрическая аномалия связана с истинной через формулу (3.60):

 

tg(v/2) = sqrt ((1 + е) / (1 - е)) * tg (E/2)

 

Таким образом, на каждый момент времени удастся получить v как функцию от t. Длина радиус-вектора от перифокуса до спутника согласно первому закону Кеплера есть уравнение эллипса в полярных координатах - см. (3.59):

 

r = p / (1 + e * cos(v))

 

где p - фокальный параметр = a * (1 - e * e).

 

Связанная с орбитой декартова система координат имеет ось х, направленную от перицентра к перигею, ось y, параллельную малой оси эллипса, и z, ортогональную плоскости орбиты. Все равно по ней ничего не растет: (r * cos(v), r * sin(v), 0) = (a * (cos(E) - e), a *sqrt(1 - e * e) * sin (E), 0) - см. (3.58). Чтобы перейти из орбитальной системы координат в земную, радиус-вектор нужно в обратном порядке отвернуть на углы -w вокруг z, -i вокруг x' и -Omega вокруг Z'' = Z, задающие ориентацию эллипса относительно Земли:

 

cos Omega	-sin Omega	0		1	0	0		cos w	-sin w	0
sin Omega	cos Omega	0		0	cos i	-sin i		sin w	cos w	0
0	0	1		0	sin i	cos i		0	0	1

то есть умножить резул��тирующую матрицу п��оизведения поворотов (3.68)

 

cos(Omega) * cos(w) - sin(Omega) * sin(w) * cos(i)	-cos(Omega) * sin(w) - sin(Omega) * cos(w) * cos(i)	sin(Omega) * sin(i)
sin(Omega) * cos(w) + cos(Omega) * sin(w) * cos(i)	-sin(Omega) * sin(w) + cos(Omega) * cos(w) * cos(i)	-cos(Omega) * sin(i)
sin(w) * sin(i)	cos(w) * sin(i)	cos(i)

 

на столбец координат радиус-вектора. Получатся трехмерные декартовы координаты в прямоугольной системе отсчета, центр которой совпадает с центром Земли, ось Z направлена на северный полюс, оси X и Y лежат в плоскости экватора. Х направлена на нулевой меридиан, Y ей ортогональна. Смещения прибавлять не нужно, т.к. центры обеих координатных систем совпадают.

Спроецировав точку с координатами (X, Y, Z) на поверхность приближающего Землю в выбранной геодезической системе, в данном случае WGC-84, референсного эллипсоида, мы получим географические координаты точки-следа от спутника на поверхности Земли. Параметры эллипсоида известны (Рис.3). Формулы пересчета между декартовыми координатами и широтой, долготой, высотой можно вывести самому, либо позаимствовать уже готовые, например, на сайте лаборатории ГАИШа - (2.11), из которых сразу получается долгота L = arctg(Y/X) и следом широта B. Нужно только помнить, что из-за вращения Земли долгота восходящего узла с каждым моментом времени уезжает в противоположном направлении. Прецессией, нутацией, мутацией и прочей флуктуацией для простоты пренебрегаем. Получается оч.красивая функция. Используя параметры орбиты по id спутника, она возвращает географическую точку на поверхности Земли (широта/долгота), над которой он в данный момент находится:

 

if object_id('ПроекцияПоложенияНаОрбите', 'FN') is not null drop function ПроекцияПоложенияНаОрбите

go

create function ПроекцияПоложенияНаОрбите(@номер_спутника tinyint, @t float) returns geography as

begin

 --Константы функции

 declare @Gz float = dbo.Константа_WGS84('Gz') --геоцентрическая гравитационная постоянная

 declare @Ra float = dbo.Константа_WGS84('Ra') --большая полуось эллипсоида Земли

 declare @eZ float = dbo.Константа_WGS84('eZ') --эксцентриситет эллипсоида Земли

 declare @a float = dbo.Параметр_орбиты(@номер_спутника, 'a') --большая полуось орбиты

 declare @e float = dbo.Параметр_орбиты(@номер_спутника, 'e') --эксцентриситет орбиты

 declare @Omega float = dbo.Параметр_орбиты(@номер_спутника, 'Omega') --долгота восходящего узла

 set @Omega -= dbo.Константа_WGS84('Wz') * @t --Земля вращается

 declare @i float = dbo.Параметр_орбиты(@номер_спутника, 'i') --наклонение

 declare @w float = dbo.Параметр_орбиты(@номер_спутника, 'w') --аргумент перигея

 --Находим истинную аномалию @v

 declare @n float = sqrt(@Gz / @a / @a / @a) --среднее движение

 declare @M float = @n * (@t - dbo.Параметр_орбиты(@номер_спутника, 'T0')) --средняя аномалия

 declare @Ea float, @Ea0 float = @M --эксцентрическая аномалия рассчитывается итеративно

 declare @eps float = 1e-8

 while 1 = 1 begin

  set @Ea = @M + @e * sin(@Ea0)

  if abs(@Ea - @Ea0) < @eps break

  set @Ea0 = @Ea

 end

 --Полярные координаты в СК орбиты (ось - линия апсид)

 --declare @v float = 2 * atan(sqrt((1 + @e) / (1 - @e)) * tan(@Ea / 2))

 --declare @r float = @a * (1 - @e * @e) / (1 + @e * cos(@v))

 declare @rcosv float = @a * (cos(@Ea) - @e)

 declare @rsinv float = @a * sqrt(1 - @e * @e) * sin(@Ea)

 --Преобразуем их в декартову геоцентрическую СК.

 --Поскольку по 3-й координате радиус-вектор имеет 0, 3-й столбец матрицы преобразования не понадобится.

 declare @a11 float = cos(@Omega) * cos(@w) - sin(@Omega) * cos(@i) * sin(@w)

 declare @a12 float = -cos(@Omega) * sin(@w) - sin(@Omega) * cos(@i) * cos(@w)

 declare @a21 float = sin(@Omega) * cos(@w) + cos(@Omega) * cos(@i) * sin(@w)

 declare @a22 float = -sin(@Omega) * sin(@w) + cos(@Omega) * cos(@i) * cos(@w)

 declare @a31 float = sin(@i) * sin(@w)

 declare @a32 float = sin(@i) * cos(@w)

 --В такие моменты особенно остро ощущаешь отсутствие в Т-SQL массивов.

 declare @X float = @a11 * @rcosv + @a12 * @rsinv

 declare @Y float = @a21 * @rcosv + @a22 * @rsinv

 declare @Z float = @a31 * @rcosv + @a32 * @rsinv

  --Переводим декартовы координаты в геодезические

 declare @B float /*широта*/, @L float /*долгота*/

 if @X = 0 and @Y = 0 --точка находится на полюсе, долгота может быть любая

  select @B = 0, @L = sign(@Z) * pi() / 2

 else begin

  set @L = atn2(@Y, @X) --долгота

  declare @B0 float = 0 --широту считаем итеративно аналогично эксцентрической аномалии

  declare @ro float --радиус кривизны эллипсоида в первом вертикале

  while 1 = 1 begin

   set @ro = @Ra / sqrt(1 - @eZ * @eZ * sin(@B0) * sin(@B0))

   set @B = atn2( (@Z + @ro * @eZ * @eZ * sin(@B0)), sqrt(@X * @X + @Y * @Y) )

   if abs(@B - @B0) < @eps break

   set @B0 = @B

  end

 end

 return geography::Point(@B * 180 / pi(), @L * 180 / pi(), 4326)

end

Рис.5

 

Нижеприводимый скрипт отмечает N промежуточных точек с равными интервалами в течение периода обращения. Для видимости они обозначены кругами радиусом в 100 км. Точки соединены отрезками в линию при помощи метода, который мы использовали в посте Превращение последовательности точек в геометрическую фигуру\Скрипт 8. Это позволяет нарисовать траекторию витка на поверхности Земли. Можно видеть, что из-за того, что вращение спутника происходит в направлении вращения Земли, ему не удается за виток описать полный круг на ее поверхности. Она от него как бы убегает.

 

declare @Period float = 2 * pi() * sqrt( power(dbo.Параметр_орбиты(1, 'a'), 3) / dbo.Константа_WGS84('Gz') )

declare @i int = 0, @N int = 8, @t0 float = 0, @t1 float = @Period

declare @g geography = geography::STGeomFromText('GEOMETRYCOLLECTION EMPTY', 4326), @b varbinary(max) = cast('' as varbinary), @p geography

while @i <= @N begin

 set @p = dbo.ПроекцияПоложенияНаОрбите(1, @t0 + (@t1 - @t0) * @i / @N)

 set @g = @g.STUnion(@p.STBuffer(100000) )

 set @b += substring(@p.STAsBinary(), 6, 16)

 set @i += 1

end

select @g union all select geography::STGeomFromWKB(0x01 + 0x02000000 + cast(reverse(cast(@N + 1 as binary(4))) as binary(4)) + @b, 4326)

 

Рис.6

 

Можно для интереса сделать @t1 = 2 * @Period и отрисовать траекторию второго витка:

 

Рис.7

 

Далее предполагается подложить на координатную сетку политическую карту мира, собрать орбитальную группировку типа GPS или ГЛОНАСС, покрутить ее во времени и посмотреть, насколько плотно она кого покрывает.

 

Продолжение следует.

 

Алексей Шуленин