Оценка аналеммы Солнца с использованием планетарных эфемеридов и ECI в преобразование AER

Открыть скрипт

Этот пример показывает, как оценить аналемму Солнца. Аналемма является кривой, которая представляет изменение углового смещения Солнца от его среднего положения на небесной сфере относительно определенной геолокации на поверхности Земли. В этом примере аналемма оценивается относительно Королевской обсерватории в Гринвиче, Великобритания. После оценки пример строит график аналеммы.

В этом примере используются данные, которые можно загрузить с помощью команды aeroDataPackage.

Определите даты года, за которые вычислить аналемму Солнца

Укажите даты, для которых нужно вычислить аналемму. В этом примере эти даты варьируются для 2014 года с 1 января по 31 декабря в полдень UTC.

dv = datetime(2014,1,1:365,12,0,0);
dvUTC = [dv.Year' dv.Month' dv.Day' dv.Hour' dv.Minute' dv.Second'];

Вычислим положение Солнца

Используйте функцию planetEphemeris, чтобы вычислить положение Солнца. В этом примере:

Функция tdbjuliandate вычисляет юлианскую дату для динамического барицентрического времени (TDB).
Функция tdbjuliandate требует наземного времени (TT).

Для вычисления наземного времени в секундах от UTC требуется различие в скоординированном универсальном времени (UTC) и международном атомном времени (TAI).

На 2014 год это различие (dAT) составляет 35 секунд.
Приблизительное наземное время (secTT) составляет dAT + 32,184 секунд.
Наземное время в году, месяце, дне, часе, минутах и секундах содержится в массиве dvTT.

dAT = 35;
secTT = dAT + 32.184;
dvTT = dv + secTT/86400;

Оцените юлианскую дату для динамического барицентрического времени на основе наземного времени с помощью массива dvTT:

jdTDB = tdbjuliandate([dvTT.Year' dvTT.Month' dvTT.Day' dvTT.Hour' dvTT.Minute' dvTT.Second']);

Определите положение Солнца для этих дат:

posSun = planetEphemeris(jdTDB,'Earth','Sun')*1000;

Вычислите различие между UTC и основным универсальным временем (UT1)

Вычислите различие между UTC и UT1, deltaUT1, с помощью измененных юлианских дат для UTC.

mjdUTC = mjuliandate(dvUTC);
dUT1 = deltaUT1(mjdUTC);

Вычисление полярного движения и смещения небесного промежуточного полюса (CIP)

Вычислите полярное движение и перемещение CIP с помощью измененных юлианских дат для UTC.

PM = polarMotion(mjdUTC);
dCIP = deltaCIP(mjdUTC);

Укажите геопотенциальное положение Королевской обсерватории в Гринвиче, Великобритания

Укажите геопотенциальное положение положения, по которому можно оценить аналемму. В этом примере это место является широтой, долготой и высотой для Королевской обсерватории в Гринвиче (51,48 степени севера, 0,0015 степени запада, 0 метров высоты).

LLAGreenwich = [51.48,-0.0015,0];
aer = eci2aer(posSun,dvUTC,repmat(LLAGreenwich,length(jdTDB),1),...
              'deltaAT',dAT*ones(length(jdTDB),1),'deltaUT1',dUT1,...
              'PolarMotion',PM,'dCIP',dCIP);

Задайте дни в течение года аналеммы, которую вы хотите построить

На аналемме можно построить интересующие дни в течение года аналеммы. Этот пример строит графики:

Первый день каждого месяца в 2014 году.
Летнее и зимнее солнцестояния.
Весеннее и осеннее равноденствия.

Чтобы получить первый день каждого месяца в 2014 году:

aerFirstMonth = aer(dvUTC(:,3)==1,:);

Получить солнцестояния и равноденствия (на 2014 год - 3/20, 6/21, 9/22, 12/21):

solsticeEquinox = [ aer(and(dvUTC(:,2)==3,dvUTC(:,3)==20),1) aer(and(dvUTC(:,2)==3,dvUTC(:,3)==20),2); ...
                    aer(and(dvUTC(:,2)==6,dvUTC(:,3)==21),1) aer(and(dvUTC(:,2)==6,dvUTC(:,3)==21),2); ...
                    aer(and(dvUTC(:,2)==9,dvUTC(:,3)==22),1) aer(and(dvUTC(:,2)==9,dvUTC(:,3)==22),2); ...
                    aer(and(dvUTC(:,2)==12,dvUTC(:,3)==21),1) aer(and(dvUTC(:,2)==12,dvUTC(:,3)==21),2)];

Графическое изображение результатов

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

for ii = 12:-1:1
	firstDays{ii} = [num2str(ii) '/' num2str(1)];
end

f = figure;
plot(aer(:,1),aer(:,2),'.',...
    solsticeEquinox(:,1),solsticeEquinox(:,2),'ks',...
    aerFirstMonth(:,1),aerFirstMonth(:,2),'ko',...
    'MarkerSize',8,'MarkerFaceColor','k');
title('Analemma observed at Greenwich Observatory');
xlabel('Azimuth [deg]');
ylabel('Elevation [deg]');
axis([176,185,10,70])
text(aerFirstMonth(:,1)+.1, aerFirstMonth(:,2)+1.2, firstDays, 'Color', 'k','HorizontalAlignment', 'Left');
text(solsticeEquinox(1,1)+.2, solsticeEquinox(1,2)-1.5, 'Spring Equinox', 'Color', 'k','HorizontalAlignment', 'Left');
text(solsticeEquinox(2,1), solsticeEquinox(2,2)+2.5, 'Summer Solstice', 'Color', 'k','HorizontalAlignment', 'Left');
text(solsticeEquinox(3,1)+.1, solsticeEquinox(3,2)+1.2, 'Fall Equinox', 'Color', 'k','HorizontalAlignment', 'Left');
text(solsticeEquinox(4,1)+.1, solsticeEquinox(4,2)-2.5, 'Winter Solstice', 'Color', 'k','HorizontalAlignment', 'Left');

Документация