HourAngle Fixed, and sun pos right

4 files changed, 80 insertions, 183 deletions

diff --git a/src/.clang-format b/src/.clang-format
deleted file mode 100644
index a1054a5..0000000
--- a/src/.clang-format
+++ /dev/null
@@ -1,121 +0,0 @@
----
-Language:        Cpp
-# BasedOnStyle:  LLVM
-AccessModifierOffset: -2
-AlignAfterOpenBracket: Align
-AlignConsecutiveAssignments: false
-AlignConsecutiveDeclarations: false
-AlignEscapedNewlines: Right
-AlignOperands:   true
-AlignTrailingComments: true
-AllowAllParametersOfDeclarationOnNextLine: true
-AllowShortBlocksOnASingleLine: false
-AllowShortCaseLabelsOnASingleLine: false
-AllowShortFunctionsOnASingleLine: All
-AllowShortIfStatementsOnASingleLine: false
-AllowShortLoopsOnASingleLine: false
-AlwaysBreakAfterDefinitionReturnType: None
-AlwaysBreakAfterReturnType: None
-AlwaysBreakBeforeMultilineStrings: false
-AlwaysBreakTemplateDeclarations: MultiLine
-BinPackArguments: true
-BinPackParameters: true
-BraceWrapping:   
-  AfterClass:      false
-  AfterControlStatement: false
-  AfterEnum:       false
-  AfterFunction:   false
-  AfterNamespace:  false
-  AfterObjCDeclaration: false
-  AfterStruct:     false
-  AfterUnion:      false
-  AfterExternBlock: false
-  BeforeCatch:     false
-  BeforeElse:      false
-  IndentBraces:    false
-  SplitEmptyFunction: true
-  SplitEmptyRecord: true
-  SplitEmptyNamespace: true
-BreakBeforeBinaryOperators: None
-BreakBeforeBraces: Attach
-BreakBeforeInheritanceComma: false
-BreakInheritanceList: BeforeColon
-BreakBeforeTernaryOperators: true
-BreakConstructorInitializersBeforeComma: false
-BreakConstructorInitializers: BeforeColon
-BreakAfterJavaFieldAnnotations: false
-BreakStringLiterals: true
-ColumnLimit:     80
-CommentPragmas:  '^ IWYU pragma:'
-CompactNamespaces: false
-ConstructorInitializerAllOnOneLineOrOnePerLine: false
-ConstructorInitializerIndentWidth: 4
-ContinuationIndentWidth: 4
-Cpp11BracedListStyle: true
-DerivePointerAlignment: false
-DisableFormat:   false
-ExperimentalAutoDetectBinPacking: false
-FixNamespaceComments: true
-ForEachMacros:   
-  - foreach
-  - Q_FOREACH
-  - BOOST_FOREACH
-IncludeBlocks:   Preserve
-IncludeCategories: 
-  - Regex:           '^"(llvm|llvm-c|clang|clang-c)/'
-    Priority:        2
-  - Regex:           '^(<|"(gtest|gmock|isl|json)/)'
-    Priority:        3
-  - Regex:           '.*'
-    Priority:        1
-IncludeIsMainRegex: '(Test)?$'
-IndentCaseLabels: false
-IndentPPDirectives: None
-IndentWidth:     4
-IndentWrappedFunctionNames: false
-JavaScriptQuotes: Leave
-JavaScriptWrapImports: true
-KeepEmptyLinesAtTheStartOfBlocks: true
-MacroBlockBegin: ''
-MacroBlockEnd:   ''
-MaxEmptyLinesToKeep: 1
-NamespaceIndentation: None
-ObjCBinPackProtocolList: Auto
-ObjCBlockIndentWidth: 2
-ObjCSpaceAfterProperty: false
-ObjCSpaceBeforeProtocolList: true
-PenaltyBreakAssignment: 2
-PenaltyBreakBeforeFirstCallParameter: 19
-PenaltyBreakComment: 300
-PenaltyBreakFirstLessLess: 120
-PenaltyBreakString: 1000
-PenaltyBreakTemplateDeclaration: 10
-PenaltyExcessCharacter: 1000000
-PenaltyReturnTypeOnItsOwnLine: 60
-PointerAlignment: Right
-ReflowComments:  true
-SortIncludes:    true
-SortUsingDeclarations: true
-SpaceAfterCStyleCast: false
-SpaceAfterTemplateKeyword: true
-SpaceBeforeAssignmentOperators: true
-SpaceBeforeCpp11BracedList: false
-SpaceBeforeCtorInitializerColon: true
-SpaceBeforeInheritanceColon: true
-SpaceBeforeParens: ControlStatements
-SpaceBeforeRangeBasedForLoopColon: true
-SpaceInEmptyParentheses: false
-SpacesBeforeTrailingComments: 1
-SpacesInAngles:  false
-SpacesInContainerLiterals: true
-SpacesInCStyleCastParentheses: false
-SpacesInParentheses: false
-SpacesInSquareBrackets: false
-Standard:        Cpp11
-StatementMacros: 
-  - Q_UNUSED
-  - QT_REQUIRE_VERSION
-TabWidth:        8
-UseTab:          Never
-...
-
diff --git a/src/calendar.cpp b/src/calendar.cpp
index 8af1851..6731692 100644
--- a/src/calendar.cpp
+++ b/src/calendar.cpp
@@ -7,36 +7,47 @@
 /*以下为"calendar.h规定之量*/
 double day = 86400;
 double delta = 1e-11;
-
-/*计算需要用的类的实例化*/
 Julian julian;
 parameter p;
 
-void radToDMS(const char *message, int lowerBound, int upperBound,
-              double radians) {
-    // Convert radians to degrees
-    double degrees = radians * 180.0 / M_PI;
+pair<double, double> getSunPos(time_t time) {
+    double t = julian.getJulianKiloYear(time);
 
-    while (degrees < lowerBound)
-        degrees += (upperBound - lowerBound);
-    while (degrees >= upperBound)
-        degrees -= (upperBound - lowerBound);
+    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);
 
-    // 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;
+    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);
 
-    // Output the result
-    printf("%s: %d°%d'%06.3f''\n", message, deg, min, sec);
+    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);
 }
 
 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: ");
+               "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 {
@@ -48,34 +59,6 @@ int main(int argc, char *argv[]) {
     date.tm_isdst = -1;
 
     time_t time = mktime(&date);
-    double t = julian.getJulianKiloYear(time);
-
-    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);
+    getSunPos(time);
     return 0;
 }
diff --git a/src/calendar.h b/src/calendar.h
index 810fa0c..248026f 100644
--- a/src/calendar.h
+++ b/src/calendar.h
@@ -10,8 +10,8 @@
 using namespace std;
 
 // 哈尔滨南岗地方经纬度
-#define LONGITUDE 126.68 * M_PI / 180
-#define LATITUDE 45.75 * M_PI / 180
+#define LONGITUDE 126.533 * M_PI / 180
+#define LATITUDE 45.800 * M_PI / 180
 
 extern double delta;
 extern char jieqi[25][10];
@@ -45,15 +45,15 @@ struct point {
     }
 };
 
+void radToDMS(const char *message, int lowerBound, int upperBound,
+              double radians);
+
 class Julian {
   private:
     static double d[23][5];
 
     double dt_ext(int y, double jsd);
 
-    //计算力学时与世界时之差,传入年份
-    double delta_t(int y);
-
   public:
     // 计算儒略日
     double getJulianDay(time_t time);
@@ -63,6 +63,9 @@ class Julian {
 
     //儒略千年数转时间戳
     time_t kiloYearToTime(double t, int year);
+
+    //计算力学时与世界时之差,传入年份
+    double delta_t(int y);
 };
 
 class parameter {
@@ -95,8 +98,13 @@ class parameter {
 
     double moon_longitude(double t);
 
-    double getHourAngle(double julianKiloTime, double longitude,
+    double get_epsilon(double T);
+
+    // 计算时角，返回弧度制
+    double getHourAngle(double julianKiloYear, double longitude,
                         double latitude, double alpha);
 };
 
+extern Julian julian;
+extern parameter p;
 #endif
\ No newline at end of file
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