HourAngle Fixed, and sun pos right

1 files changed, 36 insertions, 9 deletions

diff --git a/src/parameters.cpp b/src/parameters.cpp
index f08736b..fdcefd1 100644
--- a/src/parameters.cpp
+++ b/src/parameters.cpp
@@ -1,9 +1,32 @@
+// 凡无特殊说明的函数，全部返回弧度制
+
 #include "parameters.h"
 #include "calendar.h"
 #include <utility>
 using namespace std;
 
-//凡无特殊说明的函数，全部返回弧度制
+void radToDMS(const char *message, int lowerBound, int upperBound,
+              double radians) {
+    // Convert radians to degrees
+    double degrees = radians * 180.0 / M_PI;
+
+    while (degrees < lowerBound)
+        degrees += (upperBound - lowerBound);
+    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;
+
+    // Output the result
+    printf("%s: %c%d°%d'%06.3f''\n", message, flag ? '-' : ' ', deg, min, sec);
+}
 
 // 计算地球日心黄经
 double parameter::get_longitude(vector<double> l, double t) {
@@ -222,13 +245,17 @@ double parameter::moon_longitude(double t) {
     return L;
 }
 
-double parameter::getHourAngle(double julianKiloTime, double longitude,
+double parameter::getHourAngle(double julianKiloYear, double longitude,
                                double latitude, double alpha) {
-    double t = 100 * 365.25 * julianKiloTime;
-    // 该时刻的格林尼治恒星时
-    double GST = 6.697374558 + 2400.051336 * t + 0.000025862333 * t * t -
-                 0.000000001722222 * t * t * t;
-    return GST * M_PI / 12 + longitude - alpha;
-    // double GST = 241.654320 + 360.9856473840 * t * 1000 * 36525;
-    // return GST * M_PI / 180 + longitude - alpha;
+    double T = julianKiloYear * 10;
+    double GST = 280.46061837 + 360.98564736629 * 36525 * T +
+                 0.000387933 * T * T - T * T * T / 38710000;
+    GST *= M_PI / 180;
+    return GST + longitude - alpha; // 返回弧度制
+}
+
+double parameter::get_epsilon(double T) {
+    double epsilon = 23.439291111 - 0.01300416667 * T -
+                     (1.638888889e-7) * T * T + (5.036111111e-7) * T * T * T;
+    return epsilon * M_PI / 180;
 }
\ No newline at end of file