Положение Солнца на небе

[xd]

Интересно, писатели в этой теме собираются что-нибудь читать? Хотя бы эту же тему с самого начала? Возьмите любую книгу, касающуюся хотя бы поверхностно эфемеридной астрономии. Там всё это будет.

[kikiriki0]

склоенине Delta и уравнение времени можно посмотреть тут
http://aa.usno.navy.mil/faq/docs/SunApprox.php

[Kurkuma]

не мог бы кто выложить страницы из астрономического ежегодника или кинуть на него ссылку: нужны значения звездного времени в среднюю гринвичскую полночь соответствующей календарной даты за август, сентябрь 2009г

[IOPA4]

Sunrise/Sunset Algorithm Example

Source:
Almanac for Computers, 1990
published by Nautical Almanac Office
United States Naval Observatory
Washington, DC 20392

Inputs:
day, month, year: date of sunrise/sunset
latitude, longitude: location for sunrise/sunset
zenith: Sun's zenith for sunrise/sunset
offical = 90 degrees 50'
civil = 96 degrees
nautical = 102 degrees
astronomical = 108 degrees

NOTE: longitude is positive for East and negative for West

Worked example (from book):
June 25, 1990:   25, 6, 1990
Wayne, NJ: 40.9, -74.3
Office zenith: 90 50' cos(zenith) = -0.01454


1. first calculate the day of the year

N1 = floor(275 * month / 9)
N2 = floor((month + 9) / 12)
N3 = (1 + floor((year - 4 * floor(year / 4) + 2) / 3))
N = N1 - (N2 * N3) + day - 30

Example:
N1 = 183
N2 = 1
N3 = 1 + floor((1990 - 4 * 497 + 2) / 3)
= 1 + floor((1990 - 1988 + 2) / 3)
= 1 + floor((1990 - 1988 + 2) / 3)
= 1 + floor(4 / 3)
= 2
N = 183 - 2 + 25 - 30 = 176

2. convert the longitude to hour value and calculate an approximate time

lngHour = longitude / 15

if rising time is desired:
t = N + ((6 - lngHour) / 24)
if setting time is desired:
t = N + ((18 - lngHour) / 24)

Example:
lngHour = -74.3 / 15 = -4.953
t = 176 + ((6 - -4.953) / 24)
= 176.456

3. calculate the Sun's mean anomaly

M = (0.9856 * t) - 3.289

Example:
M = (0.9856 * 176.456) - 3.289
= 170.626

4. calculate the Sun's true longitude
[Note throughout the arguments of the trig functions
(sin, tan) are in degrees. It will likely be necessary to
convert to radians. eg sin(170.626 deg) =sin(170.626*pi/180
radians)=0.16287]

L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634
NOTE: L potentially needs to be adjusted into the range [0,360) by adding/subtracting 360

Example:
L = 170.626 + (1.916 * sin(170.626)) + (0.020 * sin(2 * 170.626)) + 282.634
= 170.626 + (1.916 * 0.16287) + (0.020 * -0.32141) + 282.634
= 170.626 + 0.31206 + -0.0064282 + 282.634
= 453.566 - 360
= 93.566

5a. calculate the Sun's right ascension

RA = atan(0.91764 * tan(L))
NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/subtracting 360

Example:
RA = atan(0.91764 * -16.046)
= atan(0.91764 * -16.046)
= atan(-14.722)
= -86.11412

5b. right ascension value needs to be in the same quadrant as L

Lquadrant = (floor( L/90)) * 90
RAquadrant = (floor(RA/90)) * 90
RA = RA + (Lquadrant - RAquadrant)

Example:
Lquadrant = (floor(93.566/90)) * 90
= 90
RAquadrant = (floor(-86.11412/90)) * 90
= -90
RA = -86.11412 + (90 - -90)
= -86.11412 + 180
= 93.886

5c. right ascension value needs to be converted into hours

RA = RA / 15

Example:
RA = 93.886 / 15
= 6.259

6. calculate the Sun's declination

sinDec = 0.39782 * sin(L)
cosDec = cos(asin(sinDec))

Example:
sinDec = 0.39782 * sin(93.566)
= 0.39782 * 0.99806
= 0.39705
cosDec = cos(asin(0.39705))
= cos(asin(0.39705))
= cos(23.394)
= 0.91780

7a. calculate the Sun's local hour angle

cosH = (cos(zenith) - (sinDec * sin(latitude))) / (cosDec * cos(latitude))

if (cosH > 1)
the sun never rises on this location (on the specified date)
if (cosH < -1)
the sun never sets on this location (on the specified date)

Example:
cosH = (-0.01454 - (0.39705 * sin(40.9))) / (0.91780 * cos(40.9))
= (-0.01454 - (0.39705 * 0.65474)) / (0.91780 * 0.75585)
= (-0.01454 - 0.25996) / 0.69372
= -0.2745 / 0.69372
= -0.39570

7b. finish calculating H and convert into hours

if if rising time is desired:
H = 360 - acos(cosH)
if setting time is desired:
H = acos(cosH)

H = H / 15

Example:
H = 360 - acos(-0.39570)
= 360 - 113.310 [ note result of acos converted to degrees]
= 246.690
H = 246.690 / 15
= 16.446

8. calculate local mean time of rising/setting

T = H + RA - (0.06571 * t) - 6.622

Example:
T = 16.446 + 6.259 - (0.06571 * 176.456) - 6.622
= 16.446 + 6.259 - 11.595 - 6.622
= 4.488

9. adjust back to UTC

UT = T - lngHour
NOTE: UT potentially needs to be adjusted into the range [0,24) by adding/subtracting 24

Example:
UT = 4.488 - -4.953
= 9.441
= 9h 26m

10. convert UT value to local time zone of latitude/longitude

localT = UT + localOffset

Example:
localT = 9h 26m + -4
= 5h 26m
= 5:26 am EDT

Ваше мнение

[kuksha]

Попробовал забить алгоритм Tau, который дан на первой странице темы.
Не сходится у меня арифметика. Кроме того алгоритм подправить бы - наверняка ещё кому-то понадобится...

Сначала - что имеем:
Извините за цитирование, но удобнее привести цитаты до третьего пункта, до которого я дошёл:

Дано: координаты наблюдателя – долгота Lon, геодезическая широта Lat.
Момент наблюдения – год Year, месяц Mon, день Day, всемирное время в часах UT.

Найти: Азимут Az, высоту над горизонтом H

1.Вычисление модифицированной  юлианской даты на начало суток

 Море готовых процедур на всех языках.
Например,
 Var1 = 10000 * Year + 100 * Mon + Day
 Если Mon <= 2    Mon = Mon + 12, Year = Year - 1
 Если Var1 <= 15821004, то  Var2 = -2 + (Year + 4716) \ 4 - 1179
 Иначе  Var2 = Year \ 400 - Year \ 100 + Year \ 4
 Var3 = 365 * Year - 679004
MD = Var3 + Var2 + 306001 * (Mon + 1) \ 10000 + Day
\ - деление нацело

2. Вычисление местного звездного времени

T0 = (MD - 51544.5) / 36525 ‘ мод.юл.дата на начало суток в юлианских столетиях
a1 = 24110.54841
a2 = 8640184.812
a3 = 0.093104
a4 = 0.0000062
S0 = a1 + a2 * T0 + a3 * T0 ^ 2 - a4 * T0 ^ 3 ' звездное время в Гринвиче на начало суток в секундах
Nsec =UT * 3600 ‘ количество секунд, прошедших  от начала суток до момента наблюдения
UT - всемирное время в часах, момент расчета
NsecS = Nsec * 366.2422 / 365.2422 ‘ количество звездных секунд
SG = (S0 + NsecS) /3600 * 15 ‘ гринвическое среднее звездное время в градусах
ST = SG + Lon ‘ местное звездное время
Lon – долгота наблюдателя

3. Вычисление эклиптических координат Солнца

T0 = (MD - 51544.5) / 36525 ‘ мод.юл.дата на начало суток в юлианских столетиях
UT – время в часах от полуночи даты
M = 357.528 + 35999.05 * T0 + 0.04107 * UT ‘ средняя аномалия
L0 = 280.46 + 36000.772 * T0 + 0.04107 * UT
L = L0 + (1.915 - 0.0048 * T0) * Sin(M) + 0.02 * Sin(2 *M)  ‘ долгота Солнца

X = Cos(L)  ' вектор
Y = Sin(L)   '  в эклиптической
Z = 0         '    системе координат

Ну и ещё немного - пример вычислений с той же страницы:

Lon = 41.43147  долгота
Lat = 43.64986  широта
Year = 2009, Mon = 6, Day = 19, UT = 14.0

1.Вычисление модифицированной юлианской даты на начало суток
var1 =   20090619
var2 =   487
var3 =   54281
MD =    55001

2. Вычисление местного звездного времени
Входные параметры: MD = 55001, UT = 14.0
T0 =   9.46338124572211E-02
S0 = 841764.178338331
Nsec =  50400
NsecS =  50537.9906264939
SG =  117.925704020104
ST = 159.357174020104

3. Вычисление эклиптических координат Солнца
Входные параметры: MD = 55001, UT = 14.0
T0=    9.46338124572211E-02
M =   164.830326338125
L0=   3687.92528576318
L =   88.4161785572869 ' приведение угла в диапазон 0-360 градусов
X =    2.76393789023711E-02
Y =    0.999617959389432
Z =    0

У меня сошлось вплоть до SG и М:
SG=3717.9257040201
M=3764.83032633812
Вычел 3600 и сошлось. Ясно, что цикличность, но в алгоритме этого всё-же нет...
Как тут правильно приведение сделать? Или где-то выше размерность поправить?

Но, с вычислением L всё хуже - получаю
L =   89.833730218558
вместо
L =   88.4161785572869
Вот этого я уже не понимаю...

Помогите пожалуйста.

[Sergey007]

Но, с вычислением L всё хуже - получаю
L =   89.833730218558
вместо
L =   88.4161785572869
Вот этого я уже не понимаю...

Видимо какие то округления при вычислениях, подсчитал калькулятором получилось 88.4161785572869 а в программе 89.833730218558.
Никак не могу понять эклиптические координаты Солнца. ( Подсчитал по этому примеру для своей долготы результат такой
X =  0.934708654090602
Y =  0.355414873026067
L = 123.283242893464
Как определить по этим данным где сейчас солнце? )
Что от чего измерять, если не затруднит объясните по простому. На пальцах

[xd]

По-простому читаем Астрономический календарь, постоянную часть. Есть в интернетах

[Sergey007]

Я не силен в астрономии (  Если не сложно объясните?

[Upsilon]

Вычел 3600 и сошлось. Ясно, что цикличность, но в алгоритме этого всё-же нет...

В конечном итоге все угловые величины являются аргументами тригонометрических функций, а в программировании тригонометрическим функциям нет разницы - сколько дополнительных 2*пи прибавлено к аргументу.

вместо
L =   88.4161785572869
Вот этого я уже не понимаю...

В программировании аргумент тригонометрической функции должен быть в радианах,у вас при вычислении L средняя аномалия М является аргументом синуса и вы не перевели ее в радианы, вот и получили 89.83373

Как определить по этим данным где сейчас солнце? )

Перевести  эклиптические координаты в экваториальные, а экваториальные - в азимутальные. Deimos дал ссылку. Формулы настолько распространенные, что перепечатывать их тут все равно что букварь цитировать.

[Sergey007]

Спасибо это исправил нужно аргумент Sin (x * pi/180)
Сейчас получился такой результат
X=-0.523441862791698
Y=0.852061392316984
Z=0
L=121.563409331487

Из примера с первой страницы

Код: [Выделить]

Eps = 23.439281   ‘ наклон эклиптики к экватору
 X’ =  X                                    '  вектор
 Y’ = Y * Cos(Eps) - Z * Sin(Eps) '   в экваториальной
 Z’ = Y * Sin(Eps) + Z * Cos(Eps) '    системе координат

5. Экваториальные геоцентрические координаты Солнца

tg (Ra) = Y’ /X’
tg (Dec)= Z’ /  Sqrt(X’ ^ 2 + Y’ ^ 2)
Ra  и Dec  это какие то константы как их получить?
Что  означает операция апостроф (верхняя запятая)?
Как эту операцию заменить на с++?

 

[xd]

Читайте книги, а?
Там написано лучше чем мы по-быстрому на пальцах расскажем.

[Sergey007]

Книги читаем. Но когда профи объясняют это понятнее.

[Sergey007]

L=121.563409331487 в эклиптических координатах это угол на восток от точки весеннего равноденствия, правильно ?

[xd]

Против часовой стрелки, если смотреть с севера.

[xd]

Книги тоже профессионалы писали. Которые ещё и излагать умеют.

[Sergey007]

Да согласен, без книг никуда. Но я ничего не понимаю, основу все равно нужно объяснить. 
L=121.563409331487 в эклиптических координатах это угол на восток от точки весеннего равноденствия, правильно я понял ?

[xd]

Да нету востока на эклиптике!

[Sergey007]

На вике прочел http://ru.wikipedia.org/wiki/%D0%AD%D0%BA%D0%BB%D0%B8%D0%BF%D1%82%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%B0%D1%8F_%D1%81%D0%B8%D1%81%D1%82%D0%B5%D0%BC%D0%B0_%D0%BA%D0%BE%D0%BE%D1%80%D0%B4%D0%B8%D0%BD%D0%B0%D1%82
Как представить эту систему что бы понять от какой точки считать?
А обязательно переводить из эклиптические координаты в экваториальные, а экваториальные - в азимутальные что бы определить где солнце относительно наблюдателя?

[Upsilon]

Что  означает операция апостроф (верхняя запятая)?

По-моему, это не операция, это переменные так обозначаются. Они ранее вычислялись.

Как представить эту систему что бы понять от какой точки считать?
А обязательно переводить из эклиптические координаты в экваториальные, а экваториальные - в азимутальные что бы определить где солнце относительно наблюдателя?

Если у вас есть пространственное воображение, то смотрите сюда:
Представьте себе 3d-пространство, в нем есть звезда по имени Солнце. Проведите через центр Солнца плоскость и назовите ее плоскостью эклиптики. В этой плоскости вообразите окружность с центром в Солнце. По этой окружности Земля бегает вокруг Солнца, делая один оборот за год. Нарисуйте в плоскости вектор от Солнца в любое направление, пересечение этого вектора с окружностью даст точку - начало отсчета измерения эклиптической долготы. Назовите эту точку точкой весеннего равноденствия. Это эклиптическая гелиоцентрическая (с центром в Солнце) система координат. В ней положение небесного тела (в данном случае Земли) характеризуется двумя величинами  - эклиптической долготой и широтой. Долгота - это угол между вектором "Солнце-точка весеннего равноденствия" и вектором "Солнце-Земля (или другое тело)". Широта - это угол между плоскостью эклиптики и телом. Орбита Земли почти лежит в плоскости эклиптики, так что считаем эклиптическую широту нулевой.

Нам нужны не координаты Земли относительно Солнца, а координаты Солнца относительно Земли, поэтому переходим к эклиптической геоцентрической (с  центром в Земле) системе координат. Для этого вычисляем координаты не вектора "Солнце-Земля", а противоположно направленного вектора "Земля-Солнце". Это вы уже научились как-то делать и получили эклиптическую долготу Солнца.  А широту считаем нулевой.

Если бы ось Земли была перпендикулярна плоскости эклиптики, то наша наблюдаемая с Земли небесная система координат совпадала бы с эклиптической системой координат. Но ось Земли наклонена и поэтому от эклиптических координат нужно перейти к экваториальным, к тем координатам, которые связаны с ориентацией Земли в пространстве. Экваториальных систем координат две - первая и вторая. Первая связана с невращающейся вокруг оси Землей, вторая - с вращающейся.  Мы переходим ко второй экваториальной системе координат, в которой положение тела характеризуется двумя углами  -  прямым восхождением RA и склонением небесного тела DEC.

После этого нужно перейти к первой экваториальной системе координат, "остановить Землю", вычислить часовой угол тела (склонение то же, что и во второй экваториальной системе координат). А после этого учесть различие между экваториальной и азимутальной системами координат, потому что эти системы координат ориентированы в пространстве по-разному. Только тогда вычисляются азимут тела и высота над горизонтом.

Короче, да, переводить координаты надо обязательно. Это по-разному ориентированные в пространстве системы координат. Координаты Солнца вычисляют легко, именно потому что они вычисляются в эклиптической системе координат. Но у нас есть еще Земля, которая наклонена к плоскости эклиптики и крутится вокруг своей оси. И еще есть наблюдатель в какой-то точке Земли, а местная азимутальная система координат для каждой точки Земли своя.

[Sergey007]

Большое спасибо за подробное объяснение.
Для геоцентрической гелиоцентрической есть константы начала отсчета цифровые ?
По примеру с первой страницы вроде все совпадает, кроме вычисления Ra Dec
В примере так:

tg (Ra) = Y’ /X’
tg (Dec)= Z’ /  Sqrt(X’ ^ 2 + Y’ ^ 2)

XY совпадают они правильные. вычисляю Ra Dec так

Код: [Выделить]

Ra = Math.Tan(Y1 / X1);   результат -5,0539826731987 а должно быть 88.2738131051079
Dec = Math.Tan( Z1 / Math.Sqrt(Math.Pow(X1, 2) + Math.Pow(Y1, 2)) );  результат 0,462689303803096 а должно быть  23.4297911659417
В чем может быть дело?