commit b9605183de1412d3d24f3c6c5a17eebbe0c1338c Author: liuyihui Date: Mon Nov 7 00:18:23 2022 +0800 config example of cmake(ComSim) diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..a3bd1c0 --- /dev/null +++ b/.clang-format @@ -0,0 +1,208 @@ +--- +Language: Cpp +# BasedOnStyle: Google +AccessModifierOffset: -4 +AlignAfterOpenBracket: Align +AlignArrayOfStructures: None +AlignConsecutiveMacros: None +AlignConsecutiveAssignments: None +AlignConsecutiveBitFields: None +AlignConsecutiveDeclarations: None +AlignEscapedNewlines: Left +AlignOperands: Align +AlignTrailingComments: true +AllowAllArgumentsOnNextLine: true +AllowAllConstructorInitializersOnNextLine: true +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortEnumsOnASingleLine: true +AllowShortBlocksOnASingleLine: Never +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: All +AllowShortLambdasOnASingleLine: All +AllowShortIfStatementsOnASingleLine: WithoutElse +AllowShortLoopsOnASingleLine: true +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: true +AlwaysBreakTemplateDeclarations: Yes +AttributeMacros: + - __capability +BinPackArguments: true +BinPackParameters: true +BraceWrapping: + AfterCaseLabel: false + AfterClass: false + AfterControlStatement: Never + AfterEnum: false + AfterFunction: false + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false + AfterExternBlock: false + BeforeCatch: false + BeforeElse: false + BeforeLambdaBody: false + BeforeWhile: false + IndentBraces: false + SplitEmptyFunction: true + SplitEmptyRecord: true + SplitEmptyNamespace: true +BreakBeforeBinaryOperators: None +BreakBeforeConceptDeclarations: true +BreakBeforeBraces: Attach +BreakBeforeInheritanceComma: false +BreakInheritanceList: BeforeColon +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: false +BreakConstructorInitializers: BeforeColon +BreakAfterJavaFieldAnnotations: false +BreakStringLiterals: true +ColumnLimit: 120 +CommentPragmas: '^ IWYU pragma:' +CompactNamespaces: false +ConstructorInitializerAllOnOneLineOrOnePerLine: true +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: true +DeriveLineEnding: true +DerivePointerAlignment: true +DisableFormat: false +EmptyLineAfterAccessModifier: Never +EmptyLineBeforeAccessModifier: LogicalBlock +ExperimentalAutoDetectBinPacking: false +FixNamespaceComments: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IfMacros: + - KJ_IF_MAYBE +IncludeBlocks: Regroup +IncludeCategories: + - Regex: '^[<"].*\.h[>"]' + Priority: 0 + SortPriority: 0 + CaseSensitive: false + - Regex: '^[<"].*\.hh[>"]' + Priority: 1 + SortPriority: 1 + CaseSensitive: false + - Regex: '.*' + Priority: 2 + SortPriority: 2 + CaseSensitive: false +IncludeIsMainRegex: '([-_](test|unittest))?$' +IncludeIsMainSourceRegex: '' +IndentAccessModifiers: false +IndentCaseLabels: true +IndentCaseBlocks: false +IndentGotoLabels: true +IndentPPDirectives: None +IndentExternBlock: AfterExternBlock +IndentRequires: false +IndentWidth: 4 +IndentWrappedFunctionNames: false +InsertTrailingCommas: None +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLinesAtTheStartOfBlocks: false +LambdaBodyIndentation: Signature +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +ObjCBinPackProtocolList: Never +ObjCBlockIndentWidth: 2 +ObjCBreakBeforeNestedBlockParam: true +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: true +PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 1 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyBreakTemplateDeclaration: 10 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 200 +PenaltyIndentedWhitespace: 0 +PointerAlignment: Left +PPIndentWidth: -1 +RawStringFormats: + - Language: Cpp + Delimiters: + - cc + - CC + - cpp + - Cpp + - CPP + - 'c++' + - 'C++' + CanonicalDelimiter: '' + BasedOnStyle: google + - Language: TextProto + Delimiters: + - pb + - PB + - proto + - PROTO + EnclosingFunctions: + - EqualsProto + - EquivToProto + - PARSE_PARTIAL_TEXT_PROTO + - PARSE_TEST_PROTO + - PARSE_TEXT_PROTO + - ParseTextOrDie + - ParseTextProtoOrDie + - ParseTestProto + - ParsePartialTestProto + CanonicalDelimiter: pb + BasedOnStyle: google +ReferenceAlignment: Pointer +ReflowComments: true +ShortNamespaceLines: 1 +SortIncludes: CaseSensitive +SortJavaStaticImport: Before +SortUsingDeclarations: true +SpaceAfterCStyleCast: false +SpaceAfterLogicalNot: false +SpaceAfterTemplateKeyword: true +SpaceBeforeAssignmentOperators: true +SpaceBeforeCaseColon: false +SpaceBeforeCpp11BracedList: false +SpaceBeforeCtorInitializerColon: true +SpaceBeforeInheritanceColon: true +SpaceBeforeParens: ControlStatements +SpaceAroundPointerQualifiers: Default +SpaceBeforeRangeBasedForLoopColon: true +SpaceInEmptyBlock: false +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 2 +SpacesInAngles: Never +SpacesInConditionalStatement: false +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInLineCommentPrefix: + Minimum: 1 + Maximum: -1 +SpacesInParentheses: false +SpacesInSquareBrackets: false +SpaceBeforeSquareBrackets: false +BitFieldColonSpacing: Both +Standard: Auto +StatementAttributeLikeMacros: + - Q_EMIT +StatementMacros: + - Q_UNUSED + - QT_REQUIRE_VERSION +TabWidth: 8 +UseCRLF: false +UseTab: Never +WhitespaceSensitiveMacros: + - STRINGIZE + - PP_STRINGIZE + - BOOST_PP_STRINGIZE + - NS_SWIFT_NAME + - CF_SWIFT_NAME +... + diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json new file mode 100644 index 0000000..6e92226 --- /dev/null +++ b/.vscode/c_cpp_properties.json @@ -0,0 +1,22 @@ +{ + "configurations": [ + { + "name": "Linux", + "includePath": [ + "${default}", + "${workspaceFolder}/include", + "/home/fox/.local/include", + "/home/fox/miniconda3/include", + "/home/fox/miniconda3/lib/python3.9/site-packages/numpy/core/include" + ], + "defines": [ + "_DEBUG", + "UNICODE", + "_UNICODE" + ], + "compilerPath": "/usr/bin/gcc", + "configurationProvider": "ms-vscode.cmake-tools" + } + ], + "version": 4 +} \ No newline at end of file diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 0000000..540bf89 --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,30 @@ +{ + // 使用 IntelliSense 了解相关属性。 + // 悬停以查看现有属性的描述。 + // 欲了解更多信息，请访问: https://go.microsoft.com/fwlink/?linkid=830387 + "version": "0.2.0", + "configurations": [ + { + "name": "CMake Debug", + "type": "cppdbg", + "request": "launch", + "program": "${workspaceFolder}/build/ComSim", + "args": [ + "6", + "16", + "300", + "--eps", + "0.005", + "-W", + "3000" + ], + "stopAtEntry": false, + "cwd": "${workspaceFolder}", + "environment": [], + "externalConsole": false, + "MIMode": "gdb", + "miDebuggerPath": "/usr/bin/gdb", + "preLaunchTask": "CMake Build" + } + ] +} \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..e9922c9 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,3 @@ +{ + "C_Cpp.default.configurationProvider": "ms-vscode.cmake-tools", +} \ No newline at end of file diff --git a/.vscode/tasks.json b/.vscode/tasks.json new file mode 100644 index 0000000..b954e3b --- /dev/null +++ b/.vscode/tasks.json @@ -0,0 +1,21 @@ +{ + "tasks": [ + { + "type": "shell", + "label": "CMake Build", + "command": "cmake", + "args": [ + "--build", + "build", + "--config", + "Release", + "--target", + "all", + "-j", + "18", + "-DCMAKE_BUILD_TYPE=Debug" + ] + } + ], + "version": "2.0.0" +} \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..edffa80 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,21 @@ +cmake_minimum_required(VERSION 3.16) + +project(ComSim) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +# set(CMAKE_CXX_FLAGS_DEBUG "-O0 -Wall -g -ggdb") + +include_directories(${EIGEN3_INCLUDE_DIR}) +include_directories(${PROJECT_SOURCE_DIR}/include) +include_directories("/home/fox/.local/include") + +file(GLOB sources ${PROJECT_SOURCE_DIR}/src/*.cpp) +file(GLOB headers ${PROJECT_SOURCE_DIR}/include/*.h) +add_executable(ComSim main.cpp ${sources} ${headers}) + +find_package(Eigen3 REQUIRED) +find_package(Python3 COMPONENTS Development NumPy REQUIRED) +target_link_libraries(ComSim Eigen3::Eigen Python3::Python Python3::NumPy) + +install(TARGETS ComSim DESTINATION bin) diff --git a/include/constant.h b/include/constant.h new file mode 100644 index 0000000..61de820 --- /dev/null +++ b/include/constant.h @@ -0,0 +1,10 @@ +#ifndef __constant__ +#define __constant__ + +#include "gsl/gsl_const_cgsm.h" +#include "gsl/gsl_const_mksa.h" +#include "gsl/gsl_const_num.h" + +#define GSL_CONST_NUM_PI (3.1415926535) + +#endif \ No newline at end of file diff --git a/include/gsl/gsl_const_cgsm.h b/include/gsl/gsl_const_cgsm.h new file mode 100644 index 0000000..2047e5e --- /dev/null +++ b/include/gsl/gsl_const_cgsm.h @@ -0,0 +1,122 @@ +/* const/gsl_const_cgsm.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, + * 2006, 2007, 2008, 2009 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#ifndef __GSL_CONST_CGSM__ +#define __GSL_CONST_CGSM__ + +#define GSL_CONST_CGSM_SPEED_OF_LIGHT (2.99792458e10) /* cm / s */ +#define GSL_CONST_CGSM_GRAVITATIONAL_CONSTANT (6.673e-8) /* cm^3 / g s^2 */ +#define GSL_CONST_CGSM_PLANCKS_CONSTANT_H (6.62606896e-27) /* g cm^2 / s */ +#define GSL_CONST_CGSM_PLANCKS_CONSTANT_HBAR (1.05457162825e-27) /* g cm^2 / s */ +#define GSL_CONST_CGSM_ASTRONOMICAL_UNIT (1.49597870691e13) /* cm */ +#define GSL_CONST_CGSM_LIGHT_YEAR (9.46053620707e17) /* cm */ +#define GSL_CONST_CGSM_PARSEC (3.08567758135e18) /* cm */ +#define GSL_CONST_CGSM_GRAV_ACCEL (9.80665e2) /* cm / s^2 */ +#define GSL_CONST_CGSM_ELECTRON_VOLT (1.602176487e-12) /* g cm^2 / s^2 */ +#define GSL_CONST_CGSM_MASS_ELECTRON (9.10938188e-28) /* g */ +#define GSL_CONST_CGSM_MASS_MUON (1.88353109e-25) /* g */ +#define GSL_CONST_CGSM_MASS_PROTON (1.67262158e-24) /* g */ +#define GSL_CONST_CGSM_MASS_NEUTRON (1.67492716e-24) /* g */ +#define GSL_CONST_CGSM_RYDBERG (2.17987196968e-11) /* g cm^2 / s^2 */ +#define GSL_CONST_CGSM_BOLTZMANN (1.3806504e-16) /* g cm^2 / K s^2 */ +#define GSL_CONST_CGSM_MOLAR_GAS (8.314472e7) /* g cm^2 / K mol s^2 */ +#define GSL_CONST_CGSM_STANDARD_GAS_VOLUME (2.2710981e4) /* cm^3 / mol */ +#define GSL_CONST_CGSM_MINUTE (6e1) /* s */ +#define GSL_CONST_CGSM_HOUR (3.6e3) /* s */ +#define GSL_CONST_CGSM_DAY (8.64e4) /* s */ +#define GSL_CONST_CGSM_WEEK (6.048e5) /* s */ +#define GSL_CONST_CGSM_INCH (2.54e0) /* cm */ +#define GSL_CONST_CGSM_FOOT (3.048e1) /* cm */ +#define GSL_CONST_CGSM_YARD (9.144e1) /* cm */ +#define GSL_CONST_CGSM_MILE (1.609344e5) /* cm */ +#define GSL_CONST_CGSM_NAUTICAL_MILE (1.852e5) /* cm */ +#define GSL_CONST_CGSM_FATHOM (1.8288e2) /* cm */ +#define GSL_CONST_CGSM_MIL (2.54e-3) /* cm */ +#define GSL_CONST_CGSM_POINT (3.52777777778e-2) /* cm */ +#define GSL_CONST_CGSM_TEXPOINT (3.51459803515e-2) /* cm */ +#define GSL_CONST_CGSM_MICRON (1e-4) /* cm */ +#define GSL_CONST_CGSM_ANGSTROM (1e-8) /* cm */ +#define GSL_CONST_CGSM_HECTARE (1e8) /* cm^2 */ +#define GSL_CONST_CGSM_ACRE (4.04685642241e7) /* cm^2 */ +#define GSL_CONST_CGSM_BARN (1e-24) /* cm^2 */ +#define GSL_CONST_CGSM_LITER (1e3) /* cm^3 */ +#define GSL_CONST_CGSM_US_GALLON (3.78541178402e3) /* cm^3 */ +#define GSL_CONST_CGSM_QUART (9.46352946004e2) /* cm^3 */ +#define GSL_CONST_CGSM_PINT (4.73176473002e2) /* cm^3 */ +#define GSL_CONST_CGSM_CUP (2.36588236501e2) /* cm^3 */ +#define GSL_CONST_CGSM_FLUID_OUNCE (2.95735295626e1) /* cm^3 */ +#define GSL_CONST_CGSM_TABLESPOON (1.47867647813e1) /* cm^3 */ +#define GSL_CONST_CGSM_TEASPOON (4.92892159375e0) /* cm^3 */ +#define GSL_CONST_CGSM_CANADIAN_GALLON (4.54609e3) /* cm^3 */ +#define GSL_CONST_CGSM_UK_GALLON (4.546092e3) /* cm^3 */ +#define GSL_CONST_CGSM_MILES_PER_HOUR (4.4704e1) /* cm / s */ +#define GSL_CONST_CGSM_KILOMETERS_PER_HOUR (2.77777777778e1) /* cm / s */ +#define GSL_CONST_CGSM_KNOT (5.14444444444e1) /* cm / s */ +#define GSL_CONST_CGSM_POUND_MASS (4.5359237e2) /* g */ +#define GSL_CONST_CGSM_OUNCE_MASS (2.8349523125e1) /* g */ +#define GSL_CONST_CGSM_TON (9.0718474e5) /* g */ +#define GSL_CONST_CGSM_METRIC_TON (1e6) /* g */ +#define GSL_CONST_CGSM_UK_TON (1.0160469088e6) /* g */ +#define GSL_CONST_CGSM_TROY_OUNCE (3.1103475e1) /* g */ +#define GSL_CONST_CGSM_CARAT (2e-1) /* g */ +#define GSL_CONST_CGSM_UNIFIED_ATOMIC_MASS (1.660538782e-24) /* g */ +#define GSL_CONST_CGSM_GRAM_FORCE (9.80665e2) /* cm g / s^2 */ +#define GSL_CONST_CGSM_POUND_FORCE (4.44822161526e5) /* cm g / s^2 */ +#define GSL_CONST_CGSM_KILOPOUND_FORCE (4.44822161526e8) /* cm g / s^2 */ +#define GSL_CONST_CGSM_POUNDAL (1.38255e4) /* cm g / s^2 */ +#define GSL_CONST_CGSM_CALORIE (4.1868e7) /* g cm^2 / s^2 */ +#define GSL_CONST_CGSM_BTU (1.05505585262e10) /* g cm^2 / s^2 */ +#define GSL_CONST_CGSM_THERM (1.05506e15) /* g cm^2 / s^2 */ +#define GSL_CONST_CGSM_HORSEPOWER (7.457e9) /* g cm^2 / s^3 */ +#define GSL_CONST_CGSM_BAR (1e6) /* g / cm s^2 */ +#define GSL_CONST_CGSM_STD_ATMOSPHERE (1.01325e6) /* g / cm s^2 */ +#define GSL_CONST_CGSM_TORR (1.33322368421e3) /* g / cm s^2 */ +#define GSL_CONST_CGSM_METER_OF_MERCURY (1.33322368421e6) /* g / cm s^2 */ +#define GSL_CONST_CGSM_INCH_OF_MERCURY (3.38638815789e4) /* g / cm s^2 */ +#define GSL_CONST_CGSM_INCH_OF_WATER (2.490889e3) /* g / cm s^2 */ +#define GSL_CONST_CGSM_PSI (6.89475729317e4) /* g / cm s^2 */ +#define GSL_CONST_CGSM_POISE (1e0) /* g / cm s */ +#define GSL_CONST_CGSM_STOKES (1e0) /* cm^2 / s */ +#define GSL_CONST_CGSM_STILB (1e0) /* cd / cm^2 */ +#define GSL_CONST_CGSM_LUMEN (1e0) /* cd sr */ +#define GSL_CONST_CGSM_LUX (1e-4) /* cd sr / cm^2 */ +#define GSL_CONST_CGSM_PHOT (1e0) /* cd sr / cm^2 */ +#define GSL_CONST_CGSM_FOOTCANDLE (1.076e-3) /* cd sr / cm^2 */ +#define GSL_CONST_CGSM_LAMBERT (1e0) /* cd sr / cm^2 */ +#define GSL_CONST_CGSM_FOOTLAMBERT (1.07639104e-3) /* cd sr / cm^2 */ +#define GSL_CONST_CGSM_CURIE (3.7e10) /* 1 / s */ +#define GSL_CONST_CGSM_ROENTGEN (2.58e-8) /* abamp s / g */ +#define GSL_CONST_CGSM_RAD (1e2) /* cm^2 / s^2 */ +#define GSL_CONST_CGSM_SOLAR_MASS (1.98892e33) /* g */ +#define GSL_CONST_CGSM_BOHR_RADIUS (5.291772083e-9) /* cm */ +#define GSL_CONST_CGSM_NEWTON (1e5) /* cm g / s^2 */ +#define GSL_CONST_CGSM_DYNE (1e0) /* cm g / s^2 */ +#define GSL_CONST_CGSM_JOULE (1e7) /* g cm^2 / s^2 */ +#define GSL_CONST_CGSM_ERG (1e0) /* g cm^2 / s^2 */ +#define GSL_CONST_CGSM_STEFAN_BOLTZMANN_CONSTANT (5.67040047374e-5) /* g / K^4 s^3 */ +#define GSL_CONST_CGSM_THOMSON_CROSS_SECTION (6.65245893699e-25) /* cm^2 */ +#define GSL_CONST_CGSM_BOHR_MAGNETON (9.27400899e-21) /* abamp cm^2 */ +#define GSL_CONST_CGSM_NUCLEAR_MAGNETON (5.05078317e-24) /* abamp cm^2 */ +#define GSL_CONST_CGSM_ELECTRON_MAGNETIC_MOMENT (9.28476362e-21) /* abamp cm^2 */ +#define GSL_CONST_CGSM_PROTON_MAGNETIC_MOMENT (1.410606633e-23) /* abamp cm^2 */ +#define GSL_CONST_CGSM_FARADAY (9.64853429775e3) /* abamp s / mol */ +#define GSL_CONST_CGSM_ELECTRON_CHARGE (1.602176487e-20) /* abamp s */ + +#endif /* __GSL_CONST_CGSM__ */ diff --git a/include/gsl/gsl_const_mksa.h b/include/gsl/gsl_const_mksa.h new file mode 100644 index 0000000..7e209e2 --- /dev/null +++ b/include/gsl/gsl_const_mksa.h @@ -0,0 +1,126 @@ +/* const/gsl_const_mksa.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, + * 2006, 2007, 2008, 2009 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#ifndef __GSL_CONST_MKSA__ +#define __GSL_CONST_MKSA__ + +#define GSL_CONST_MKSA_SPEED_OF_LIGHT (2.99792458e8) /* m / s */ +#define GSL_CONST_MKSA_GRAVITATIONAL_CONSTANT (6.67428e-11) /* m^3 / kg s^2 */ +#define GSL_CONST_MKSA_PLANCKS_CONSTANT_H (6.62606896e-34) /* kg m^2 / s */ +#define GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR (1.05457162825e-34) /* kg m^2 / s */ +#define GSL_CONST_MKSA_ASTRONOMICAL_UNIT (1.49597870691e11) /* m */ +#define GSL_CONST_MKSA_LIGHT_YEAR (9.46053620707e15) /* m */ +#define GSL_CONST_MKSA_PARSEC (3.08567758135e16) /* m */ +#define GSL_CONST_MKSA_GRAV_ACCEL (9.80665e0) /* m / s^2 */ +#define GSL_CONST_MKSA_ELECTRON_VOLT (1.602176487e-19) /* kg m^2 / s^2 */ +#define GSL_CONST_MKSA_MASS_ELECTRON (9.10938188e-31) /* kg */ +#define GSL_CONST_MKSA_MASS_MUON (1.88353109e-28) /* kg */ +#define GSL_CONST_MKSA_MASS_PROTON (1.67262158e-27) /* kg */ +#define GSL_CONST_MKSA_MASS_NEUTRON (1.67492716e-27) /* kg */ +#define GSL_CONST_MKSA_RYDBERG (2.17987196968e-18) /* kg m^2 / s^2 */ +#define GSL_CONST_MKSA_BOLTZMANN (1.3806504e-23) /* kg m^2 / K s^2 */ +#define GSL_CONST_MKSA_MOLAR_GAS (8.314472e0) /* kg m^2 / K mol s^2 */ +#define GSL_CONST_MKSA_STANDARD_GAS_VOLUME (2.2710981e-2) /* m^3 / mol */ +#define GSL_CONST_MKSA_MINUTE (6e1) /* s */ +#define GSL_CONST_MKSA_HOUR (3.6e3) /* s */ +#define GSL_CONST_MKSA_DAY (8.64e4) /* s */ +#define GSL_CONST_MKSA_WEEK (6.048e5) /* s */ +#define GSL_CONST_MKSA_INCH (2.54e-2) /* m */ +#define GSL_CONST_MKSA_FOOT (3.048e-1) /* m */ +#define GSL_CONST_MKSA_YARD (9.144e-1) /* m */ +#define GSL_CONST_MKSA_MILE (1.609344e3) /* m */ +#define GSL_CONST_MKSA_NAUTICAL_MILE (1.852e3) /* m */ +#define GSL_CONST_MKSA_FATHOM (1.8288e0) /* m */ +#define GSL_CONST_MKSA_MIL (2.54e-5) /* m */ +#define GSL_CONST_MKSA_POINT (3.52777777778e-4) /* m */ +#define GSL_CONST_MKSA_TEXPOINT (3.51459803515e-4) /* m */ +#define GSL_CONST_MKSA_MICRON (1e-6) /* m */ +#define GSL_CONST_MKSA_ANGSTROM (1e-10) /* m */ +#define GSL_CONST_MKSA_HECTARE (1e4) /* m^2 */ +#define GSL_CONST_MKSA_ACRE (4.04685642241e3) /* m^2 */ +#define GSL_CONST_MKSA_BARN (1e-28) /* m^2 */ +#define GSL_CONST_MKSA_LITER (1e-3) /* m^3 */ +#define GSL_CONST_MKSA_US_GALLON (3.78541178402e-3) /* m^3 */ +#define GSL_CONST_MKSA_QUART (9.46352946004e-4) /* m^3 */ +#define GSL_CONST_MKSA_PINT (4.73176473002e-4) /* m^3 */ +#define GSL_CONST_MKSA_CUP (2.36588236501e-4) /* m^3 */ +#define GSL_CONST_MKSA_FLUID_OUNCE (2.95735295626e-5) /* m^3 */ +#define GSL_CONST_MKSA_TABLESPOON (1.47867647813e-5) /* m^3 */ +#define GSL_CONST_MKSA_TEASPOON (4.92892159375e-6) /* m^3 */ +#define GSL_CONST_MKSA_CANADIAN_GALLON (4.54609e-3) /* m^3 */ +#define GSL_CONST_MKSA_UK_GALLON (4.546092e-3) /* m^3 */ +#define GSL_CONST_MKSA_MILES_PER_HOUR (4.4704e-1) /* m / s */ +#define GSL_CONST_MKSA_KILOMETERS_PER_HOUR (2.77777777778e-1) /* m / s */ +#define GSL_CONST_MKSA_KNOT (5.14444444444e-1) /* m / s */ +#define GSL_CONST_MKSA_POUND_MASS (4.5359237e-1) /* kg */ +#define GSL_CONST_MKSA_OUNCE_MASS (2.8349523125e-2) /* kg */ +#define GSL_CONST_MKSA_TON (9.0718474e2) /* kg */ +#define GSL_CONST_MKSA_METRIC_TON (1e3) /* kg */ +#define GSL_CONST_MKSA_UK_TON (1.0160469088e3) /* kg */ +#define GSL_CONST_MKSA_TROY_OUNCE (3.1103475e-2) /* kg */ +#define GSL_CONST_MKSA_CARAT (2e-4) /* kg */ +#define GSL_CONST_MKSA_UNIFIED_ATOMIC_MASS (1.660538782e-27) /* kg */ +#define GSL_CONST_MKSA_GRAM_FORCE (9.80665e-3) /* kg m / s^2 */ +#define GSL_CONST_MKSA_POUND_FORCE (4.44822161526e0) /* kg m / s^2 */ +#define GSL_CONST_MKSA_KILOPOUND_FORCE (4.44822161526e3) /* kg m / s^2 */ +#define GSL_CONST_MKSA_POUNDAL (1.38255e-1) /* kg m / s^2 */ +#define GSL_CONST_MKSA_CALORIE (4.1868e0) /* kg m^2 / s^2 */ +#define GSL_CONST_MKSA_BTU (1.05505585262e3) /* kg m^2 / s^2 */ +#define GSL_CONST_MKSA_THERM (1.05506e8) /* kg m^2 / s^2 */ +#define GSL_CONST_MKSA_HORSEPOWER (7.457e2) /* kg m^2 / s^3 */ +#define GSL_CONST_MKSA_BAR (1e5) /* kg / m s^2 */ +#define GSL_CONST_MKSA_STD_ATMOSPHERE (1.01325e5) /* kg / m s^2 */ +#define GSL_CONST_MKSA_TORR (1.33322368421e2) /* kg / m s^2 */ +#define GSL_CONST_MKSA_METER_OF_MERCURY (1.33322368421e5) /* kg / m s^2 */ +#define GSL_CONST_MKSA_INCH_OF_MERCURY (3.38638815789e3) /* kg / m s^2 */ +#define GSL_CONST_MKSA_INCH_OF_WATER (2.490889e2) /* kg / m s^2 */ +#define GSL_CONST_MKSA_PSI (6.89475729317e3) /* kg / m s^2 */ +#define GSL_CONST_MKSA_POISE (1e-1) /* kg m^-1 s^-1 */ +#define GSL_CONST_MKSA_STOKES (1e-4) /* m^2 / s */ +#define GSL_CONST_MKSA_STILB (1e4) /* cd / m^2 */ +#define GSL_CONST_MKSA_LUMEN (1e0) /* cd sr */ +#define GSL_CONST_MKSA_LUX (1e0) /* cd sr / m^2 */ +#define GSL_CONST_MKSA_PHOT (1e4) /* cd sr / m^2 */ +#define GSL_CONST_MKSA_FOOTCANDLE (1.076e1) /* cd sr / m^2 */ +#define GSL_CONST_MKSA_LAMBERT (1e4) /* cd sr / m^2 */ +#define GSL_CONST_MKSA_FOOTLAMBERT (1.07639104e1) /* cd sr / m^2 */ +#define GSL_CONST_MKSA_CURIE (3.7e10) /* 1 / s */ +#define GSL_CONST_MKSA_ROENTGEN (2.58e-4) /* A s / kg */ +#define GSL_CONST_MKSA_RAD (1e-2) /* m^2 / s^2 */ +#define GSL_CONST_MKSA_SOLAR_MASS (1.98892e30) /* kg */ +#define GSL_CONST_MKSA_BOHR_RADIUS (5.291772083e-11) /* m */ +#define GSL_CONST_MKSA_NEWTON (1e0) /* kg m / s^2 */ +#define GSL_CONST_MKSA_DYNE (1e-5) /* kg m / s^2 */ +#define GSL_CONST_MKSA_JOULE (1e0) /* kg m^2 / s^2 */ +#define GSL_CONST_MKSA_ERG (1e-7) /* kg m^2 / s^2 */ +#define GSL_CONST_MKSA_STEFAN_BOLTZMANN_CONSTANT (5.67040047374e-8) /* kg / K^4 s^3 */ +#define GSL_CONST_MKSA_THOMSON_CROSS_SECTION (6.65245893699e-29) /* m^2 */ +#define GSL_CONST_MKSA_BOHR_MAGNETON (9.27400899e-24) /* A m^2 */ +#define GSL_CONST_MKSA_NUCLEAR_MAGNETON (5.05078317e-27) /* A m^2 */ +#define GSL_CONST_MKSA_ELECTRON_MAGNETIC_MOMENT (9.28476362e-24) /* A m^2 */ +#define GSL_CONST_MKSA_PROTON_MAGNETIC_MOMENT (1.410606633e-26) /* A m^2 */ +#define GSL_CONST_MKSA_FARADAY (9.64853429775e4) /* A s / mol */ +#define GSL_CONST_MKSA_ELECTRON_CHARGE (1.602176487e-19) /* A s */ +#define GSL_CONST_MKSA_VACUUM_PERMITTIVITY (8.854187817e-12) /* A^2 s^4 / kg m^3 */ +#define GSL_CONST_MKSA_VACUUM_PERMEABILITY (1.25663706144e-6) /* kg m / A^2 s^2 */ +#define GSL_CONST_MKSA_DEBYE (3.33564095198e-30) /* A s^2 / m^2 */ +#define GSL_CONST_MKSA_GAUSS (1e-4) /* kg / A s^2 */ + +#endif /* __GSL_CONST_MKSA__ */ diff --git a/include/gsl/gsl_const_num.h b/include/gsl/gsl_const_num.h new file mode 100644 index 0000000..385a660 --- /dev/null +++ b/include/gsl/gsl_const_num.h @@ -0,0 +1,43 @@ +/* const/gsl_const_num.h + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, + * 2006, 2007, 2008, 2009 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#ifndef __GSL_CONST_NUM__ +#define __GSL_CONST_NUM__ + +#define GSL_CONST_NUM_FINE_STRUCTURE (7.297352533e-3) /* 1 */ +#define GSL_CONST_NUM_AVOGADRO (6.02214199e23) /* 1 / mol */ +#define GSL_CONST_NUM_YOTTA (1e24) /* 1 */ +#define GSL_CONST_NUM_ZETTA (1e21) /* 1 */ +#define GSL_CONST_NUM_EXA (1e18) /* 1 */ +#define GSL_CONST_NUM_PETA (1e15) /* 1 */ +#define GSL_CONST_NUM_TERA (1e12) /* 1 */ +#define GSL_CONST_NUM_GIGA (1e9) /* 1 */ +#define GSL_CONST_NUM_MEGA (1e6) /* 1 */ +#define GSL_CONST_NUM_KILO (1e3) /* 1 */ +#define GSL_CONST_NUM_MILLI (1e-3) /* 1 */ +#define GSL_CONST_NUM_MICRO (1e-6) /* 1 */ +#define GSL_CONST_NUM_NANO (1e-9) /* 1 */ +#define GSL_CONST_NUM_PICO (1e-12) /* 1 */ +#define GSL_CONST_NUM_FEMTO (1e-15) /* 1 */ +#define GSL_CONST_NUM_ATTO (1e-18) /* 1 */ +#define GSL_CONST_NUM_ZEPTO (1e-21) /* 1 */ +#define GSL_CONST_NUM_YOCTO (1e-24) /* 1 */ + +#endif /* __GSL_CONST_NUM__ */ diff --git a/include/property.h b/include/property.h new file mode 100644 index 0000000..90cd4dc --- /dev/null +++ b/include/property.h @@ -0,0 +1,58 @@ +#ifndef __property__ +#define __property__ + +#define READ_ONLY 1 +#define WRITE_ONLY 2 +#define READ_WRITE 3 + +template +class property { +public: + property() { + m_cObject = NULL; + Set = NULL; + Get = NULL; + } + //-- This to set a pointer to the class that contain the + // property -- + void setContainer(Container* cObject) { m_cObject = cObject; } + //-- Set the set member function that will change the value -- + void setter(void (Container::*pSet)(ValueType value)) { + if ((nPropType == WRITE_ONLY) || (nPropType == READ_WRITE)) + Set = pSet; + else + Set = NULL; + } + //-- Set the get member function that will retrieve the value -- + void getter(ValueType (Container::*pGet)()) { + if ((nPropType == READ_ONLY) || (nPropType == READ_WRITE)) + Get = pGet; + else + Get = NULL; + } + //-- Overload the '=' sign to set the value using the set + // member -- + ValueType operator=(const ValueType& value) { + assert(m_cObject != NULL); + assert(Set != NULL); + (m_cObject->*Set)(value); + return value; + } + //-- To make possible to cast the property class to the + // internal type -- + operator ValueType() { + assert(m_cObject != NULL); + assert(Get != NULL); + return (m_cObject->*Get)(); + } + +private: + Container* m_cObject; //-- Pointer to the module that + // contains the property -- + void (Container::*Set)(ValueType value); + //-- Pointer to set member function -- + ValueType (Container::*Get)(); + //-- Pointer to get member function -- +}; + +#endif \ No newline at end of file diff --git a/include/utils.h b/include/utils.h new file mode 100644 index 0000000..8371839 --- /dev/null +++ b/include/utils.h @@ -0,0 +1,43 @@ +#ifndef __utils__ +#define __utils__ + +#include +#include + +template +std::vector cutArrs(std::vector &Arrs, int begin, int end) { + std::vector result; + result.assign(Arrs.begin() + begin, Arrs.begin() + end); + return result; +} + +template +double average(T x, int len) { + double sum = 0; + for (int i = 0; i < len; i++) sum += x[i]; + return sum / len; +} + +template +double variance(T x, int len) { + double sum = 0; + double avg = average(x, len); + for (int i = 0; i < len; i++) sum += pow(x[i] - avg, 2); + return sum / len; +} + +template +double standardDev(T x, int len) { + double var = variance(x, len); + return sqrt(var); +} + +template +double maxBias(T x, int len) { + double res = -1; + double avg = average(x, len); + for (int i = 0; i < len; i++) res = std::max(res, abs(avg - x[i])); + return res; +} + +#endif \ No newline at end of file diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..ad75f90 --- /dev/null +++ b/main.cpp @@ -0,0 +1,286 @@ +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +using namespace std; +using Eigen::Vector2d; +namespace plt = matplotlibcpp; + +#define PI GSL_CONST_NUM_PI +#define kb GSL_CONST_MKSA_BOLTZMANN + +#define sigma 3.405e-10 +#define epsilon 1.65e-21 +#define mass (6.63382e-26 / 39.95) +#define tem (epsilon / kb) +#define vel sqrt(epsilon / mass) +#define tim (sigma * sqrt(mass / epsilon)) +#define M 39.95 + +bool debug; +int N, W, Epoch; +double L, T, dt, eps; + +class Molecule { +public: + Molecule(){}; + Molecule(double m, Vector2d r0, Vector2d v0) { + r = r0; + v = v0; + cnt = 0; + length = 0; + lastCrash = 0; + this->m = m; + }; + ~Molecule(){}; + +public: + int cnt, lastCrash; + double m, length; + Vector2d a, aPre, r, v; + +public: + void warp() { + r[0] -= L * floor(r[0] / L); + r[1] -= L * floor(r[1] / L); + } + void crash(int index) { + cnt += 1; + length += v.norm() * (index - lastCrash) * dt; + lastCrash = index; + } + void stretch(double factor) { v *= factor; } +}; + +Vector2d separateDistance(Vector2d ra, Vector2d rb) { + double dx = ra[0] - rb[0]; + double dy = ra[1] - rb[1]; + if (dx > L / 2) + dx -= L; + else if (dx < -L / 2) + dx += L; + if (dy > L / 2) + dy -= L; + else if (dy < -L / 2) + dy += L; + return Vector2d(dx, dy); +} + +double lennardJones(Vector2d r) { return 4 * (pow(1. / r.norm(), 12) - pow(1. / r.norm(), 6)); } + +Vector2d lennardJonesFource(Vector2d r) { + double fm = 48 * (pow(1. / r.norm(), 14) - 0.5 * pow(1. / r.norm(), 8)); + return fm * r; +} + +double getKinetic(vector mol) { + double Ek = 0; + for (int i = 0; i < N; i++) Ek += 0.5 * mol[i].m * mol[i].v.squaredNorm(); + return Ek; +} + +double getPotential(vector mol) { + double Vr = 0; + for (int i = 0; i < N; i++) + for (int j = i + 1; j < N; j++) Vr += lennardJones(separateDistance(mol[i].r, mol[j].r)); + return Vr; +} + +double getEnergy(vector mol) { + double Ek = getKinetic(mol); + double Vr = getPotential(mol); + return Ek + Vr; +} + +double getTemperature(vector mol) { + double Ek = getKinetic(mol); + return Ek / (N - 1); +} + +double getPressure(vector mol) { + double Te = getTemperature(mol); + double p = N * Te; + + Vector2d r, f; + for (int i = 1; i < N; i++) + for (int j = 0; j < i; j++) { + r = separateDistance(mol[i].r, mol[j].r); + f = lennardJonesFource(r); + p += r.dot(f) / 2; + } + return p / pow(L, 2); +} + +int main(int argc, char const *argv[]) { + auto args = util::argparser("Argon MD Simulate"); + args.set_program_name("ComSim") + .add_help_option() + .add_argument("length", "length of box in sigma") + .add_argument("number", "number of particles") + .add_argument("temperature", "temperature in Kelvin") + .add_option("-D", "--debug", "debug mode") + .add_option("-E", "--epoch", "maximum number of rounds simulated,\ndefault is 10000", 10000) + .add_option("-W", "--width", "width of window to get average,\ndefault is 1000", 1000) + .add_option("", "--dt", "delta t to simualte,\ndefault is 0.01", 0.01) + .add_option("", "--eps", "error t0 end simualte,\ndefault is 0.01", 0.01) + .parse(argc, argv); + + L = args.get_argument_double("length"); + N = args.get_argument_int("number"); + T = args.get_argument_double("temperature") / tem; + debug = args.get_option_bool("-D"); + Epoch = args.get_option_int("-E"); + W = min(args.get_option_int("-W"), Epoch); + dt = args.get_option_double("--dt"); + eps = args.get_option_double("--eps"); + + int cnt = 0; + double fac, eqT, eqE, eqP; + vector argon(N); + Vector2d Vc = Vector2d(0, 0); + vector x(N), y(N); + vector Ti, Te, MTe, En, MEn, Pr, MPr; + + default_random_engine seed(time(NULL)); + uniform_real_distribution u01(0, 1); + normal_distribution gauss(0, sqrt(T * tem * kb / (M * mass))); + + int NM = sqrt(N); + double H = (L - 1.) / (NM - 1); + for (int i = 0; i < N; i++) { + double vx = gauss(seed) / vel; + double vy = gauss(seed) / vel; + + argon[i] = Molecule(M, Vector2d((i / NM) * H + 0.5, (i % NM) * H + 0.5), Vector2d(vx, vy)); + } + + for (int i = 0; i < N; i++) Vc += argon[i].v; + Vc /= N; + for (int i = 0; i < N; i++) argon[i].v -= Vc; + eqT = getTemperature(argon); + + do { + fac = sqrt(T / eqT); + for (int i = 0; i < N; i++) { + argon[i].stretch(fac); + argon[i].a = Vector2d(0, 0); + for (int j = i + 1; j < N; j++) { + Vector2d f = lennardJonesFource(separateDistance(argon[i].r, argon[j].r)); + argon[i].a += f; + argon[j].a -= f; + } + argon[i].a /= argon[i].m; + } + + Ti.push_back(Epoch * cnt * dt); + Te.push_back(getTemperature(argon)); + En.push_back(getEnergy(argon)); + Pr.push_back(getPressure(argon)); + MTe.push_back(Te[Epoch * cnt]); + MEn.push_back(En[Epoch * cnt]); + MPr.push_back(Pr[Epoch * cnt]); + + for (int i = 1; i < Epoch; i++) { + for (int j = 0; j < N; j++) { + argon[j].r += argon[j].v * dt + argon[j].a * pow(dt, 2) / 2; + argon[j].aPre = argon[j].a; + argon[j].a = Vector2d(0, 0); + } + for (int j = 0; j < N; j++) { + for (int k = j + 1; k < N; k++) { + Vector2d f = lennardJonesFource(separateDistance(argon[j].r, argon[k].r)); + argon[j].a += f; + argon[k].a -= f; + } + argon[j].a /= argon[j].m; + argon[j].v += (argon[j].aPre + argon[j].a) * dt / 2; + argon[j].warp(); + } + + int k = i + Epoch * cnt; + Ti.push_back(k * dt); + Te.push_back(getTemperature(argon)); + En.push_back(getEnergy(argon)); + Pr.push_back(getPressure(argon)); + + if (i >= W) { + MTe.push_back((MTe[k - 1] * W + Te[k] - Te[k - W]) / W); + MEn.push_back((MEn[k - 1] * W + En[k] - En[k - W]) / W); + MPr.push_back((MPr[k - 1] * W + Pr[k] - Pr[k - W]) / W); + } else { + MTe.push_back((MTe[k - 1] * i + Te[k]) / (i + 1)); + MEn.push_back((MEn[k - 1] * i + En[k]) / (i + 1)); + MPr.push_back((MPr[k - 1] * i + Pr[k]) / (i + 1)); + } + } + + cnt += 1; + eqT = average(cutArrs(MTe, Epoch * cnt - W, Epoch * cnt - 1), W); + eqE = average(cutArrs(MEn, Epoch * cnt - W, Epoch * cnt - 1), W); + eqP = average(cutArrs(MPr, Epoch * cnt - W, Epoch * cnt - 1), W); + if (debug) cout << eqT * tem << ", " << abs(eqT - T) / T << ", " << eqE << ", " << eqP << endl; + + if (abs(eqT - T) / T < eps) break; + } while (1); + cout << eqT * tem << ", " << abs(eqT - T) / T << ", " << eqE << ", " << eqP << endl; + + plt::figure_size(900, 600); + plt::subplot(3, 1, 1); + plt::named_plot("Instant", Ti, Te); + plt::named_plot("Average", Ti, MTe); + plt::legend(); + plt::title("Equilibrium Temperature : " + to_string(eqT)); + plt::xlabel("Time / ($\\sigma\\sqrt{m/\\epsilon}$)"); + plt::ylabel("Temperature / ($\\epsilon/k_b$)"); + + plt::subplot(3, 1, 2); + plt::named_plot("Instant", Ti, En); + plt::named_plot("Average", Ti, MEn); + plt::legend(); + plt::title("Equilibrium Energy : " + to_string(eqE)); + plt::xlabel("Time / ($\\sigma\\sqrt{m/\\epsilon}$)"); + plt::ylabel("Energy / ($\\epsilon$)"); + + plt::subplot(3, 1, 3); + plt::named_plot("Instant", Ti, Pr); + plt::named_plot("Average", Ti, MPr); + plt::legend(); + plt::title("Equilibrium Pressure : " + to_string(eqP)); + plt::xlabel("Time / ($\\sigma\\sqrt{m/\\epsilon}$)"); + plt::ylabel("Pressure / ($\\epsilon/\\sigma^2$)"); + + plt::suptitle("Argon N = " + to_string(N)); + plt::tight_layout(); + plt::save("result.png"); + plt::show(); + + return 0; +} + +/* + Unit of quantities: + Length : σ (Lennard-Jones parameters) + Energy : ε (Lennard-Jones parameters) + Masses : m (mass of the atom) + Velocity : sqrt(ε / m) + Time : σ * sqrt(m / ε) + Temperature : ε / k (k is Boltzmann's constant) +*/ + +/* + 2D Maxwell Distribution + normal_distribution gauss(0, sqrt(T * kb / m)); + vector x(100000); + for (int i = 0; i < 100000; i++) x[i] = sqrt(pow(gauss(seed), 2) + pow(gauss(seed), 2)); + plt::hist(x, 100); + plt::show(); +*/