УДК 681.3:527(075.8)
ГРАФИЧЕСКОЕ МОДЕЛИРОВАНИЕ АСТРОНОМИЧЕСКОЙ ОБСЕРВАЦИИ
А.А. Барабаш, Дальрыбвтуз, Владивосток
В статье приводится описание методики реализации графоаналитического метода астрономической обсервации по данным наблюдений нескольких светил в электронной таблице MS Excel. Аналитическая обработка исходных данных ведется в модулях VBA Project, а графические построения и результаты вычислений отображаются на листе электронной таблицы. Статья полезна для подготовки судоводителей по программе дисциплины «Мореходная астрономия».
В основу графоаналитического метода обработки астронавигационных наблюдений положен метод высотных линий положения (ВЛП), когда элементы линии - переносы n и азимуты A -вычисляются аналитически, а прокладка линий выполняется графическим путём на бланке или навигационной карте. Для получения обсервованного места требуются данные наблюдений от двух до пяти светил, на основании которых получают линейные уравнения ВЛП:
где А<р - приращение широты относительно счислимого места; Aw -отшествие или Aw/cos <р = АЛ - приращение долготы; А - азимут светила; п - перенос, n = h0 - hc - разность между обсервованной и счислимой высотами светила.
Поскольку в любых измерениях присутствуют случайные погрешности, то для их уравнивания при решении системы уравнений (1) применяют метод наименьших квадратов (МНК): его графический или аналитический варианты. При аналитическом решении систему уравнений из нескольких ВЛП приводят к нормальному виду из двух линейных уравнений (2):
(1)
A#?£cos2/\ + Aw]rcos/\sin/\ = £ncosA; A<pEcos/\sin/\ + Aw^sin2 A = XnsinA
и решая эту систему обычным путём, получают вероятнейшие значения приращений координат А <р и А и/ (3)
:
^соб2А^]пбіпА -^біп АсобА^псобА 2]соб2 А^]біп2 А-^(біпАсобА)2 ’ соб2 А^псобА - ^біп АсобА^пбіп А соб2 А^біп2 А - Л(біп А соб А)2
а далее сами искомые координаты ср0 и Ло (4):
<ро=<рс+а<р:
Я = Л + ДИ//СОБ(з|
(4)
Аналитический вариант МНК в виде системы формул (3) удобен для машинного решения, но страдает полным отсутствием наглядности, столь необходимой в учебной практике. Более приемлемым вариантом в этом случае может служить центрографический приём отыскания координат вероятнейшего места, когда уравниваются не сами линии положения, а координаты точек попарного пересечения линий в фигуре погрешностей. Реализация центрографического приёма выполняется в три этапа:
1. Каждой точке пересечения двух линий положения приписывается вес Р/.
Ріі = к вігі2 (А/ - А ) = к вігі2 (ДА,у),
(5)
где к = Р)Р] - произведение весов линий положения; при равноточных наблюдениях можно принять к = 1.
2. Вычисляются координаты всех точек пересечения линий положения:
п вігі А -п вігі А
А(р..=-----------І------І-------
и
віп ДА..
ч
п сов А -п сов А
Ш = —!■--------------------------------1-1---------і.
Ч віп ДА
іі
(6)
количество точек N пересечения зависит от количества линий положения п и определяется из выражения: N = п (п -1) /2.
3. Рассчитываются вероятнейшие значения приращения координат
А(рв и Ли/е:
(7)
А<Ре = Т.Рц^<Рц /ЦРу',
Ам/е
Вероятнейшие значения приращений координат добавляются к счислимым координатам так же, как и в формуле (4), и получают искомые обсервованные координаты. Оценка точности обсервации в
этом случае производится удвоенным значением радиальной погрешности, величину которой можно рассчитать по формуле
м = +2 ■yjl.PilTPr (8)
Приведённая здесь методика получения координат обсервованного места по формулам (5-8) использована нами для моделирования графической прокладки линий положения на рабочем листе электронной таблицы MS Excel. Аналитические расчёты при этом выполняются в модуле VBA Project, а конечные результаты обработки данных и «рисунок» обсервации отображаются на рабочем листе MS Excel. Программа вычислений в модуле содержит три процедуры: Sub Рисунок(); Sub Место_МНК() и Sub Рис_Чистить(). Для обеспечения возможности обмена данными между процедурами вначале объявляются массивы переменных типа Global.
Процедура Sub Рисунок () обеспечивает управление всем процессом моделирования обсервации по командам пользователя, реализуемых с помощью двух элементов управления - экранные кнопки «Рисунок» и «Стереть». Активизация первой кнопки обеспечивает вывод на экран графической модели обсервации, один из примеров которой показан на рисунках 1, 2.
А В С D
1 Азимуты 231,7° 179,7° 151*
2 Переносы 1,1' 1,3’ 2,7'
3
4 Д <р - 2’ kS
5 ДЛ = 0,6' kE
6 <о0- 41° 59,1' N
7 \ ~ 132° 6' Е
8 2Мо = ±1,8'
9
10 нрнав ■
11
12 ■
13 Стереть ■
Рис. 1. Данные для построения высотных линий положения (азимуты, переносы), координаты места (<ра Л0) и оценка точности обсервации (удвоенное значение радиальной погрешности - 2М0)
Мс ВЛП 3
г' \\
* л
A1 N 1 г k 1 \ \ 1
\ t \
% \ ^ ч ВЛП 2
Мо\(®) \ /
уЛ
\ АЗ
ВЛП 1
Рис. 2. Пример прокладки линий положений на листе электронной таблицы: Мс - счислимое место; Мо - обсервованное место; А1 - азимуты;
ВЛП1,2,3 - высотные линии положения
Итоги обсервации, показанные на рис. 1, формируются в соответствии с результатами прокладки ВЛП (см. рис. 2) и их обработки по формулам (5-7). Для упрощения обработки чисел, заданных в градусной и часовой мере, их ввод и вывод производится в несколько ячеек электронной таблицы, например, значение широты места ср0 = 41 °59,1 'N записывается в ячейки вб - 41°, сб - 59,1' и D6 - N. Для вычислительной обработки подобная запись преобразовывается в десятичное число со знаком. Знак определяется по наименованию (N -"+" и S - "-") и присваивается числу - в нашем примере: 41 + 59,1/60 = 41,985.
В процедуре Sub Рисунок() после объявления внутренних переменных и массивов встроен блок защиты от повторного наложения рисунков друг на друга, когда пользователь, не удалив прежние построения, повторно щёлкает по кнопке «Рисунок». В этом случае производится проверка контрольной ячейки, а с помощью функции MgBox () выводится сообщение, предлагающее ему удалить прежний рисунок, щёлкнув по кнопке «Стереть». Программный код, фрагмент которого показан ниже, включает управляющую структуру, If условие
Then переход, обеспечивающую проверку содержимого ячейки F1 и передачу управления следующей строке кода или переход на метку metka0 в зависимости от содержания ячейки.
Range("F1'').Select
If ActiveCell.FormulaR1C1 <> "" Then GoTo metka0
metka0:
MsgBox ("Сотрите прежний рисунок, щёлкнув по кнопке Стереть"), , "Предупреждение"
Exit Sub
Проверка правильности заданных пользователем исходных параметров обсервации ведётся на начальном этапе. Неверно заданные параметры лишь комментируются программой и не сопровождаются принудительным запретом на дальнейшую работу.
За начало графических построений выбрано постоянное место в координатах xO, уО, численно равных центру ячейки H12 (приблизительно середина экрана при его масштабе 100 %). На рис. 2 это место обозначено символом Мс. Для построения линий использован метод AddLine с установкой набора необходимых свойств, например:
ActiveSheet.Lines.Add(xo, yo - 5, xo, yo + 5).Select With Selection.Border .LineStyle = xlContinuous .ColorIndex = 5 .Weight = xlMedium End With
Для вывода форматированных текстов в заданные ячейки таблицы использовался метод Range, например, в следущем примере в середину ячейки H11 выводится текст «Mc» синего цвета:
Range("H11").Select Selection.Font.ColorIndex = 5 ActiveCell.FormulaR1C1 = " Mc"
Поскольку тип шрифта и его размер явно не заданы, то параметры шрифта устанавливаются по умолчанию: «Arial Cyr, 10».
Масштаб рисунка. При прокладке ВЛП на бумаге пользователь выбирает масштаб для графических построений, руководствуясь двумя соображениями: рисунок должен полностью разместиться на листе и быть как можно крупнее. То же реализовано и в программном коде: с помощью цикла For - Next вычисляется среднеарифметическое значение из всех переносов и присваивается переменной Crpr; а с помощью цикла Do While
- Loop подбирается рациональный коэффициент К для вывода графики. For i = 1 To N Spr = Spr + Abs(npr(i))
Next
Crpr = (Spr / N) * K Do While Crpr > 100 K = K - 2
Crpr = (Spr / N) * K Loop
Линии азимутов1 (на рис. 2 - пунктирные линии-стрелки, обозначенные символами Ai). Экранные координаты линий азимутов (начало - xo + dx, yo - dy; конец - x(i) + dx, y(i) - dy) рассчитываются в структуре цикла For - Next.
For i = 1 To N Az = Ac(i)
If npr(i) < 0 Then Az = Az + 180 If Az > 360 Then Az = Az - 360 yg = Az * pr
dx = R * Sin(yg) * 0.93: dy = R * Cos(yg) x(i) = xo + K * Abs(npr(i)) * Sin(yg) * 0.93 y(i) = yo - K * Abs(npr(i)) * Cos(yg)
' стрелки азимутов
ActiveSheet.Lines.Add(xo + dx, yo - dy, x(i) + dx, y(i) - dy).Select
Next
Направление рисования определяется проверкой знака переноса npr(i) в структурах принятия решений If Then: при отрицательном переносе стрелка азимута выводится в противоположную сторону Az = Az + 180. Координаты конечных точек линий азимутов x(i) и y(i) вычисляются по формулам:
Xj = x0+K- Abc (nprj ) • sin( yg )■ 0,93 y,=y0+K- Abc(npr,■) ■ cos( yg)
с учётом экранных искажений (коэффициент 0,93).
Для вывода линий на экран использован метод. Lines.Add (начало, конец). Select, позволяющий форматирование графического объекта в структуре With Selection. ... End
With.
Название азимутов. Стрелки азимутов на рис. 2 обозначены символами Ai, вывод которых происходит в адресные ячейки электронной таблицы. Вычисление координат ячеек (Row - строка, Col
- колонка) выполняется пересчётом экранных координат в координаты ячеек с накоплением их названий в массиве NAz(i^ последующим форматированием размера Font.Size и цвета Font.ColorIndex индексов:
' имена азимутов
Row = Int((y(i) - dy) / 13): Col = Int((x(i) + dx) / 42)
1 Для рисования линий может использоваться метод AddLine: Shapes.AddLine(10, 10, 250, 250).Line.
NAr(i) = Row: NAc(i) = Col Cells(Row, Col).Value = NAz(i)
Cells(Row, Col).Select Selection.Font.Size = 10 Selection.Font.ColorIndex = 14
Линии положения (на рис. 2 - ВЛП1,2,3). Координаты линий положения в отличие от линий азимутов вычисляются относительно координат определяющей точки на линии азимута (x(i), y(i)) в два этапа: вначале вычисляются координаты начала линии x1, y1 и кооректируются, если их значения выходят за границы экрана; затем вычисляются координаты конца линии x2, y2 и вновь корректируются, если их значения выходят за границы зоны рисунка. Для проверки этих условий использованы блочные структуры If Then End If и операторы принудительного перехода GoTo Метка:
'линии положения L1 = l: L2 = l 1:
dlx = L1 * Cos(yg) * 0.93: dly = L1 * Sin(yg) x1 = x(i) - dlx: y1 = y(i) - dly
If x1 < 1 Or y1 < 1 Then
L1 = L1 - 5 GoTo 1 End If 2:
dlx = L2 * Cos(yg) * 0.93: dly = L2 * Sin(yg)
x2 = x(i) + dlx: y2 = y(i) + dly
If x2 < 1 Or y2 < 1 Then
L2 = L2 - 5
GoTo 2
End If
Сама линия положения на основании вычисленных координат начала - x1, y1 и конца - x2, y2 выводится на экран с форматированием прежним методом:
ActiveSheet.Lines.Add(x1, y1, x2, y2).Select With Selection.Border .LineStyle = xlContinuous .ColorIndex = 11 .Weight = xlThin End With
Названия линий положения (на рис. 2 - ВЛП1,2,3). Адреса ячеек Cells(Row, Col), в которые выводятся названия линий положения, вычисляются относительно начальных координат линий (х1, у1) и затем корректируются, если выходят за границы заданной области рисования:
Row = Int(y1 * 0.088): Col = Int(x1 * 0.02)
If Row < 1 Then Row = 1 'корректировка адресов
If Col < 1 Then Col = 1 'по строкам и колонкам
If Row > 20 Then Row = 20
If Col > 17 Then Col = 17
Ncr(i) = Row: Ncc(i) = Col
Cells(Row, Col).Value = NL(i) 'имена линий положения Cells(Row, Col).Select 'форматирование текста
Selection.Font.Size = 10 Selection.Font.ColorIndex = 11
Отыскание обсервованного места (Мо на рис. 2). Координаты обсервованного места вычисляются в процедуре Sub место_мнк() по методу наименьших квадратов (5-8).
В этой процедуре вначале резервируются локальные массивы (двухмерные и одномерные) координат вершин фигуры погрешностей и их весов:
Dim df(1 To 10), dw(1 To 10) As Single Dim P(1 To 10)
Dim Sp, Sdf, Sdw, srA, Ar(1 To 10), а затем выполняется заполнение этих массивов с помощью процедуры из двух вложенных циклов:
For i = 1 To N - 1 For j = i + 1 To N srA = Sin(Ar(j) - Ar(i))
P(Ct)= srA Л 2 'разность широт вершины
df(Ct)=(npr(i)*Sin(Ar(j))-npr(j)*Sin(Ar(i)))/srA ' отшествие вершины
dw(Ct)=(npr(j)*Cos(Ar(i))-npr(i)*Cos(Ar(j)))/srA
Ct = Ct + 1 ' счётчик
Next
Next
Средневзвешенные координаты (Atp = Sdf и Aw = Sdw) искомого места вычисляются для координат на равносторонней плоскости в цикле из N = n (n -1) /2 повторений:
For i = 1 To Int(N * (N - 1) / 2)
Sp = Sp + P(i)
Sdf = Sdf + P(i) * df(i) 'разность широт места
Sdw = Sdw + P(i) * dw(i) 'отшествие места
Next
Полученные в процедуре Sub место_мнк() средневзвешенные приращения координат возвращаются в процедуру Sub Рисунок (), на их основе рассчитываются экранные координаты обсервованного места:
dx = wd * K * 0.93: dy = -fd * K
CenterX = xo + dx: CenterY = yo + dy 'координаты
графического места
Относительно точки в координатах CenterX, CenterY на экране рисуется символ астрономической обсервации из трёх вложенных
окружностей разных радиусов. Для этого использован метод
Ovals.Add(X, Y, Rx, Ry) .
ActiveSheet.Ovals.Add(CenterX-7,CenterY-
7,15,4).Select
ActiveSheet.Ovals.Add(CenterX-4,CenterY-4,9,8).Select ActiveSheet.Ovals.Add(CenterX-1,CenterY-1,3,2).Select Символ обсервованного места Мо заносится в ячейку электронной таблицы, параметры которой вычисляются пересчётом координат точек в номера строки Ncr(N + 1) и столбца Ncc(N + 1):
Ncr(N + 1) = Int(CenterY / 14) + 1 Ncc(N + 1) = Int(CenterX / 45)
Cells(Ncr(N + 1), Ncc(N + 1)).Value = " Mo" Результаты обсервации. Числовые параметры обсервованного места (приращение координат Atp и ДЯ и сами координаты <р и Я), показанные на рис. 1, формируются на основании анализа наименований (N или S; E или W) самих координат и их приращений. При совпадении их наименований выполняется сложение, а при несовпадении - вычитание:
Namef = " N" ' наименование
координат
If FoGr < 0 Then Namef = " S"
Namel = " E"
If LoGr < 0 Then Namel = " W"
Namedf = "k N" ' наименование
приращений
If fd < 0 Then Namedf = "k S"
Namedl = "k E"
If ld < 0 Then Namedl = "k W"
Затем полученные символы заносятся в фиксированные ячейки электронной таблицы.
Стирание рисунка. Выполняется с помощью процедуры Sub Рис_Чистить(), которая запускается экранной кнопкой "Стереть". Очистка элементов графических построений организована в циклах по группам однородных объектов. Например, для удаления линий положения (графических объектов) использован цикл, содержащий строку поиска линии, ее выделения и стирания - метод Delete:
For i = 1 To N * 2 + 2
ActiveSheet.Lines(1).Select
Selection.Delete
Next
Необходимость удвоенного значения повторений в цикле (N * 2 + 2) вызвана тем обстоятельством, что каждая линия положения
2 В VBA 6.3 и выше использован более унифицированный метод Shapes.AddShape (msoShapeOval, x, y, Rx, Ry).Select.
рассчитывалась и строилась в два приёма: влево и вправо
относительно определяющей точки на линии азимута.
Для удаления текстовых объектов, например, названий азимутов, использован простой цикл, приводящий содержимое ячеек в состояние
«пусто»:
For i = 1 To N
Cells(NAr(i), NAc(i)).Value = ""
Next
Для удаления названий линий положения использована система из двух вложенных циклов, также приводящая содержимое заполненных ячеек в состояние «пусто»:
For i = 2 To 7 For j = 2 To 3 Cells(i, j).Value = ""
Next
Next
Библиографический список
1. Красавцев Б.И. Мореходная астрономия: Учеб. для вузов. М.: Транспорт, 1986. С. 255.
2. Немцев О.В. Приборы и методы астронавигации: Учеб. пособие. Владивосток: Дальрыбвтуз, 2001. С. 165.