From 505d07a7abccfc707e2d27c2deb2fe460854e751 Mon Sep 17 00:00:00 2001
From: We-unite <3205135446@qq.com>
Date: Thu, 7 Mar 2024 21:50:22 +0800
Subject: Confused about hourAngle

---
 src/calendar.cpp   | 24 ++++++++++++++++++++++++
 src/calendar.h     |  9 ++++++++-
 src/parameters.cpp | 13 ++++++++++++-
 3 files changed, 44 insertions(+), 2 deletions(-)

diff --git a/src/calendar.cpp b/src/calendar.cpp
index 9ad8413..8af1851 100644
--- a/src/calendar.cpp
+++ b/src/calendar.cpp
@@ -53,5 +53,29 @@ int main(int argc, char *argv[]) {
     pair<double, double> sun = p.sun_longitude(t);
     radToDMS("Sun longitude: ", 0, 360, sun.first);
     radToDMS("Sun latitude: ", -90, 90, sun.second);
+
+    double alpha, delta;                          // 太阳赤经赤纬
+    double lambda = sun.first, beta = sun.second; // 太阳黄经黄纬
+    // 黄赤交角epsilon，23°26'21.448''转换为弧度
+    double epsilon = (23.0 + 26.0 / 60.0 + 21.448 / 3600.0) * M_PI / 180.0;
+
+    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(epsilon) * cos(beta)));
+    radToDMS("Sun right ascension: ", 0, 360, alpha);
+    radToDMS("Sun declination: ", -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: ", 0, 360, H);
+    radToDMS("Azimuth: ", 0, 360, A);
+    radToDMS("Elevation: ", 0, 90, h);
     return 0;
 }
diff --git a/src/calendar.h b/src/calendar.h
index ef785c0..810fa0c 100644
--- a/src/calendar.h
+++ b/src/calendar.h
@@ -9,6 +9,10 @@
 #include <vector>
 using namespace std;
 
+// 哈尔滨南岗地方经纬度
+#define LONGITUDE 126.68 * M_PI / 180
+#define LATITUDE 45.75 * M_PI / 180
+
 extern double delta;
 extern char jieqi[25][10];
 
@@ -90,6 +94,9 @@ class parameter {
     pair<double, double> sun_longitude(double t);
 
     double moon_longitude(double t);
+
+    double getHourAngle(double julianKiloTime, double longitude,
+                        double latitude, double alpha);
 };
 
-#endif
+#endif
\ No newline at end of file
diff --git a/src/parameters.cpp b/src/parameters.cpp
index 2e2758a..f08736b 100644
--- a/src/parameters.cpp
+++ b/src/parameters.cpp
@@ -1,5 +1,5 @@
-#include "calendar.h"
 #include "parameters.h"
+#include "calendar.h"
 #include <utility>
 using namespace std;
 
@@ -221,3 +221,14 @@ double parameter::moon_longitude(double t) {
         L -= 360;
     return L;
 }
+
+double parameter::getHourAngle(double julianKiloTime, 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;
+}
\ No newline at end of file
-- 
cgit v1.2.3-70-g09d2