Тангаж, рысканье в норме. Выводим на орбиту ИСЗ.
Мои поздравления по случаю состоявшейся вчера 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 или ГЛОНАСС, покрутить ее во времени и посмотреть, насколько плотно она кого покрывает.
Продолжение следует.
Алексей Шуленин