//#include "StdAfx.h"
//#include "GSanSystem.h"
//#include "math.h"
#include <iostream>
#include <cmath>
#define M_PI 3.14159265358979323846
#define Rad(x) x*M_PI/180.0
using namespace std;
//GSanSystem::GSanSystem(void){}
//GSanSystem::~GSanSystem(void){}
/*int round(double n) {
double t;
t = n - double(n);
n*=10;
if (t >= 0.5) {
ceil(n);
} else {
floor(n);
}
n/=10;
return (int)n;
}*/
int main() {
//1.Вычисление модифицированной юлианской даты на начало суток
double lon, //Долгота
lat,//Широта
Year,//год
Mon,// месяцы
Day,// дни
hour,//часы
min,//минуты
sec,// секунды
zona;// Часовой пояс
cin >> lon >> lat >> Year >> Mon >> Day >> hour >> min >> sec >> zona;
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 = ((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;// звездное время в Гринвиче на начало суток в секундах
cout << S0 << endl;
//UT - всемирное время в часах, момент расчета
double UT = hour-zona + min/60 + sec/3600;
if (UT > 24) UT -= 24;
if (UT < 0) UT += 24;
double Nsec = UT * 3600; // количество секунд, прошедших от начала суток до момента наблюдения
double NsecS = Nsec * 366.2422 / 365.2422; //количество звездных секунд
cout << NsecS << endl;
double GMT = (S0 + NsecS) /3600 * 15;//гринвическое среднее звездное время в градусах SG
while (GMT > 360) 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 -= 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 -= 360;
double X = cos(Rad(L)) ; // вектор
double Y = sin(Rad(L)) ; // в эклиптической
double Z = 0 ; // системе координат
cout << X << endl << Y << endl << Z << endl;
// Координаты Cолнца в прямоугольной экваториальной системе координат
double Eps = 23.439281 ; //наклон эклиптики к экватору
double X1 = X ; // вектор
double Y1 = Y * cos(Rad(Eps)) - Z * sin(Rad(Eps)); // в экваториальной
double Z1 = Y * sin(Rad(Eps)) + Z * cos(Rad(Eps)) ;// системе координат
// Экваториальные геоцентрические координаты Солнца
// RA - прямое восхождение Солнца на нужный момент времени
//DEC - склонение Солнца на нужный момент времени
double Ra = atan2(Y1 ,X1)*180/M_PI;
double Dec = atan2 (Z1 , sqrt(X1*X1 + Y1*Y1))*180/M_PI;
cout << Ra << endl << Dec << endl;
// Азимутальные координаты Солнца
//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;
Az = Az + 180; // Переводим в геодезический азимут
// получаем подсолнечную точку
// Долгота Солнца
double LonSan = Ra - GST;
// Широта Солнца
double LatSan = Dec;
cout << H << endl << Az << endl;
}