Calculate sunrise and sunset, showing time and azimuthdev

Add sun trace table, and put all data into file

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 <cstdio>
+#include <cstdlib>
 #include <ctime>
 #include <iomanip>
 
-/*以下为"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<double, double> getSunPos(time_t time) {
-    double t = julian.getJulianKiloYear(time);
+FILE *fp;
 
+pair<double, double> getSunPos(double t, double flag = false) {
     pair<double, double> 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 <YYYY-MM-DD,HH:MM:SS> "
-               "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<double, double> 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<double, double> 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 <utility>
 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<int>(degrees);
     double min_sec = (degrees - deg) * 60.0;
     int min = static_cast<int>(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<double, double> parameter::sun_longitude(double t) {
     return make_pair(L, -B);
 }
 
-// 计算月球地心黄经，外部调用，返回角度制
-double parameter::moon_longitude(double t) {
-    double T = 10 * t;
-    vector<double> 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; // 返回弧度制
 }