Этот пример показывает, как оценить аналемму Солнца. Аналемма является кривой, которая представляет изменение углового смещения Солнца от его среднего положения на небесной сфере относительно определенной геолокации на поверхности Земли. В этом примере аналемма оценивается относительно Королевской обсерватории в Гринвиче, Великобритания. После оценки пример строит график аналеммы.
В этом примере используются данные, которые можно загрузить с помощью команды 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, deltaUT1, с помощью измененных юлианских дат для UTC.
mjdUTC = mjuliandate(dvUTC); dUT1 = deltaUT1(mjdUTC);
Вычислите полярное движение и перемещение 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');
