1 files changed, 90 insertions, 30 deletions

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