1 files changed, 24 insertions, 0 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;
 }