A A A A Автор Тема: Расчет положения солнца и луны (Прочитано 18172 раз)

StormBS · « : 08 Окт 2010 [19:29:46] »

Уважаемые форумчани, помогите решить новичку в астрономии задачу. Исходные данные: Дата, время, часовой пояс, координаты точки наблюдения (широта, долгота). Необходимо вычислить азимут, высоту (угол места), прямое восхождение, склонение, долготу, широту, и дистанцию( расстояние от точки наблюдения) для Солнца и Луны. Мне необходимо это для написание программы, которая должна на карте земли (2D) отражать положение данных объектов. На данный момент удалось получить азимут и высоту Солнца, а получить его широту и долготу не могу ни как, что принципиально для программы.

Tau · « **Ответ #1 :** 08 Окт 2010 [19:54:54] »

Цитата

а получить его широту и долготу не могу ни как

Прямое восхождение и склонение Солнца, видимо, известны, раз вы до азимута дошли.
Широта подсолнечной точки равна склонению Солнца.
Долгота подсолнечной точки L:
L = RA - S
RA - прямое восхождение Солнца
S - звездное гринвическое время

AAV · « **Ответ #2 :** 08 Окт 2010 [20:04:23] »

https://astronomy.ru/forum/index.php/topic,21684.msg395573.html#msg395573

StormBS · « **Ответ #3 :** 08 Окт 2010 [20:17:53] »

Цитата: Tau от 08 Окт 2010 [19:54:54]

Цитата
а получить его широту и долготу не могу ни как
Прямое восхождение и склонение Солнца, видимо, известны, раз вы до азимута дошли.
Широта подсолнечной точки равна склонению Солнца.
Долгота подсолнечной точки L:
L = RA - S
RA - прямое восхождение Солнца
S - звездное гринвическое время

На соседней ветки нашел такую формулу:
_________________________________________________ ___________________
RA - прямое восхождение Солнца на нужный момент времени
DEC - склонение Солнца на нужный момент времени
ST - гринвическое звездное время

Долгота подсолнечной точки L:
L = atn2(Y,X) ' перевести в градусы
где
X = cos(ST)*cos(RA)*cos(DEC) + sin(ST)*sin(RA)*cos(DEC)
Y = -sin(ST)*cos(RA)*cos(DEC) + cos(ST)*sin(RA)*cos(DEC)
atn2(Y,X) - функция арктангенса, имеющаяся практически в каждом языке программирования и дающая значение в диапазоне от 0 до 360 градусов, в отличие от обычной функции арктангенса atn(X), которая дает значение в диапазоне от -90 до 90 градусов.

Широта подсолнечной точки F:
F = DEC
_________________________________________________ _____________
Если не сложно объясните в чем разница

Tau · « **Ответ #4 :** 08 Окт 2010 [20:28:13] »

Цитата

Если не сложно объясните в чем разница

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

.
Просто некоторые дамы питают слабость к аналитической геометрии и иногда в состоянии думать только в определенных категориях.

StormBS · « **Ответ #5 :** 08 Окт 2010 [20:37:25] »

Спасибо за общую помощь завтра со свежей головой буду разбираться.

StormBS · « **Ответ #6 :** 09 Окт 2010 [13:47:30] »

Вроде бы получилось. Кому интересен код на VC2010 выкладываю ниже, но без заголовочного файла.

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

#include "StdAfx.h"
#include "GSanSystem.h"
#include "math.h"

#define M_PI       3.14159265358979323846
#define Rad(x)		x*M_PI/180.0
GSanSystem::GSanSystem(void)
{
}

GSanSystem::~GSanSystem(void)
{
}

int GSanSystem::round (double n)
{
	double t;
	t=n-double(n);
	if(t>=0.5)
	{
		n*=10;
		ceil(n);
		n/=10;
	}
	else
	{
		n*=10;
		floor(n);
		n/=10;
	}
	return (int)n;
}


void GSanSystem::SunPos(
			double  lon, //Долгота	
			double  lat,//Широта	
			double  Year,//год
			double  Mon,// месяцы
			double  Day,// дни
			double  hous,//часы
			double  min,//минуты
			double sec,// секунды
			double zona// Часовой пояс
		   )

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

	 double Var1,Var2,Var3;
	 Var1 = 10000 * Year + 100 * Mon + Day;
	 if(Mon <= 2 )
	 {
		Mon = Mon + 12;
		Year = Year - 1;
	 }
	 if( Var1 <= 15821004) Var2 = -2 + floor((Year + 4716) / 4) - 1179;
	 else  Var2 = floor(Year /400) - floor(Year / 100) + floor(Year / 4);
	 Var3 = 365 * Year - 679004;
	 // MJD - Модифицированная Юлианская дата
	 double MJD = Var3 + Var2 + floor(306001 * (Mon + 1)/ 10000) + Day;
	
	
	// Вычисление Гринвеческого звездного времени

	double T0 = (MJD - 51544.5) / 36525; // мод.юл.дата на начало суток в юлианских столетиях
	double a1 = 24110.54841;
	double a2 = 8640184.812;
	double a3 = 0.093104;
	double a4 = 0.0000062;
	double S0 = a1 + a2 * T0 + a3 * T0*T0- a4 * T0 *T0*T0;// звездное время в Гринвиче на начало суток в секундах
	//UT - всемирное время в часах, момент расчета
	double UT = hous-zona + min/60 + sec/3600;
	if(UT>24)UT=UT-24;
	if(UT<0)UT=UT+24;
	double Nsec = UT * 3600; // количество секунд, прошедших  от начала суток до момента наблюдения
	
	double NsecS = Nsec * 366.2422 / 365.2422; //количество  звездных секунд
	double GMT = (S0 + NsecS) /3600 * 15;//гринвическое среднее звездное время в градусах SG
	while(GMT>360)GMT=GMT-360;
	double GST = GMT + lon;// местное звездное время ST
	//Lon – долгота наблюдателя

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

	T0 = (MJD - 51544.5) / 36525; // мод.юл.дата на начало суток в юлианских столетиях
	double M = 357.528 + 35999.05 * T0 + 0.04107 * UT;// средняя аномалия
	while(M>360)M=M-360;
	double L0 = 280.46 + 36000.772 * T0 + 0.04107 * UT;
	double L = L0 + (1.915 - 0.0048 * T0) * sin(Rad(M)) + 0.02 * sin(Rad(2 *M));//долгота Солнца
	while(L>360)L=L-360;

	double X = cos(Rad(L)) ; // вектор
	double Y = sin(Rad(L)) ; //  в эклиптической
	double Z = 0      ; //  системе координат

	

	// Координаты Cолнца в прямоугольной экваториальной системе координат

	double Eps = 23.439281   ; //наклон эклиптики к экватору
	double  X’ =  X ;                          // вектор
	double  Y’ = Y * cos(Rad(Eps)) - Z * sin(Rad(Eps)); //   в экваториальной
	double  Z’ = Y * sin(Rad(Eps)) + Z * cos(Rad(Eps)) ;//    системе координат

	
	// Экваториальные геоцентрические координаты Солнца
	// RA - прямое восхождение Солнца на нужный момент времени
	//DEC - склонение Солнца на нужный момент времени

	double Ra = atan2(Y’ ,X’)*180/M_PI;
	double Dec = atan2 (Z’ , sqrt(X’*X’ + Y’*Y’))*180/M_PI;

	
	//  Азимутальные координаты Солнца

	
	//Lat - широта

	double Th = GST - Ra  ;//часовой угол
	double z  = acos(sin(Rad(lat)) * sin(Rad(Dec)) + cos(Rad(lat)) * cos(Rad(Dec)) * cos(Rad(Th)))*180/M_PI;// косинус зенитного угла
	double H = 90 - z;
	double Az = atan2( sin(Rad(Th)) * cos(Rad(Dec)) * cos(Rad(lat)),sin(Rad(H)) * sin(Rad(lat)) - sin(Rad(Dec)))*180/M_PI;

	// получаем подсолнечную точку
	// Долгота Солнца
	double LonSan = Ra - GST;
	// Широта Солнца
	double LatSan = Dec;

}

Кто подскажет как вычислить тоже самое для Луны, а также Дистанцию до этих объектов.

Tau · « **Ответ #7 :** 09 Окт 2010 [21:13:07] »

Цитата

вычислить тоже самое для Луны

https://astronomy.ru/forum/index.php/topic,78687.0.html
Ответ#15, прицепленный файл.

Цитата

а также Дистанцию до этих объектов.

Если под дистанцией вы понимаете расстояние от центра Земли до Луны/Солнца, то
в случае Луны она равна R*ER (в том doc-е выше найдете это место),
в случае Солнца AE*( 1.000140612 - 0.016708617*cos(M) - 0.000139589*cos(2*M)), где
AE = 149597892.1111 км (величина астрономической единицы)
M - средняя аномалия Солнца которую вы уже вычислили

StormBS · « **Ответ #8 :** 13 Окт 2010 [13:57:21] »

Уважаемая Tau Вы бы не могли привести более полный листинг функции расчета "Вычисление подлунной точки." из прикрепленного файла (Word). Пишу на с++ и данные с Вашими очень разнятся.

Tau · « **Ответ #9 :** 13 Окт 2010 [14:27:32] »

Подозреваю, что они разнятся на этапе вычисления звездного времени, там юлианскую дату округлить надо.
T = (int(MD)-51544.5)/36525
Прошу прощения

StormBS · « **Ответ #10 :** 13 Окт 2010 [14:54:43] »

Очень большая разница в вычисление эклиптической долготы, эклиптической широты и
горизонтального параллакса

T0 = 0.10764145562400176

      // Вычисление эклиптических координат Луны и вектора (X, Y, Z) в эклиптической системе координат.
      // эклиптическая долгота
      double Ekl_Lon = 218.32 + 481267.883*T0 + 6.29*sin(Rad(134.9 + 477198.85*T0)) -
         1.27*sin(Rad(259.2-413335.38*T0)) + 0.66*sin(Rad(235.7 + 890534.23*T0)) +
         0.21*sin(Rad(269.9 + 954397.7*T0)) - 0.19*sin(Rad(357.5 + 35999.05*T0)) -
         0.11*sin(Rad(186.6 + 966404.05*T0));
      // эклиптическая широта
      double Ekl_Lat = 5.13*sin(Rad(93.3 + 483202.03*T0)) + 0.28*sin(Rad(228.2 + 960400.87*T0)) -
         0.28*sin(Rad(318.3 + 6003.18*T0)) - 0.17*sin(Rad(217.6 - 407332.2*T0));
// Горизонтальный паралакс
      double p = 0.9508 + 0.0518*cos(Rad(134.9 + 477198.85*T0)) + 0.0095*cos(Rad(259.2-413335.38*T0)) +
         0.0078*cos(Rad(235.7 + 890534.23*T0)) + 0.0028*cos(Rad(269.9 + 954397.7*T0));

      Ekl_Lon   52028.770816762059
      Ekl_Lat   4.6130842559697580
      p   0.97568441824445629
________________________
#define Rad(x) x*M_PI/180.0

Tau · « **Ответ #11 :** 13 Окт 2010 [15:29:39] »

А у меня эти ваши строчки правильно посчитались

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

#include <iostream.h>
#include <cmath>
void main(void)    
{
  double T0 = 0.10764145562400176;
  double q=3.14159265358979/180.0;
  double Ekl_Lon = 218.32 + 481267.883*T0 + 6.29*sin(q*(134.9 + 477198.85*T0)) -
         1.27*sin(q*(259.2-413335.38*T0)) + 0.66*sin(q*(235.7 + 890534.23*T0)) +
         0.21*sin(q*(269.9 + 954397.7*T0)) - 0.19*sin(q*(357.5 + 35999.05*T0)) -
         0.11*sin(q*(186.6 + 966404.05*T0));
    cout << "Lon=" << Ekl_Lon << endl;
  double Ekl_Lat = 5.13*sin(q*(93.3 + 483202.03*T0)) + 0.28*sin(q*(228.2 + 960400.87*T0)) -
         0.28*sin(q*(318.3 + 6003.18*T0)) - 0.17*sin(q*(217.6 - 407332.2*T0));
	  cout << "Lat=" << Ekl_Lat << endl;
  double  p = 0.9508 + 0.0518*cos(q*(134.9 + 477198.85*T0)) + 0.0095*cos(q*(259.2-413335.38*T0)) +
         0.0078*cos(q*(235.7 + 890534.23*T0)) + 0.0028*cos(q*(269.9 + 954397.7*T0));
	  cout << "p=" << p << endl;
}

Не понимаю

StormBS · « **Ответ #12 :** 13 Окт 2010 [15:50:16] »

Спасибо за помощь, разобрался. Страшно ступил, директива #define не поддерживает арифметические операции. Пришлось подправить код следующим образом вместо #define Rad(x) x*M_PI/180 подставил функцию
inline double Rad(double x){return x*M_PI/180; }; И все сошлось

StormBS · « **Ответ #13 :** 15 Окт 2010 [16:57:45] »

Осталось еще несколько вопросов.
1. Как посчитать расстояние от точки наблюдения до Луны
2. Как вычислить из ранее полученных данных Пр. восх. и склонение
3. Азимут и Угол Места(высота) от точки наблюдения на Луну.

Tau · « **Ответ #14 :** 15 Окт 2010 [18:24:00] »

https://astronomy.ru/forum/index.php/topic,78687.0.html
Ответ# 17 прицепленный файл (расстояние от точки наблюдения до Луны - это модуль топоцентрического вектора из п.5)
азимут: ответ#49

Tau · « **Ответ #15 :** 15 Окт 2010 [20:53:18] »

Смотрите схему - как связаны между собой задача вычисления подлунной точки и задача вычисления азимута и высоты, что у вас уже вычислено и что надо еще вычислить.

StormBS · « **Ответ #16 :** 16 Окт 2010 [09:29:24] »

Спасибо Tau, сегодня буду разбираться.

StormBS · « **Ответ #17 :** 16 Окт 2010 [12:34:56] »

Цитата: Tau от 15 Окт 2010 [20:53:18]

Смотрите схему - как связаны между собой задача вычисления подлунной точки и задача вычисления азимута и высоты, что у вас уже вычислено и что надо еще вычислить.

По схеме написал функцию, но к сожалению большая погрешность в расстоянии от точки наблюдения до луны (сравнивал с орбитроном), а также в расчете пр. восхождения и склонения. Код функции приложен снизу.

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

void GSanSystem::LunePos(
			double  lon, //Долгота	
			double  lat,//Широта	
			double  Year,//год
			double  Mon,// месяцы
			double  Day,// дни
			double  hous,//часы
			double  min,//минуты
			double sec,// секунды
			double zona// Часовой пояс
		   )

{
	
		//UT - всемирное время в часах, момент расчета
		double UT = hous-zona + min/60 + sec/3600;
		if(UT>24)UT=UT-24;
		if(UT<0)UT=UT+24;
 		
		//Вычисление модифицированной  юлианской даты на начало суток

		 double Var1,Var2,Var3;
		 Var1 = 10000 * Year + 100 * Mon + Day;
		 if(Mon <= 2 )
		 {
			Mon = Mon + 12;
			Year = Year - 1;
		 }
		 if( Var1 <= 15821004) Var2 = -2 + floor((Year + 4716) / 4) - 1179;
		 else  Var2 = floor(Year /400) - floor(Year / 100) + floor(Year / 4);
		 Var3 = 365 * Year - 679004;
		 // MJD - Модифицированная Юлианская дата
		double MJD = Var3 + Var2 + floor(306001 * (Mon + 1)/ 10000) + Day+UT/24;
		double T0 = (MJD - 51544.5) / 36525; // мод.юл.дата на начало суток в юлианских столетиях
		// Вычисление эклиптических координат Луны и вектора (X, Y, Z) в эклиптической системе координат.
		// эклиптическая долгота
		double Ekl_Lon = 218.32 + 481267.883*T0 + 6.29*sin(Rad(134.9 + 477198.85*T0)) - 
			1.27*sin(Rad(259.2-413335.38*T0)) + 0.66*sin(Rad(235.7 + 890534.23*T0)) +
			0.21*sin(Rad(269.9 + 954397.7*T0)) - 0.19*sin(Rad(357.5 + 35999.05*T0)) - 
			0.11*sin(Rad(186.6 + 966404.05*T0));
		// эклиптическая широта
		double Ekl_Lat = 5.13*sin(Rad(93.3 + 483202.03*T0)) + 0.28*sin(Rad(228.2 + 960400.87*T0)) - 
			0.28*sin(Rad(318.3 + 6003.18*T0)) - 0.17*sin(Rad(217.6 - 407332.2*T0));
		// Горизонтальный паралакс
		double p = 0.9508 + 0.0518*cos(Rad(134.9 + 477198.85*T0)) + 0.0095*cos(Rad(259.2-413335.38*T0)) + 
			0.0078*cos(Rad(235.7 + 890534.23*T0)) + 0.0028*cos(Rad(269.9 + 954397.7*T0));
		// Растояние до Луны в радиусах Земли
		double R = 1/sin(Rad(p));
		double ER = 6378.136 ; // км радиус земли
		double Rast = ER*R; // растояние до Луны
		double X = Rast*cos(Rad(Ekl_Lon))*cos(Rad(Ekl_Lat));
		double Y = Rast*sin(Rad(Ekl_Lon))*cos(Rad(Ekl_Lat));
		double Z = Rast*sin(Rad(Ekl_Lat));
		
		//Вычисление вектора Луны (X`,Y`,Z`) в экваториальной системе координат.
		double e = (84381.488 - 46.815*T0 - 0.00059*T0*T0 + 0.001813*T0*T0*T0)/3600;
		double X’ = X;
		double Y’ = Y*cos(Rad(e)) - Z*sin(Rad(e));
		double Z’ = Y*sin(Rad(e)) + Z*cos(Rad(e));

				
		//Вычисление среднего звездного времени на дату MD и время UT.
		T0 = (int(MJD)-51544.5)/36525;
		double a1 = 24110.54841;
		double a2 = 8640184.812;
		double a3 = 0.093104;
		double a4 = 0.0000062;
		double S0 = a1 + a2 * T0 + a3 * T0*T0- a4 * T0 *T0*T0;// звездное время в Гринвиче на начало суток в секундах
		double w = 1.002737909350795;
		double Ss = S0 + UT*86400*w/24; // звездное время в секундах
		double S = Ss/240 ; // звездное время в градусах
		double Sm =  S + lon; // Местное звезное время

		//Вычисление геоцентрического вектора наблюдателя (Xn, Yn, Zn) на момент звездного времени S.
		double H = 0; //высота точки наблюдения над уровнем моря
		double Rz = 6378.14 ;// Радиус земли, км.
		double f = 1/298.26;// Сжатие Земли
		double C = 1/(sqrt(1 + f*(f-2)*sin(Rad(lat))*sin(Rad(lat))));
		double W = C*(1-f)*(1-f);

		double Xn = (R*C + H)*cos(Rad(S))*cos(Rad(lat));
		double Yn = (R*C + H)*sin(Rad(S))*cos(Rad(lat));
		double Zn = (R*W + H)*sin(Rad(lat));

		// Вычисление топоцентрического вектора луны
		double Xt = X’ - Xn;
		double Yt = Y’ - Yn;
		double Zt = Z’ - Zn;

		// Вычисление растояния от точки наблюдения до луны
		double Dist = sqrt(Xt*Xt + Yt*Yt + Zt*Zt);

		// Вычисление топоцентрическое прямое восхождение и склонение
		// RA_L - прямое восхождение Луны  на нужный момент времени
		//DEC_L - склонение Луны на нужный момент времени
		double RA_L = atan2(Yt,Xt)*180/M_PI;
		double DEC_L = atan2(Zt,sqrt(Xt*Xt + Yt*Yt))*180/M_PI;
		
		// Вычисление Азимута и Угла места лун
		double Th = Sm - RA_L  ;//часовой угол
		double z  = acos(sin(Rad(lat)) * sin(Rad(DEC_L)) + cos(Rad(lat)) * cos(Rad(DEC_L)) * cos(Rad(Th)))*180/M_PI;// косинус зенитного угла
		H = 90 - z; // угол места 
		Az = atan2( sin(Rad(Th)) * cos(Rad(DEC_L)) * cos(Rad(lat)),sin(Rad(H)) * sin(Rad(lat)) - sin(Rad(DEC_L)))*180/M_PI+180; // азимут

	
		//Перевод вектора Луны (X’, Y’, Z’) в гринвическую систему координат.
		double Xg = X’*cos(Rad(S)) + Y’*sin(Rad(S));
		double Yg = -X’*sin(Rad(S)) + Y’*cos(Rad(S));
		double Zg = Z’;
		//S – звездное гринвическое время

		//Разложение вектора Луны в гринвической системе координат на углы.
		LonLune = atan2(Yg,Xg)*180/M_PI;
		LatLune = atan2(Zg,sqrt(Xg*Xg + Yg*Yg))*180/M_PI;
		/////////////////////////////////////////////////////////////////////////////////////////



}

Sleepwalker · « **Ответ #18 :** 16 Окт 2010 [13:50:35] »

Вызывает сомнения вот этот блок:

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

double Xn = (R*C + H)*cos(Rad(S))*cos(Rad(lat));
double Yn = (R*C + H)*sin(Rad(S))*cos(Rad(lat));
double Zn = (R*W + H)*sin(Rad(lat));

// Вычисление топоцентрического вектора луны
double Xt = X’ - Xn;
double Yt = Y’ - Yn;
double Zt = Z’ - Zn;

UPD:
Не R а Rz

StormBS · « **Ответ #19 :** 16 Окт 2010 [14:07:47] »

Замену сделал. Результат уже ближе к истине.

Новости:

A A A A Автор Тема: Расчет положения солнца и луны (Прочитано 18172 раз)

Tau

Tau

Tau

Tau

Tau

Tau

Tau