From b1be91bc7562571bcb156d20ec3f9c8b69340722 Mon Sep 17 00:00:00 2001 From: We-unite <3205135446@qq.com> Date: Fri, 8 Mar 2024 22:53:11 +0800 Subject: Calculate sunrise and sunset, showing time and azimuth Add sun trace table, and put all data into file --- src/CMakeLists.txt | 6 +-- src/Julian.cpp | 3 +- src/calendar.cpp | 120 +++++++++++++++++++++++++++++++++++++++-------------- src/calendar.h | 42 +++---------------- src/parameters.cpp | 56 ++++--------------------- 5 files changed, 107 insertions(+), 120 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 53ddadc..14401b8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -5,10 +5,6 @@ add_definitions(-w) message("Project root: ${PROJECT_SOURCE_DIR}") set(EXECUTABLE_OUTPUT_PATH ${PROJECT_SOURCE_DIR}/bin) -set(LIBRARY_OUTPUT_PATH ${PROJECT_SOURCE_DIR}/lib) -# 将Julian.cpp与parameters.cpp编译为一个动态库 -add_library(parameters SHARED Julian.cpp parameters.cpp) # 将calendar.cpp与List.cpp编译为可执行文件，链接动态库 -add_executable(calendar calendar.cpp) -target_link_libraries(calendar PRIVATE parameters) +add_executable(calendar calendar.cpp Julian.cpp parameters.cpp) \ No newline at end of file diff --git a/src/Julian.cpp b/src/Julian.cpp index 5a0522b..ecfe391 100644 --- a/src/Julian.cpp +++ b/src/Julian.cpp @@ -86,6 +86,5 @@ double Julian::getJulianKiloYear(time_t time) { //儒略千年数转时间戳 time_t Julian::kiloYearToTime(double t, int year) { double jd = t * 365250.0 + 2451545.0; - time_t time = (time_t)((jd - 2440587.5) * 86400.0 - delta_t(year)); - return time; + return (time_t)((jd - 2440587.5) * 86400.0 - delta_t(year)); } diff --git a/src/calendar.cpp b/src/calendar.cpp index 6731692..7e04857 100644 --- a/src/calendar.cpp +++ b/src/calendar.cpp @@ -1,64 +1,124 @@ // Encoding: UTF-8 #include "calendar.h" #include +#include #include #include -/*以下为"calendar.h规定之量*/ -double day = 86400; -double delta = 1e-11; +// 常量 +#define day_secs 86400 +#define delta_time 1e-11 + +/* calendar.h规定之量 */ Julian julian; parameter p; +double epsilon; -pair getSunPos(time_t time) { - double t = julian.getJulianKiloYear(time); +FILE *fp; +pair getSunPos(double t, double flag = false) { pair sun = p.sun_longitude(t); - radToDMS("Sun longitude:\t\t", 0, 360, sun.first); - radToDMS("Sun latitude:\t\t", -90, 90, sun.second); double alpha, delta; // 太阳赤经赤纬 double lambda = sun.first, beta = sun.second; // 太阳黄经黄纬 // 黄赤交角epsilon - double epsilon = p.get_epsilon(t); - radToDMS("Ecliptic Obliquity:\t\t", 0, 90, epsilon); + epsilon = p.get_epsilon(t); delta = asin(sin(epsilon) * sin(lambda) * cos(beta) + cos(epsilon) * sin(beta)); alpha = atan2( (cos(epsilon) * sin(lambda) * cos(beta) - sin(epsilon) * sin(beta)), (cos(lambda) * cos(beta))); - radToDMS("Sun right ascension:\t\t", 0, 360, alpha); - radToDMS("Sun declination:\t\t", -90, 90, delta); // 时角 double H = p.getHourAngle(t, LONGITUDE, LATITUDE, alpha); double A, h; //方位角和高度角 h = asin(sin(LATITUDE) * sin(delta) + cos(LATITUDE) * cos(delta) * cos(H)); A = atan2(-cos(delta) * sin(H), cos(LATITUDE) * sin(delta) - sin(LATITUDE) * cos(delta) * cos(H)); - A += M_PI; - // 转换为角度制输出 - radToDMS("Hour angle:\t\t", 0, 360, H); - radToDMS("Azimuth:\t\t", 0, 360, A); - radToDMS("Elevation:\t\t", -90, 90, h); + if (flag) { + fprintf( + fp, "%s\t%s\t%s\t%s\t%s\t%s\n", radToDMS(0, 360, lambda).c_str(), + radToDMS(-90, 90, beta).c_str(), radToDMS(0, 360, alpha).c_str(), + radToDMS(-90, 90, delta).c_str(), radToDMS(0, 360, A).c_str(), + radToDMS(-90, 90, h).c_str()); + } + return make_pair(A, h + M_PI / 216); } -int main(int argc, char *argv[]) { - Date date; - if (argc != 2) { - printf("Input the time you want to calculate in " - "format:\t\t"); - scanf("%d-%d-%d,%d:%d:%d", &date.tm_year, &date.tm_mon, &date.tm_mday, - &date.tm_hour, &date.tm_min, &date.tm_sec); - } else { - sscanf(argv[1], "%d-%d-%d,%d:%d:%d", &date.tm_year, &date.tm_mon, - &date.tm_mday, &date.tm_hour, &date.tm_min, &date.tm_sec); +void binarySearch(string message, time_t begin, Date date) { + double t = julian.getJulianKiloYear(begin); + double t_end = julian.getJulianKiloYear(begin + day_secs / 2); + double middle; + time_t res; + pair sunPos; + while (t_end - t >= delta_time) { + middle = (t + t_end) / 2; + sunPos = getSunPos(middle, false); + if (sunPos.second * getSunPos(t).second > 0) { + t = middle; + } else { + t_end = middle; + } + } + res = julian.kiloYearToTime(middle, date.tm_year + 1900); + Date *tm = localtime(&res); + printf("%s\t%4d-%02d-%02d %02d:%02d:%02d\tAzimuth\t%s\n", message.c_str(), + tm->tm_year + 1900, tm->tm_mon + 1, tm->tm_mday, tm->tm_hour, + tm->tm_min, tm->tm_sec, radToDMS(0, 360, sunPos.first).c_str()); + fprintf(fp, "%s\t%4d-%02d-%02d %02d:%02d:%02d\tAzimuth\t%s\n", + message.c_str(), tm->tm_year + 1900, tm->tm_mon + 1, tm->tm_mday, + tm->tm_hour, tm->tm_min, tm->tm_sec, + radToDMS(0, 360, sunPos.first).c_str()); +} + +void getSunTime(Date date) { + time_t time = mktime(&date); + binarySearch("Sunrize", time, date); + binarySearch("Sunset", time + day_secs / 2, date); +} + +void showSunPos(Date date) { + time_t time = mktime(&date); + double t; + fprintf(fp, "\n\nSun Trace Table\n\n"); + fprintf(fp, "--------------------------------------------------------------" + "--------------------------------------------------------------" + "----------\n"); + fprintf(fp, "Time |\tlongitude\t\tlatitude\t\tright " + "ascension\tdeclination\tAzimuth\t\tAltitude\n"); + fprintf(fp, "--------------------------------------------------------------" + "--------------------------------------------------------------" + "----------\n"); + for (int i = 0; i < 24; i++, time += 3600, date.tm_hour++) { + t = julian.getJulianKiloYear(time); + fprintf(fp, "%02d:%02d:%02d |\t", date.tm_hour, date.tm_min, + date.tm_sec); + getSunPos(t, true); } + fprintf(fp, "--------------------------------------------------------------" + "--------------------------------------------------------------" + "----------\n"); +} + +int main() { + Date date; + cout << "Input a day you want(in YYYY-MM-DD): "; + scanf("%d-%d-%d", &date.tm_year, &date.tm_mon, &date.tm_mday); date.tm_year -= 1900; date.tm_mon -= 1; - date.tm_isdst = -1; + date.tm_hour = date.tm_min = date.tm_sec = 0; + + fp = fopen("sunPos.txt", "w"); + if (fp == NULL) { + printf("File open failed\n"); + exit(1); + } + fprintf(fp, "Sun trace File (Date: %4d-%02d-%02d)\n\n\n", + date.tm_year + 1900, date.tm_mon + 1, date.tm_mday); + getSunTime(date); + + showSunPos(date); + fclose(fp); - time_t time = mktime(&date); - getSunPos(time); return 0; -} +} \ No newline at end of file diff --git a/src/calendar.h b/src/calendar.h index 248026f..87ca413 100644 --- a/src/calendar.h +++ b/src/calendar.h @@ -10,43 +10,13 @@ using namespace std; // 哈尔滨南岗地方经纬度 -#define LONGITUDE 126.533 * M_PI / 180 -#define LATITUDE 45.800 * M_PI / 180 - -extern double delta; -extern char jieqi[25][10]; +#define M_PI 3.14159265358979323846 +#define LONGITUDE 126.66854 * M_PI / 180 +#define LATITUDE 45.75996 * M_PI / 180 typedef struct tm Date; -struct point { - int year, mon, day, hour, min, sec; - bool isShuo, isJieqi, isZhongqi; - int JieqiIndex, MonthIndex; - bool RunYue; - double time; - point *next; - - point(int year, int mon, int day, int hour, int min, int sec, bool isShuo, - bool isJieqi, bool isZhongqi, int JieqiIndex, int MonthIndex, - bool RunYue, double time) { - this->year = year; - this->mon = mon; - this->day = day; - this->hour = hour; - this->min = min; - this->sec = sec; - this->isShuo = isShuo; - this->isJieqi = isJieqi; - this->isZhongqi = isZhongqi; - this->JieqiIndex = JieqiIndex; - this->MonthIndex = MonthIndex; - this->RunYue = RunYue; - this->time = time; - this->next = NULL; - } -}; -void radToDMS(const char *message, int lowerBound, int upperBound, - double radians); +string radToDMS(int lowerBound, int upperBound, double radians); class Julian { private: @@ -96,8 +66,6 @@ class parameter { public: pair sun_longitude(double t); - double moon_longitude(double t); - double get_epsilon(double T); // 计算时角，返回弧度制 @@ -107,4 +75,6 @@ class parameter { extern Julian julian; extern parameter p; +extern double epsilon; + #endif \ No newline at end of file diff --git a/src/parameters.cpp b/src/parameters.cpp index fdcefd1..ff57305 100644 --- a/src/parameters.cpp +++ b/src/parameters.cpp @@ -5,8 +5,7 @@ #include using namespace std; -void radToDMS(const char *message, int lowerBound, int upperBound, - double radians) { +string radToDMS(int lowerBound, int upperBound, double radians) { // Convert radians to degrees double degrees = radians * 180.0 / M_PI; @@ -15,17 +14,19 @@ void radToDMS(const char *message, int lowerBound, int upperBound, while (degrees >= upperBound) degrees -= (upperBound - lowerBound); - bool flag = degrees < 0; - if (flag) - degrees = -degrees; // Calculate degrees, minutes, and seconds int deg = static_cast(degrees); double min_sec = (degrees - deg) * 60.0; int min = static_cast(min_sec); double sec = (min_sec - min) * 60.0; + if (degrees < 0) { + min = -min; + sec = -sec; + } - // Output the result - printf("%s: %c%d°%d'%06.3f''\n", message, flag ? '-' : ' ', deg, min, sec); + char buf[100]; + sprintf(buf, "%4d°%2d'%06.3f''", deg, min, sec); + return string(buf); } // 计算地球日心黄经 @@ -205,51 +206,12 @@ pair parameter::sun_longitude(double t) { return make_pair(L, -B); } -// 计算月球地心黄经，外部调用，返回角度制 -double parameter::moon_longitude(double t) { - double T = 10 * t; - vector angles = get_angle(T); - double a, b, c, d, moon_longitude_a; - - double L = 218.3164591 + 481267.88134236 * T - 0.0013268 * T * T + - T * T * T / 538841 - T * T * T * T / 65194000; //月球平黄经 - double EI = 0, sita = 0; - double E = 1 - 0.002516 * T - 0.0000074 * T * T; //地球离心率 - - for (int i = 0; i < moon_size; i++) { - sita = moon_longitude_parameters[i][0] * angles[0] + - moon_longitude_parameters[i][1] * angles[1] + - moon_longitude_parameters[i][2] * angles[2] + - moon_longitude_parameters[i][3] * angles[3]; - sita *= M_PI / 180; - EI += moon_longitude_parameters[i][4] * sin(sita) * - pow(E, fabs(moon_longitude_parameters[i][1])); - } - - //地外行星影响月球地心黄经修正项 - double a1, a2; - a1 = 119.75 + 131.849 * T; - a2 = 53.09 + 479264.290 * T; - - EI += 3958.0 * sin(a1 * M_PI / 180); - EI += 1962.0 * sin((L - angles[3]) * M_PI / 180); - EI += 318.0 * sin(a2 * M_PI / 180); - - L += EI / 1000000; - L += 180 / M_PI * nutation(T); //地球章动修正 - - while (L < 0) - L += 360; - while (L >= 360) - L -= 360; - return L; -} - double parameter::getHourAngle(double julianKiloYear, double longitude, double latitude, double alpha) { double T = julianKiloYear * 10; double GST = 280.46061837 + 360.98564736629 * 36525 * T + 0.000387933 * T * T - T * T * T / 38710000; + GST += nutation(T) * cos(epsilon); GST *= M_PI / 180; return GST + longitude - alpha; // 返回弧度制 } -- cgit v1.2.3-70-g09d2