From 6fd843db21af18be0634bac770c6b74ad4e78bb7 Mon Sep 17 00:00:00 2001 From: We-unite <3205135446@qq.com> Date: Fri, 8 Mar 2024 21:24:55 +0800 Subject: HourAngle Fixed, and sun pos right --- src/.clang-format | 121 ----------------------------------------------------- src/calendar.cpp | 77 +++++++++++++--------------------- src/calendar.h | 20 ++++++--- src/parameters.cpp | 45 ++++++++++++++++---- 4 files changed, 80 insertions(+), 183 deletions(-) delete mode 100644 src/.clang-format (limited to 'src') 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 getSunPos(time_t time) { + double t = julian.getJulianKiloYear(time); - while (degrees < lowerBound) - degrees += (upperBound - lowerBound); - while (degrees >= upperBound) - degrees -= (upperBound - lowerBound); + pair 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(degrees); - double min_sec = (degrees - deg) * 60.0; - int min = static_cast(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 " - "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 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 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(degrees); + double min_sec = (degrees - deg) * 60.0; + int min = static_cast(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 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 -- cgit v1.2.3-70-g09d2