From c24116bed3b87d1599ac763d271c6a894df63011 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 10:44:29 -0500 Subject: [PATCH 01/13] Changed FOOEXE to SWIFTEST_DRIVER --- CMakeLists.txt | 2 +- src/CMakeLists.txt | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 34b05b21e..255faac9c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -54,7 +54,7 @@ ENDIF(FCNAME STREQUAL "pgf90") ############################################################ # Define the executable name -SET(FOOEXE swiftest_driver) +SET(SWIFTEST_DRIVER swiftest_driver) # Define some directories SET(SRC ${CMAKE_SOURCE_DIR}/src) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 25411a4dc..b4204732c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -104,26 +104,26 @@ SET(FOO_src ) # Define the executable in terms of the source files -ADD_EXECUTABLE(${FOOEXE} ${FOO_src}) +ADD_EXECUTABLE(${SWIFTEST_DRIVER} ${FOO_src}) ##################################################### # Add the needed libraries and special compiler flags ##################################################### # Uncomment if you need to link to BLAS and LAPACK -TARGET_LINK_LIBRARIES(${FOOEXE} ${NETCDF_LIBRARIES} ${NETCDF_FORTRAN_LIBRARIES}) +TARGET_LINK_LIBRARIES(${SWIFTEST_DRIVER} ${NETCDF_LIBRARIES} ${NETCDF_FORTRAN_LIBRARIES}) # Uncomment if you have parallization IF(USE_OPENMP) - SET_TARGET_PROPERTIES(${FOOEXE} PROPERTIES + SET_TARGET_PROPERTIES(${SWIFTEST_DRIVER} PROPERTIES COMPILE_FLAGS "${OpenMP_Fortran_FLAGS}" LINK_FLAGS "${OpenMP_Fortran_FLAGS}") ELSEIF(USE_MPI) - SET_TARGET_PROPERTIES(${FOOEXE} PROPERTIES + SET_TARGET_PROPERTIES(${SWIFTEST_DRIVER} PROPERTIES COMPILE_FLAGS "${MPI_Fortran_COMPILE_FLAGS}" LINK_FLAGS "${MPI_Fortran_LINK_FLAGS}") INCLUDE_DIRECTORIES(${MPI_Fortran_INCLUDE_PATH}) - TARGET_LINK_LIBRARIES(${FOOEXE} ${MPI_Fortran_LIBRARIES}) + TARGET_LINK_LIBRARIES(${SWIFTEST_DRIVER} ${MPI_Fortran_LIBRARIES}) ENDIF(USE_OPENMP) @@ -137,4 +137,4 @@ IF(WIN32) ELSE() SET(CMAKE_INSTALL_PREFIX /usr/local) ENDIF(WIN32) -INSTALL(TARGETS ${FOOEXE} RUNTIME DESTINATION bin) \ No newline at end of file +INSTALL(TARGETS ${SWIFTEST_DRIVER} RUNTIME DESTINATION bin) \ No newline at end of file From e6319df7091ae9a19cba17b7f9d501fe5014601b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 11:01:40 -0500 Subject: [PATCH 02/13] Fixed parallelization flag generation so code can compile with OpenMP --- CMakeLists.txt | 3 +-- cmake/Modules/FindOpenMP_Fortran.cmake | 10 ++++------ 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 255faac9c..3125fc3e2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -30,7 +30,7 @@ ENDIF(NOT CMAKE_Fortran_COMPILER_SUPPORTS_F90) # Set some options the user may choose # Uncomment the below if you want the user to choose a parallelization library OPTION(USE_MPI "Use the MPI library for parallelization" OFF) -OPTION(USE_OPENMP "Use OpenMP for parallelization" OFF) +OPTION(USE_OPENMP "Use OpenMP for parallelization" ON) # This INCLUDE statement executes code that sets the compile flags for DEBUG, # RELEASE, and TESTING. You should review this file and make sure the flags @@ -40,7 +40,6 @@ INCLUDE(${CMAKE_MODULE_PATH}/SetFortranFlags.cmake) # taken care of here, such as the fact that the FindOpenMP routine doesn't know # about Fortran. INCLUDE(${CMAKE_MODULE_PATH}/SetParallelizationLibrary.cmake) - INCLUDE(${CMAKE_MODULE_PATH}/SetUpNetCDF.cmake) # There is an error in CMAKE with this flag for pgf90. Unset it diff --git a/cmake/Modules/FindOpenMP_Fortran.cmake b/cmake/Modules/FindOpenMP_Fortran.cmake index c8e0ca2b4..32777569e 100644 --- a/cmake/Modules/FindOpenMP_Fortran.cmake +++ b/cmake/Modules/FindOpenMP_Fortran.cmake @@ -26,14 +26,14 @@ INCLUDE (${CMAKE_ROOT}/Modules/FindPackageHandleStandardArgs.cmake) SET (OpenMP_Fortran_FLAG_CANDIDATES - #Microsoft Visual Studio - "/openmp" - #Intel windows - "/Qopenmp" #Intel "-qopenmp" + #Intel windows + "/Qopenmp" #Gnu "-fopenmp" + #Portland Group + "-mp" #Empty, if compiler automatically accepts openmp " " #Sun @@ -42,8 +42,6 @@ SET (OpenMP_Fortran_FLAG_CANDIDATES "+Oopenmp" #IBM XL C/c++ "-qsmp" - #Portland Group - "-mp" ) IF (DEFINED OpenMP_Fortran_FLAGS) From 33b77d1d6a0ee2fde1526e317d162acb9c7b4e75 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 16:00:56 -0500 Subject: [PATCH 03/13] Fixed CMake flags so that the fast or strict math model is usedon specific source files where necessary. Added a PROFILE build type that turns on the optimization report in ifort. --- CMakeLists.txt | 4 +- README.md | 6 ++- cmake/Modules/SetCompileFlag.cmake | 4 +- cmake/Modules/SetFortranFlags.cmake | 62 +++++++++++++++++------------ src/CMakeLists.txt | 33 ++++++++++----- 5 files changed, 67 insertions(+), 42 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 3125fc3e2..96ce5efa9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,7 @@ # You should have received a copy of the GNU General Public License along with Swiftest. # If not, see: https://www.gnu.org/licenses. -# CMake project file for FOO +# CMake project file for SWIFTEST ################################################## # Define the project and the depencies that it has @@ -63,7 +63,7 @@ SET(BIN ${CMAKE_SOURCE_DIR}/bin) # Have the .mod files placed in the lib folder SET(CMAKE_Fortran_MODULE_DIRECTORY ${LIB}) -# The source for the FOO binary and have it placed in the bin folder +# The source for the SWIFTEST binary and have it placed in the bin folder ADD_SUBDIRECTORY(${SRC} ${BIN}) # Add a distclean target to the Makefile diff --git a/README.md b/README.md index 74c1a652c..5e935af08 100644 --- a/README.md +++ b/README.md @@ -88,10 +88,14 @@ To build Swiftest with the release flags (default), type the following: ``` $ cmake .. ``` -To buid with the debug flags, type: +To build with the debug flags, type: ``` $ cmake .. -DCMAKE_BUILD_TYPE=DEBUG ``` +To build with profiling flags, type: +``` +$ cmake .. -DCMAKE_BUILD_TYPE=PROFILE +``` Finally, to build with the testing flags, type: ``` $ cmake .. -DCMAKE_BUILD_TYPE=TESTING diff --git a/cmake/Modules/SetCompileFlag.cmake b/cmake/Modules/SetCompileFlag.cmake index 1d110ae6d..4141c4773 100644 --- a/cmake/Modules/SetCompileFlag.cmake +++ b/cmake/Modules/SetCompileFlag.cmake @@ -23,12 +23,12 @@ # "-Wall" # GNU # "-warn all" # Intel # ) -# The optin "-Wall" will be checked first, and if it works, will be +# The option "-Wall" will be checked first, and if it works, will be # appended to the CMAKE_C_FLAGS variable. If it doesn't work, then # "-warn all" will be tried. If this doesn't work then checking will # terminate because REQUIRED was given. # -# The reasong that the variable must be given twice (first as the name then +# The reasoning that the variable must be given twice (first as the name then # as the value in quotes) is because of the way CMAKE handles the passing # of variables in functions; it is difficult to extract a variable's # contents and assign new values to it from within a function. diff --git a/cmake/Modules/SetFortranFlags.cmake b/cmake/Modules/SetFortranFlags.cmake index e0b21862b..7850fbdb8 100644 --- a/cmake/Modules/SetFortranFlags.cmake +++ b/cmake/Modules/SetFortranFlags.cmake @@ -21,32 +21,36 @@ STRING(TOUPPER "${CMAKE_BUILD_TYPE}" BT) IF(BT STREQUAL "RELEASE") SET(CMAKE_BUILD_TYPE RELEASE CACHE STRING - "Choose the type of build, options are DEBUG, RELEASE, or TESTING." + "Choose the type of build, options are DEBUG, RELEASE, PROFILE, or TESTING." FORCE) ELSEIF(BT STREQUAL "DEBUG") SET (CMAKE_BUILD_TYPE DEBUG CACHE STRING - "Choose the type of build, options are DEBUG, RELEASE, or TESTING." + "Choose the type of build, options are DEBUG, RELEASE, PROFILE, or TESTING." FORCE) ELSEIF(BT STREQUAL "TESTING") SET (CMAKE_BUILD_TYPE TESTING CACHE STRING - "Choose the type of build, options are DEBUG, RELEASE, or TESTING." + "Choose the type of build, options are DEBUG, RELEASE, PROFILE, or TESTING." FORCE) +ELSEIF(BT STREQUAL "PROFILE") + SET (CMAKE_BUILD_TYPE PROFILE CACHE STRING + "Choose the type of build, options are DEBUG, RELEASE, PROFILE, or TESTING." + FORCE) ELSEIF(NOT BT) SET(CMAKE_BUILD_TYPE RELEASE CACHE STRING - "Choose the type of build, options are DEBUG, RELEASE, or TESTING." + "Choose the type of build, options are DEBUG, RELEASE, PROFILE, or TESTING." FORCE) MESSAGE(STATUS "CMAKE_BUILD_TYPE not given, defaulting to RELEASE") ELSE() - MESSAGE(FATAL_ERROR "CMAKE_BUILD_TYPE not valid, choices are DEBUG, RELEASE, or TESTING") + MESSAGE(FATAL_ERROR "CMAKE_BUILD_TYPE not valid, choices are DEBUG, RELEASE, PROFILE, or TESTING") ENDIF(BT STREQUAL "RELEASE") ######################################################### # If the compiler flags have already been set, return now ######################################################### -IF(CMAKE_Fortran_FLAGS_RELEASE AND CMAKE_Fortran_FLAGS_TESTING AND CMAKE_Fortran_FLAGS_DEBUG) +IF(CMAKE_Fortran_FLAGS_RELEASE AND CMAKE_Fortran_FLAGS_TESTING AND CMAKE_Fortran_FLAGS_DEBUG AND CMAKE_Fortran_FLAGS_PROFILE) RETURN () -ENDIF(CMAKE_Fortran_FLAGS_RELEASE AND CMAKE_Fortran_FLAGS_TESTING AND CMAKE_Fortran_FLAGS_DEBUG) +ENDIF(CMAKE_Fortran_FLAGS_RELEASE AND CMAKE_Fortran_FLAGS_TESTING AND CMAKE_Fortran_FLAGS_DEBUG AND CMAKE_Fortran_FLAGS_PROFILE) ######################################################################## # Determine the appropriate flags for this compiler for each build type. @@ -81,7 +85,6 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS}" ################### ### DEBUG FLAGS ### ################### - # NOTE: debugging symbols (-g or /debug:full) are already on by default # Disable optimizations @@ -163,7 +166,7 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}" # Aligns a variable to a specified boundary and offset SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG}" - Fortran "-align all" # Intel + Fortran "-align all -align array64byte" # Intel ) # Enables changing the variable and array memory layout @@ -215,7 +218,6 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_TESTING "${CMAKE_Fortran_FLAGS_TESTING}" ##################### ### RELEASE FLAGS ### ##################### - # NOTE: agressive optimizations (-O3) are already turned on by default # Unroll loops @@ -234,19 +236,6 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" "-Minline" # Portland Group ) -# Interprocedural (link-time) optimizations -SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" - Fortran "-ipo" # Intel - "/Qipo" # Intel Windows - "-flto" # GNU - "-Mipa" # Portland Group - ) - -# Single-file optimizations -SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" - Fortran "-ip" # Intel - "/Qip" # Intel Windows - ) # Allows for lines longer than 80 characters without truncation SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" @@ -299,7 +288,28 @@ SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" Fortran "-fma" # Intel ) -# Enables agressive optimixation on floating-points -SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" - Fortran "-fp-model=fast" # Intel + +##################### +### MATH FLAGS ### +##################### +# Some subroutines require more strict floating point operation optimizations for repeatability +SET_COMPILE_FLAG(STRICTMATH_FLAGS "${STRICTMATH_FLAGS}" + Fortran "-fp-model=precise -prec-div -prec-sqrt -assume protect-parens" # Intel + "/fp:precise /Qprec-div /Qprec-sqrt /assume:protect-parens" # Intel Windows + ) + +# Most subroutines can use aggressive optimization of floating point operations without problems. +SET_COMPILE_FLAG(FASTMATH_FLAGS "${FASTMATH_FLAGS}" + Fortran "-fp-model=fast" + "/fp:fast" + ) + +##################### +### PROFILE FLAGS ### +##################### +# Enables the optimization reports to be generated +SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_PROFILE "${CMAKE_Fortran_FLAGS_RELEASE}" + Fortran "-pg -qopt-report=5 -traceback -p -g3" # Intel + "/Qopt-report:5 /traceback -g3" # Windows Intel + "-pg -fbacktrace" ) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b4204732c..d467332f8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -12,7 +12,7 @@ ######################################## # Add the source files -SET(FOO_src +SET(FAST_MATH_FILES ${SRC}/modules/encounter_classes.f90 ${SRC}/modules/fraggle_classes.f90 ${SRC}/modules/helio_classes.f90 @@ -41,12 +41,10 @@ SET(FOO_src ${SRC}/gr/gr.f90 ${SRC}/helio/helio_drift.f90 ${SRC}/helio/helio_gr.f90 - ${SRC}/helio/helio_kick.f90 ${SRC}/helio/helio_setup.f90 ${SRC}/helio/helio_step.f90 ${SRC}/helio/helio_util.f90 ${SRC}/io/io.f90 - ${SRC}/kick/kick.f90 ${SRC}/netcdf/netcdf.f90 ${SRC}/obl/obl.f90 ${SRC}/operators/operator_cross.f90 @@ -55,7 +53,6 @@ SET(FOO_src ${SRC}/rmvs/rmvs_discard.f90 ${SRC}/rmvs/rmvs_encounter_check.f90 ${SRC}/rmvs/rmvs_io.f90 - ${SRC}/rmvs/rmvs_kick.f90 ${SRC}/rmvs/rmvs_setup.f90 ${SRC}/rmvs/rmvs_step.f90 ${SRC}/rmvs/rmvs_util.f90 @@ -66,7 +63,6 @@ SET(FOO_src ${SRC}/symba/symba_encounter_check.f90 ${SRC}/symba/symba_gr.f90 ${SRC}/symba/symba_io.f90 - ${SRC}/symba/symba_kick.f90 ${SRC}/symba/symba_setup.f90 ${SRC}/symba/symba_step.f90 ${SRC}/symba/symba_util.f90 @@ -96,24 +92,31 @@ SET(FOO_src ${SRC}/whm/whm_coord.f90 ${SRC}/whm/whm_drift.f90 ${SRC}/whm/whm_gr.f90 - ${SRC}/whm/whm_kick.f90 ${SRC}/whm/whm_setup.f90 ${SRC}/whm/whm_step.f90 ${SRC}/whm/whm_util.f90 ${SRC}/main/swiftest_driver.f90 ) +SET(STRICT_MATH_FILES + ${SRC}/kick/kick.f90 + ${SRC}/helio/helio_kick.f90 + ${SRC}/rmvs/rmvs_kick.f90 + ${SRC}/symba/symba_kick.f90 + ${SRC}/whm/whm_kick.f90 +) + +set(SWIFTEST_src ${FAST_MATH_FILES} ${STRICT_MATH_FILES}) # Define the executable in terms of the source files -ADD_EXECUTABLE(${SWIFTEST_DRIVER} ${FOO_src}) +ADD_EXECUTABLE(${SWIFTEST_DRIVER} ${SWIFTEST_src}) ##################################################### # Add the needed libraries and special compiler flags ##################################################### -# Uncomment if you need to link to BLAS and LAPACK +# # Uncomment if you need to link to BLAS and LAPACK TARGET_LINK_LIBRARIES(${SWIFTEST_DRIVER} ${NETCDF_LIBRARIES} ${NETCDF_FORTRAN_LIBRARIES}) -# Uncomment if you have parallization IF(USE_OPENMP) SET_TARGET_PROPERTIES(${SWIFTEST_DRIVER} PROPERTIES COMPILE_FLAGS "${OpenMP_Fortran_FLAGS}" @@ -127,7 +130,6 @@ ELSEIF(USE_MPI) ENDIF(USE_OPENMP) - ##################################### # Tell how to install this executable ##################################### @@ -137,4 +139,13 @@ IF(WIN32) ELSE() SET(CMAKE_INSTALL_PREFIX /usr/local) ENDIF(WIN32) -INSTALL(TARGETS ${SWIFTEST_DRIVER} RUNTIME DESTINATION bin) \ No newline at end of file +INSTALL(TARGETS ${SWIFTEST_DRIVER} RUNTIME DESTINATION bin) + + +#Set strict vs fast math flags +STRING(TOUPPER "${CMAKE_BUILD_TYPE}" BT) +IF(BT STREQUAL "RELEASE" OR BT STREQUAL "PROFILE") + SET_PROPERTY(SOURCE ${STRICT_MATH_FILES} APPEND_STRING PROPERTY COMPILE_FLAGS "${STRICTMATH_FLAGS}") + SET_PROPERTY(SOURCE ${FAST_MATH_FILES} APPEND_STRING PROPERTY COMPILE_FLAGS "${FASTMATH_FLAGS}") +ENDIF() + From 167c09a33744ee0f72a5d37b7da9672ca842ec81 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 16:01:36 -0500 Subject: [PATCH 04/13] Removed simd call as it isn't compatible with the function anyway --- src/helio/helio_drift.f90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index 60c6d52a8..1076532c0 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -114,11 +114,9 @@ subroutine helio_drift_linear_all(xh, pt, dt, n, lmask) ! Internals integer(I4B) :: i - !$omp parallel do simd default(shared) schedule(static) do i = 1, n if (lmask(i)) call helio_drift_linear_one(xh(1,i), xh(2,i), xh(3,i), pt(1), pt(2), pt(3), dt) end do - !$omp end parallel do simd return end subroutine helio_drift_linear_all From 40d58b3a7d41ba967e452a63ba56d3399f2d9bed Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 17:02:32 -0500 Subject: [PATCH 05/13] Fixed problem with the reading of the oldorigin info --- src/netcdf/netcdf.f90 | 59 ++++++++++++++++++++----------------------- 1 file changed, 28 insertions(+), 31 deletions(-) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 4c037a291..16333dec8 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -406,7 +406,6 @@ module subroutine netcdf_open(self, param, readonly) end if end if - end if if ((param%out_form == EL) .or. (param%out_form == XVEL)) then @@ -437,29 +436,22 @@ module subroutine netcdf_open(self, param, readonly) ! end if ! Optional Variables - - status = nf90_inq_varid(self%ncid, NPL_VARNAME, self%npl_varid) - if (status /= nf90_noerr) write(*,*) "Warning! NPL variable not set in input file. Calculating." - - status = nf90_inq_varid(self%ncid, NTP_VARNAME, self%ntp_varid) - if (status /= nf90_noerr) write(*,*) "Warning! NTP variable not set in input file. Calculating." - - if (param%integrator == SYMBA) then - status = nf90_inq_varid(self%ncid, NPLM_VARNAME, self%nplm_varid) - if (status /= nf90_noerr) write(*,*) "Warning! NPLM variable not set in input file. Calculating." - end if - if (param%lrhill_present) then status = nf90_inq_varid(self%ncid, RHILL_VARNAME, self%rhill_varid) if (status /= nf90_noerr) write(*,*) "Warning! RHILL variable not set in input file. Calculating." end if - ! Variables The User Doesn't Need to Know About - + ! Optional variables The User Doesn't Need to Know About + status = nf90_inq_varid(self%ncid, NPL_VARNAME, self%npl_varid) + status = nf90_inq_varid(self%ncid, NTP_VARNAME, self%ntp_varid) status = nf90_inq_varid(self%ncid, STATUS_VARNAME, self%status_varid) status = nf90_inq_varid(self%ncid, J2RP2_VARNAME, self%j2rp2_varid) status = nf90_inq_varid(self%ncid, J4RP4_VARNAME, self%j4rp4_varid) + if (param%integrator == SYMBA) then + status = nf90_inq_varid(self%ncid, NPLM_VARNAME, self%nplm_varid) + end if + if (param%lclose) then status = nf90_inq_varid(self%ncid, ORIGIN_TYPE_VARNAME, self%origin_type_varid) status = nf90_inq_varid(self%ncid, ORIGIN_TIME_VARNAME, self%origin_time_varid) @@ -985,7 +977,7 @@ module subroutine netcdf_read_particle_info_system(self, iu, param, plmask, tpma if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_time_varid, rtemp), "netcdf_read_particle_info_system nf90_getvar origin_time_varid" ) else - rtemp = 0.0_DP + rtemp = param%t0 end if call cb%info%set_value(origin_time=rtemp(1)) @@ -996,29 +988,31 @@ module subroutine netcdf_read_particle_info_system(self, iu, param, plmask, tpma call tp%info(i)%set_value(origin_time=rtemp(tpind(i))) end do - status = nf90_inq_varid(iu%ncid, ORIGIN_XHX_VARNAME, iu%origin_xhx_varid) if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_xhx_varid, rtemp_arr(1,:)), "netcdf_read_particle_info_system nf90_getvar origin_xhx_varid" ) - else - ! [TODO]: This doesn't work when the input mode is EL. This needs to get filled in later after xv2el has been called - ! call check( nf90_get_var(iu%ncid, iu%xhx_varid, rtemp_arr(1,:)), "netcdf_read_particle_info_system nf90_getvar xhx_varid" ) + else if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + call check( nf90_get_var(iu%ncid, iu%xhx_varid, rtemp_arr(1,:)), "netcdf_read_particle_info_system nf90_getvar xhx_varid" ) + else + rtemp_arr(1,:) = 0._DP end if status = nf90_inq_varid(iu%ncid, ORIGIN_XHY_VARNAME, iu%origin_xhy_varid) if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_xhy_varid, rtemp_arr(2,:)), "netcdf_read_particle_info_system nf90_getvar origin_xhy_varid" ) - else - ! [TODO]: This doesn't work when the input mode is EL. This needs to get filled in later after xv2el has been called - ! call check( nf90_get_var(iu%ncid, iu%xhy_varid, rtemp_arr(2,:)), "netcdf_read_particle_info_system nf90_getvar xhy_varid" ) + else if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + call check( nf90_get_var(iu%ncid, iu%xhy_varid, rtemp_arr(2,:)), "netcdf_read_particle_info_system nf90_getvar xhx_varid" ) + else + rtemp_arr(2,:) = 0._DP end if status = nf90_inq_varid(iu%ncid, ORIGIN_XHZ_VARNAME, iu%origin_xhz_varid) if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_xhz_varid, rtemp_arr(3,:)), "netcdf_read_particle_info_system nf90_getvar origin_xhz_varid" ) + else if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + call check( nf90_get_var(iu%ncid, iu%xhz_varid, rtemp_arr(3,:)), "netcdf_read_particle_info_system nf90_getvar xhz_varid" ) else - ! [TODO]: This doesn't work when the input mode is EL. This needs to get filled in later after xv2el has been called - ! call check( nf90_get_var(iu%ncid, iu%xhz_varid, rtemp_arr(3,:)), "netcdf_read_particle_info_system nf90_getvar xhz_varid" ) + rtemp_arr(3,:) = 0._DP end if do i = 1, npl @@ -1031,25 +1025,28 @@ module subroutine netcdf_read_particle_info_system(self, iu, param, plmask, tpma status = nf90_inq_varid(iu%ncid, ORIGIN_VHX_VARNAME, iu%origin_vhx_varid) if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_vhx_varid, rtemp_arr(1,:)), "netcdf_read_particle_info_system nf90_getvar origin_vhx_varid" ) + else if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + call check( nf90_get_var(iu%ncid, iu%vhx_varid, rtemp_arr(1,:)), "netcdf_read_particle_info_system nf90_getvar vhx_varid" ) else - ! [TODO]: This doesn't work when the input mode is EL. This needs to get filled in later after xv2el has been called - ! call check( nf90_get_var(iu%ncid, iu%vhx_varid, rtemp_arr(1,:)), "netcdf_read_particle_info_system nf90_getvar vhx_varid" ) + rtemp_arr(1,:) = 0._DP end if status = nf90_inq_varid(iu%ncid, ORIGIN_VHY_VARNAME, iu%origin_vhy_varid) if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_vhy_varid, rtemp_arr(2,:)), "netcdf_read_particle_info_system nf90_getvar origin_vhy_varid" ) + else if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + call check( nf90_get_var(iu%ncid, iu%vhy_varid, rtemp_arr(2,:)), "netcdf_read_particle_info_system nf90_getvar vhy_varid" ) else - ! [TODO]: This doesn't work when the input mode is EL. This needs to get filled in later after xv2el has been called - ! call check( nf90_get_var(iu%ncid, iu%vhy_varid, rtemp_arr(2,:)), "netcdf_read_particle_info_system nf90_getvar vhy_varid" ) + rtemp_arr(2,:) = 0._DP end if status = nf90_inq_varid(iu%ncid, ORIGIN_VHZ_VARNAME, iu%origin_vhz_varid) if (status == nf90_noerr) then call check( nf90_get_var(iu%ncid, iu%origin_vhz_varid, rtemp_arr(3,:)), "netcdf_read_particle_info_system nf90_getvar origin_vhz_varid" ) + else if ((param%out_form == XV) .or. (param%out_form == XVEL)) then + call check( nf90_get_var(iu%ncid, iu%vhz_varid, rtemp_arr(3,:)), "netcdf_read_particle_info_system nf90_getvar vhz_varid" ) else - ! [TODO]: This doesn't work when the input mode is EL. This needs to get filled in later after xv2el has been called - ! call check( nf90_get_var(iu%ncid, iu%vhz_varid, rtemp_arr(3,:)), "netcdf_read_particle_info_system nf90_getvar vhz_varid" ) + rtemp_arr(3,:) = 0._DP end if do i = 1, npl From 75b910691c4302b0f46859528fefc3b6bbbf3cd6 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 17:14:14 -0500 Subject: [PATCH 06/13] Added OMP_NUM_THREADS environment variable to the Jupyter notebook in the Basic Simulation example and cleared all its output --- .../Basic_Simulation/run_simulation.ipynb | 110 ++---------------- 1 file changed, 9 insertions(+), 101 deletions(-) diff --git a/examples/Basic_Simulation/run_simulation.ipynb b/examples/Basic_Simulation/run_simulation.ipynb index 8f2ab51b3..ea78c2690 100644 --- a/examples/Basic_Simulation/run_simulation.ipynb +++ b/examples/Basic_Simulation/run_simulation.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "86c845ce-1801-46ca-8a8a-1cabb266e6a6", "metadata": {}, "outputs": [], @@ -10,122 +10,30 @@ "import swiftest\n", "import xarray as xr\n", "import numpy as np\n", - "import os" + "import os\n", + "%env OMP_NUM_THREADS=8" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "d716c371-8eb4-4fc1-82af-8b5c444c831e", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Reading Swiftest file /home/daminton/git_debug/swiftest/examples/Basic_Simulation/param.in\n" - ] - } - ], + "outputs": [], "source": [ "sim = swiftest.Simulation()" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "id": "ec7452d6-4c9b-4df3-acc0-b11c32264b91", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "tstop 10.0 y\n", - "Writing parameter inputs to file /home/daminton/git_debug/swiftest/examples/Basic_Simulation/param.in\n", - "Running a Swiftest symba run from tstart=0.0 y to tstop=10.0 y\n", - "\u001b]2;cd /home/daminton/git_debug/swiftest/examples/Basic_Simulation\u0007\u001b]1;\u0007\u001b]2;/home/daminton/git_debug/swiftest/bin/swiftest_driver symba \u0007\u001b]1;\u0007 Parameter input file is /home/daminton/git_debug/swiftest/examples/Basic_Simulation/param.in\n", - " \n", - " Warning! NPLM variable not set in input file. Calculating.\n", - " *************** Main Loop *************** \n", - "Time = 1.00000E+00; fraction done = 0.100; Number of active plm, pl, tp = 13, 14, 10\n", - "Integration steps: Total wall time: 6.30684E-01; Interval wall time: 5.29890E-01;Interval wall time/step: 2.70141E-03\n", - "Time = 2.00000E+00; fraction done = 0.200; Number of active plm, pl, tp = 13, 14, 10\n", - "Integration steps: Total wall time: 1.66455E+00; Interval wall time: 5.27720E-01;Interval wall time/step: 2.99766E-03\n", - "Time = 3.00000E+00; fraction done = 0.300; Number of active plm, pl, tp = 13, 14, 10\n", - "Integration steps: Total wall time: 2.64051E+00; Interval wall time: 5.20805E-01;Interval wall time/step: 2.88832E-03\n", - "Time = 4.00000E+00; fraction done = 0.400; Number of active plm, pl, tp = 13, 14, 10\n", - "Integration steps: Total wall time: 3.60585E+00; Interval wall time: 5.24579E-01;Interval wall time/step: 2.82311E-03\n", - "Time = 5.00000E+00; fraction done = 0.500; Number of active plm, pl, tp = 13, 14, 10\n", - "Integration steps: Total wall time: 4.58823E+00; Interval wall time: 5.37595E-01;Interval wall time/step: 2.96439E-03\n", - "Time = 6.00000E+00; fraction done = 0.600; Number of active plm, pl, tp = 13, 14, 10\n", - "Integration steps: Total wall time: 5.55600E+00; Interval wall time: 5.27663E-01;Interval wall time/step: 2.83397E-03\n", - "Time = 7.00000E+00; fraction done = 0.700; Number of active plm, pl, tp = 13, 14, 10\n", - "Integration steps: Total wall time: 6.64194E+00; Interval wall time: 5.83431E-01;Interval wall time/step: 3.11589E-03\n", - "Time = 8.00000E+00; fraction done = 0.800; Number of active plm, pl, tp = 13, 14, 10\n", - "Integration steps: Total wall time: 7.63044E+00; Interval wall time: 5.52408E-01;Interval wall time/step: 2.95843E-03\n", - "Time = 9.00000E+00; fraction done = 0.900; Number of active plm, pl, tp = 13, 14, 10\n", - "Integration steps: Total wall time: 8.62339E+00; Interval wall time: 5.46944E-01;Interval wall time/step: 3.01578E-03\n", - "Time = 1.00000E+01; fraction done = 1.000; Number of active plm, pl, tp = 13, 14, 10\n", - "Integration steps: Total wall time: 9.67297E+00; Interval wall time: 6.00750E-01;Interval wall time/step: 3.20243E-03\n", - "\n", - "Normal termination of Swiftest (version 1.0)\n", - "------------------------------------------------\n", - "\n", - "Creating Dataset from NetCDF file\n", - "Successfully converted 11 output frames.\n", - "Swiftest simulation data stored as xarray DataSet .ds\n" - ] - } - ], + "outputs": [], "source": [ "sim.run(tstop=10.0)" ] }, - { - "cell_type": "code", - "execution_count": 9, - "id": "5b0a57a6-dbd5-4d34-8e6e-91fc8ad8789f", - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[,\n", - " ,\n", - " ,\n", - " ,\n", - " ,\n", - " ,\n", - " ,\n", - " ,\n", - " ,\n", - " ,\n", - " ,\n", - " ,\n", - " ,\n", - " ]" - ] - }, - "execution_count": 9, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGwCAYAAABcnuQpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAACEU0lEQVR4nOzdfXyN9f/A8de12c7uN5vZDWez2TDGRiKjkHszFaEIaygJCSWiKFnubyIhG2VuSlTfiuKbkbvITS18VWyNH2uRtsy2s51z/f6YnXZsY9M456z389H16Fyfz+f6XO9zOTvX+3yuO0VVVRUhhBBCCCtlY+4AhBBCCCH+CUlmhBBCCGHVJJkRQgghhFWTZEYIIYQQVk2SGSGEEEJYNUlmhBBCCGHVJJkRQgghhFWrYe4A7jSDwcCFCxdwdXVFURRzhyOEEEKIClBVlb/++gt/f39sbG4+9lLtk5kLFy6g1WrNHYYQQgghbsO5c+eoW7fuTdtU+2TG1dUVKNoYbm5uZo5GCCGEEBWRnZ2NVqs17sdvptonM8WHltzc3CSZEUIIIaxMRU4RkROAhRBCCGHVJJkRQgghhFWTZEYIIYQQVq3anzMjhCXT6/UUFBSYOwzxL2dnZ4etra25wxDitkkyI4QZqKpKRkYGf/75p7lDEQIADw8PfH195X5cwipZTDITHx/PlClTeO6551i0aBFQ9IU/Y8YMVq5cyZUrV2jdujXLli2jSZMm5g1WiH+oOJGpXbs2Tk5OsgMRZqOqKteuXSMzMxMAPz8/M0ckROVZRDJz+PBhVq5cSbNmzUzK58yZw4IFC1izZg0NGjRg5syZdOnShdOnT1founMhLJFerzcmMl5eXuYORwgcHR0ByMzMpHbt2nLISVgds58AfPXqVQYNGsSqVauoWbOmsVxVVRYtWsTLL79Mnz59CA8PZ+3atVy7do3169ebMWIh/pnic2ScnJzMHIkQfyv+PMo5XMIamT2ZefbZZ4mOjqZz584m5ampqWRkZNC1a1djmUajoX379uzfv7/c/vLz88nOzjaZhLBEcmhJWBL5PAprZtbDTBs3buTo0aMcPny4VF1GRgYAPj4+JuU+Pj78+uuv5fYZHx/PjBkzqjZQIYQQQlgss43MnDt3jueee45169bh4OBQbrsbfy2oqnrTXxCTJ08mKyvLOJ07d67KYhZCCCGE5THbyMyRI0fIzMzknnvuMZbp9Xr27NnD0qVLOX36NFA0QlPy7PrMzMxSozUlaTQaNBrNnQtcCCGEEBbFbMlMp06dSElJMSl78sknadSoEZMmTSI4OBhfX1927NhB8+bNAdDpdOzevZvZs2ebI+S7TlVVoGhSVcP1/xeXGa6/NlxvazBthwrX64teG0r0WXLESwFFQUEpen3DfFG7vyeT5cqsL66zMamT4/FCCCHuFLMlM66uroSHh5uUOTs74+XlZSwfN24cs2bNIjQ0lNDQUGbNmoWTkxMDBw40R8gmzv/fetLTV/2dXKjXk4brSYRKyeTi72Tj7ySkYklK9XOrJKh43ub66+LEyxZFsUFRbP+esAXFBkWpUaLu+jw2ULJtmXUlluXvdrdfZ2MSm2JTAxvFDkWpgY2N3fVl7CgsrIHBYIten4deT6ntIUmgEEJUjkXcZ6Y8L774Irm5uYwaNcp407yvvvrKIu4xU1iQTW5uurnDuAXTxODvnaRaNGhDyaRKLTF/JxWv74bSO71aC2Jj44+H+6vk5iro9RVIVm4yclZUXV5d6XnTtpRfd5ME62Z9mvZx8ytkOnToQLNmzXBwcODdd9/F3t6ekSNHMn36dAAWLFhAYmIiZ8+exdPTk5iYGObMmYOLiwsAa9asYdy4caxbt44JEyZw7tw5evbsydq1a9m8eTOvvvoqWVlZPPHEEyxatMh47xSdTsfUqVNJSkrizz//JDw8nNmzZ9OhQ4db/1sIISySRSUzycnJJvOKojB9+nTjl5sl8fV9CI+a917/Ure54dDM9VEFxebvL31FKRoRKPGF//foQ4lDMhVtV2aSUlxv849/0RePDpkmObdOgsoadeL6mBUl64yHwUovZ7psyVj0f08YUNXConYly6/XoeoxqPrrdSWW5e92xjr0qIbiuuL+Cv+uU8tZF6Z9c0NcqlqIwVBwfb7A+Bq1ZolRI+WG7VHmP0bxFrw+f+dTzqpTfiKk1+eydm0io0fHsWvXhxw6dIynn36Re+4JoVOnB9Drs5g7dyqBgVrSfj3PuOemMGHCNRYvfhNFUSgoyObatWssWjSf999/h7/+usqAAbE8/HAMHh4efPLJJlJT0xgwYAht2txD//6PAgpDh8bx66/pJCWtxd/fn48//oTu3bvz/fdHCA1tUHbcyKXLwrxUVb3+XVLy+67w7+8cgx6uf78Z1MKiowVltS81X047KrY8qgGDWoiXZztq1XrQbNtHUdXq/Zs4Ozsbd3d3srKycHNzM3c4QpCXl0dqaipBQUGlruQzTQbLSxpvnlSWbltWEslN6m623K36rLjo6Dj0ej3bt681lnXsOJAHHmjFjBnjSrXfuvUrxo+fSWrqHgCSkj5h1KhpHDv2OcHBWgDGjXudTZv+w88/J+PiUnQTuD59RhIQUIdFi6Zx9uw5WrToxalTO/Dzq23su3fvEdxzTzivvvrcTSK+cTSrqOzGpOfGkbNSo2XGMkr0U8Zr47rKqTe+LNnq9tvk5RWQnn4ed/ffsLPXXz8ca4Ni8iPq+g+M4nTa+DksLisuL/EZK1V2wzIl2t6s//KWMUn0jespbmu4vlP/+4eJ8QcJaukfNer19iV+qFD8g8XYV/GPlus/Yq73c2Nb47yxraHED56S9fq/6411NyYZhSW2g2UKDHyGkPoTq7TPyuy/LWpkRoh/u793GmBtAwEVTqiu74hsbBwID2+Eo2Ogsc7fX8uVK3k4OPizK/kb5sxeyP/+9xPZ2X9RWKgnLy+PggJHnJwdsbV1wsnJkYYNw43Jna+vH4GBWtw9vI07PB8fHy5d+hMbGw0pKT+jqir33NObkjuH/PwCPL08KLnDLuMdXv/PtL66/BzU6VQKCq7wy5l4DIYL5g5HVIKi2JU4l6/keYI1riejNVBsbAFbbIrLjecA3jjd2McNfdnU+Pu8wBLtPTxamXUbSDIjhKgSlU3EFMUWjcYZO7u/f3HZ2toDdly8eJWHH3qckSNH8sYbs/H09GTv3r0MGzYMW9taODp4YG/viZ2dPY6OAcbl7e1rYm/vhLNTsLHMzs4NGxsDLi4NsLc/hq2tLUeOHC31/CEXFxfc3HxRbxgpUEv80i9rvmSCdvP5G/oprjZJjtQb/l88EnFjfRltbtrfrdep1xdiY6PB3a0FKIGYXiFZdMi56Dz9v0eliv+vlBhxMh35uT66pJTR1mSZstuW3X/JuoosU+LiAEqerF80X7q+eETq74sOjCf+F49WGXfm1+uxKdG2ZL1StCx/Xxzw92hX8bpuqEcpkVAUX9xQRpJS4sIGIcmMEMICfffddxQWFjJ//nxsbIq+rD/44IN/3G/z5s3R6/VkZmZy//33l9lGuWHHaG0jZLfLxiYPe/sCgoJm3/RGpkJYIknphBAWp379+hQWFvLWW29x9uxZ3n//fd55551/3G+DBg0YNGgQQ4YMYcuWLaSmpnL48GFmz57NF198UQWRCyHMQZIZIYTFiYyMZMGCBcyePZvw8HCSkpKIj4+vkr4TExMZMmQIEyZMoGHDhvTu3Ztvv/0WrVZbJf0LIe4+uZpJiLvsZlczCWEu8rkUlqYy+28ZmRFCCCGEVZNkRgghhBBWTZIZIYQQQlg1SWaEEEIIYdUkmRFCCCGEVZNkRgghhBBWTZIZIYQQQlg1SWaEEEIIYdUkmRFCCCGEVZNkRghRYbGxsSiKwsiRI0vVjRo1CkVRiI2NvfuBCSH+1SSZEUJUilarZePGjeTm5hrL8vLy2LBhAwEBAbfdr6qqFBYWVkWIJvR6PQaDocr7FUJYDklmhBCV0qJFCwICAtiyZYuxbMuWLWi1Wpo3b24sU1WVOXPmEBwcjKOjIxEREWzevNlYn5ycjKIofPnll7Rs2RKNRsM333yDwWBg9uzZhISEoNFoCAgI4I033jBZ5s8//zT2c/z4cRRFIS0tDYA1a9bg4eHBZ599RuPGjY392tnZkZGRYfJeJkyYwAMPPHAHtpIQ4m6qYe4AhBBFO/7cAr1Z1u1oZ4uiKJVa5sknnyQxMZFBgwYBkJCQQFxcHMnJycY2U6dOZcuWLSxfvpzQ0FD27NnDE088gbe3N+3btze2e/HFF5k3bx7BwcF4eHgwefJkVq1axcKFC2nXrh0XL17kf//7X6Xiu3btGvHx8bz77rt4eXlRt25dgoODef/993nhhRcAKCwsZN26dbz55puV6lsIYXkkmRHCAuQW6Gn8ypdmWffJ17rhZF+5r4LBgwczefJk0tLSUBSFffv2sXHjRmMyk5OTw4IFC/j6669p06YNAMHBwezdu5cVK1aYJDOvvfYaXbp0AeCvv/5i8eLFLF26lKFDhwJQv3592rVrV6n4CgoKePvtt4mIiDCWDRs2jMTERGMy8/nnn3Pt2jX69+9fqb6FEJZHkhkhRKXVqlWL6Oho1q5di6qqREdHU6tWLWP9yZMnycvLMyYpxXQ6ncmhKICWLVsaX586dYr8/Hw6der0j+Kzt7enWbNmJmWxsbFMnTqVgwcPct9995GQkED//v1xdnb+R+sSQpifJDNCWABHO1tOvtbNbOu+HXFxcYwePRqAZcuWmdQVn3D7+eefU6dOHZM6jUZjMl8ymXB0dLzpOm1sik7zU1XVWFZQUFCqnaOjY6lDZ7Vr1yYmJobExESCg4P54osvTA6LCSGslyQzQlgARVEqfajH3Lp3745OpwOgWzfTRKz4xNv09HSTQ0q3EhoaiqOjI//9738ZPnx4qXpvb28ALl68SM2aNYGiE4Aravjw4Tz22GPUrVuX+vXr07Zt2wovK4SwXNb17SmEsBi2tracOnXK+LokV1dXJk6cyPPPP4/BYKBdu3ZkZ2ezf/9+XFxcjOfD3MjBwYFJkybx4osvYm9vT9u2bfn99985ceIEw4YNIyQkBK1Wy/Tp05k5cyY///wz8+fPr3DM3bp1w93dnZkzZ/Laa6/d/psXQlgUSWaEELfNzc2t3LrXX3+d2rVrEx8fz9mzZ/Hw8KBFixZMmTLlpn1OmzaNGjVq8Morr3DhwgX8/PyMN+mzs7Njw4YNPPPMM0RERHDvvfcyc+ZM+vXrV6F4bWxsiI2NZdasWQwZMqTib1QIYdEUteTB52ooOzsbd3d3srKybvrFK8TdkpeXR2pqKkFBQTg4OJg7nH+dESNG8Ntvv/Hpp5+aOxSLIp9LYWkqs/+WkRkhxL9CVlYWhw8fJikpiU8++cTc4QghqpAkM0KIf4WHHnqIQ4cO8fTTT5e6ZFwIYd0kmRFC/CvIZdhCVF/ybCYhhBBCWDVJZoQQQghh1cyazCxfvpxmzZrh5uaGm5sbbdq0Ydu2bcb62NhYFEUxme677z4zRiyEEEIIS2PWc2bq1q3Lm2++SUhICABr167loYce4tixYzRp0gQoustoYmKicRl7e3uzxCqEEEIIy2TWZCYmJsZk/o033mD58uUcPHjQmMxoNBp8fX3NEZ4QQgghrIDFnDOj1+vZuHEjOTk5tGnTxlienJxM7dq1adCgASNGjCAzM/Om/eTn55OdnW0yCSGEEKL6Mnsyk5KSgouLCxqNhpEjR7J161YaN24MQI8ePUhKSuLrr79m/vz5HD58mAcffJD8/Pxy+4uPj8fd3d04abXau/VWhBBCCGEGZn+cgU6nIz09nT///JOPPvqId999l927dxsTmpIuXrxIYGAgGzdupE+fPmX2l5+fb5LsZGdno9Vq5XEGwmJY623jY2JiyM3NZefOnaXqDhw4QFRUFEeOHKFFixZmiE78U9b6uRTVl1U9zsDe3t54AnDLli05fPgwixcvZsWKFaXa+vn5ERgYyM8//1xufxqNBo1Gc8fiFeLfatiwYfTp04dff/2VwMBAk7qEhAQiIyMlkRFCmIXZDzPdSFXVcg8jXb58mXPnzuHn53eXoxJC9OrVi9q1a7NmzRqT8mvXrrFp0yaGDRvG/v37eeCBB3B0dESr1TJ27FhycnKMbevVq8esWbOIi4vD1dWVgIAAVq5caaxPTk5GURT+/PNPY9nx48dRFIW0tDQAfv31V2JiYqhZsybOzs40adKEL7744k6+dSGEhTNrMjNlyhS++eYb0tLSSElJ4eWXXyY5OZlBgwZx9epVJk6cyIEDB0hLSyM5OZmYmBhq1arFI488Ys6whah6qgq6HPNMFTzSXKNGDYYMGcKaNWsoeXT6ww8/RKfTERERQbdu3ejTpw8//PADmzZtYu/evYwePdqkn/nz59OyZUuOHTvGqFGjeOaZZ/jf//5X4U317LPPkp+fz549e0hJSWH27Nm4uLhUeHkhRPVj1sNMv/32G4MHD+bixYu4u7vTrFkztm/fTpcuXcjNzSUlJYX33nuPP//8Ez8/Pzp27MimTZtwdXU1Z9hCVL2CazDL3zzrnnIB7J0r1DQuLo65c+eSnJxMx44dgaJDTH369GHVqlUMHDiQcePGARAaGsqSJUto3749y5cvN56H0bNnT0aNGgXApEmTWLhwIcnJyTRq1KhCMaSnp9O3b1+aNm0KQHBwcGXerRCiGjJrMrN69epy6xwdHfnyyy/vYjRCiFtp1KgRUVFRJCQk0LFjR86cOcM333zDV199xXPPPccvv/xCUlKSsb2qqhgMBlJTUwkLCwOgWbNmxnpFUfD19b3lLRdKGjt2LM888wxfffUVnTt3pm/fviZ9CiH+fcx+ArAQArBzKhohMde6K2HYsGGMHj2aZcuWkZiYSGBgIJ06dcJgMPD0008zduzYUssEBAT8vTo7O5M6RVEwGAwA2NgUHfkueRiroKDApP3w4cPp1q0bn3/+OV999RXx8fHMnz+fMWPGVOp9CCGqD0lmhLAEilLhQz3m1r9/f5577jnWr1/P2rVrGTFiBIqi0KJFC06cOGG8OvF2eHt7A0W3YahZsyZQdALwjbRaLSNHjmTkyJFMnjyZVatWSTIjxL+YxV3NJISwbC4uLgwYMIApU6Zw4cIFYmNjgaLzXw4cOMCzzz7L8ePH+fnnn/n0008rlWSEhISg1WqZPn06P/30E59//jnz5883aTNu3Di+/PJLUlNTOXr0KF9//bXxEJYQ4t9JkhkhRKUNGzaMK1eu0LlzZ+MhpGbNmrF7925+/vln7r//fpo3b860adMqdSsFOzs7NmzYwP/+9z8iIiKYPXs2M2fONGmj1+t59tlnCQsLo3v37jRs2JC33367St+fEMK6mP0OwHdaZe4gKMTdIHdaFZZIPpfC0lRm/y0jM0IIIYSwapLMCCGEEMKqSTIjhBBCCKsmyYwQQgghrJokM0IIIYSwapLMCCGEEMKqSTIjhBBCCKsmyYwQQgghrJokM0IIIYSwapLMCCEsUlpaGoqilPmgSSGEKEmSGSFEhcXGxqIoSqmpe/fu/7jfhx9+uGqCFEL869QwdwBCCOvSvXt3EhMTTco0Gs1t9aXX61EUpSrCEkL8i8nIjBCiUjQaDb6+viZTzZo1AViwYAFNmzbF2dkZrVbLqFGjuHr1qnHZNWvW4OHhwWeffUbjxo3RaDQ8+eSTrF27lk8++cQ40pOcnGxc5uzZs3Ts2BEnJyciIiI4cODA3X7LQggLJyMzQlgAVVXJLcw1y7odazhW2eiIjY0NS5YsoV69eqSmpjJq1ChefPFF3n77bWOba9euER8fz7vvvouXlxe+vr7k5eWRnZ1tHPHx9PTkwoULALz88svMmzeP0NBQXn75ZR5//HF++eUXatSQry8hRBH5NhDCAuQW5tJ6fWuzrPvbgd/iZOdU4fafffYZLi4uJmWTJk1i2rRpjBs3zlgWFBTE66+/zjPPPGOSzBQUFPD2228TERFhLHN0dCQ/Px9fX99S65s4cSLR0dEAzJgxgyZNmvDLL7/QqFGjCscshKjeJJkRQlRKx44dWb58uUmZp6cnALt27WLWrFmcPHmS7OxsCgsLycvLIycnB2dnZwDs7e1p1qxZhddXsq2fnx8AmZmZkswIIYwkmRHCAjjWcOTbgd+abd2V4ezsTEhISKnyX3/9lZ49ezJy5Ehef/11PD092bt3L8OGDaOgoODv9TlW7rCWnZ2d8XXxcgaDoVIxCyGqN0lmhLAAiqJU6lCPJfruu+8oLCxk/vz52NgUXVvwwQcfVGhZe3t79Hr9nQxPCFGNSTIjhKiU/Px8MjIyTMpq1KhB/fr1KSws5K233iImJoZ9+/bxzjvvVKjPevXq8eWXX3L69Gm8vLxwd3e/E6ELIaopuTRbCFEp27dvx8/Pz2Rq164dkZGRLFiwgNmzZxMeHk5SUhLx8fEV6nPEiBE0bNiQli1b4u3tzb59++7wuxBCVCeKqqqquYO4k7Kzs3F3dycrKws3NzdzhyMEeXl5pKamEhQUhIODg7nDEQKQz6WwPJXZf8vIjBBCCCGsmiQzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqZk1mli9fTrNmzXBzc8PNzY02bdqwbds2Y72qqkyfPh1/f38cHR3p0KEDJ06cMGPEQgghhLA0Zr0DcN26dXnzzTeNz3lZu3YtDz30EMeOHaNJkybMmTOHBQsWsGbNGho0aMDMmTPp0qULp0+fxtXV1ZyhWy1VVUtN5ZVbUtuS8Zf3viypvJiiKMbnCRW/Ll4mLy/PZPnynld0s+cYVXaZipbfbL4yz1USQoi7weJumufp6cncuXOJi4vD39+fcePGMWnSJKDoNuo+Pj7Mnj2bp59+ukL93amb5qWlpXH27FlUVcVgMJj835LLhPm5uLjQtm1b6tSpQ40a1v1EkZslOVVdVxWvzdXWGhJAuWmesDSV2X9bzDepXq/nww8/JCcnhzZt2pCamkpGRgZdu3Y1ttFoNLRv3579+/eXm8zk5+eTn59vnM/Ozr4j8Z47d449e/bckb4tTVmjC7eaqrJtcX3JeG4Wa2XK7+Qy5Y04aTQaatSogZ2dXbnJzM1+Y9ytUaWKqMiomTD1T5Kyf7r8zV7rdDry8vI4cuSI8W/PxsYGGxsb4+uK/k1b0nQ729Aaks877cZR8eKysl5D0TaztbW9K7GVxezJTEpKCm3atCEvLw8XFxe2bt1K48aN2b9/PwA+Pj4m7X18fPj111/L7S8+Pp4ZM2bc0ZgB/P39adWqlckf/Y1/8JZYduNruHXSIapW8S/gmjVrWuQv4JJJ2I3lcXFxvPfee4wYMYK3337bpG7MmDGsWLGCwYMH8+6775bbR1mvb1Z3O+3Kel3Rstvtq6IsNQEsLCwkLy+PAwcOcPXqVXOHY3FuJ6m83XbFbvbZv3G+qupuV7t27ejcuXOV9HU7zJ7MNGzYkOPHj/Pnn3/y0UcfMXToUHbv3m2sv/EfWFXVm+5gJ0+ezPjx443z2dnZaLXaKo+7fv361K9fv8r7FcLcbvbrVFEUtFotH3zwAYsXL8bR0REoStA2bdpEQEAANjY22Nvb39a6CwoKsLOzu/3gzayiSdKdSNb+6fI6nQ57e3saNmxIQUGByeHpkoetLWEqGc/dcqvEQpiX2ZMZe3t74wnALVu25PDhwyxevNh4nkxGRgZ+fn7G9pmZmaVGa0rSaDRoNJo7G7QQ/2ItWrTg7NmzbNmyhUGDBgGwZcsWtFotwcHBxnbbt29n5syZ/Pjjj9ja2tKmTRsWL15s/BGQlpZGUFAQmzZt4u233+bgwYMsX76cBx98kNGjR7N37150Oh316tVj7ty59OzZ0yzvtzKs+TBFXl4eTk5OdOnSxSJHDMtzq8SnZLsbl7tbdZXto6KjO7eav5ttzXmICSwgmbmRqqrk5+cTFBSEr68vO3bsoHnz5kDRL4fdu3cze/ZsM0cpRNVSVRU1N9cs61YcHSu9833yySdJTEw0JjMJCQnExcWRnJxsbJOTk8P48eNp2rQpOTk5vPLKKzzyyCMcP34cG5u/7woxadIk5s+fT2JiIhqNhqeeegqdTseePXtwdnbm5MmTuLi4VMl7FdWPHA4XYOZkZsqUKfTo0QOtVstff/3Fxo0bSU5OZvv27SiKwrhx45g1axahoaGEhoYya9YsnJycGDhwoDnDFqLKqbm5nG5xj1nW3fDoERQnp0otM3jwYCZPnkxaWhqKorBv3z7j32+xvn37miyzevVqateuzcmTJwkPDzeWjxs3jj59+hjn09PT6du3L02bNgUwGe0RQoiymDWZ+e233xg8eDAXL17E3d2dZs2asX37drp06QLAiy++SG5uLqNGjeLKlSu0bt2ar776Su4xI4SZ1apVi+joaNauXYuqqkRHR1OrVi2TNmfOnGHatGkcPHiQS5cuGW8NkJ6ebpLMtGzZ0mS5sWPH8swzz/DVV1/RuXNn+vbtS7Nmze78mxJCWC2zJjOrV6++ab2iKEyfPp3p06ffnYCEMBPF0ZGGR4+Ybd23Iy4ujtGjRwOwbNmyUvUxMTFotVpWrVqFv78/BoOB8PBwdDqdSTtnZ2eT+eHDh9OtWzc+//xzvvrqK+Lj45k/fz5jxoy5rTiFENWfxZ0zI8S/kaIolT7UY27du3c3JibdunUzqbt8+TKnTp1ixYoV3H///QDs3bu3wn1rtVpGjhzJyJEjmTx5MqtWrZJkRghRLklmhBC3xdbWllOnThlfl1SzZk28vLxYuXIlfn5+pKen89JLL1Wo33HjxtGjRw8aNGjAlStX+PrrrwkLC6vy+IUQ1YckM0KI21beLcZtbGzYuHEjY8eOJTw8nIYNG7JkyRI6dOhwyz71ej3PPvss58+fx83Nje7du7Nw4cIqjlwIUZ1Y3LOZqtqdejaTELdLnoEjLJF8LoWlqcz+2+amtUIIIYQQFk6SGSGEEEJYNUlmhBBCCGHVJJkRQgghhFWTZEYIIYQQVk2SGSGEEEJYNUlmhBBCCGHVJJkRQgghhFWTZEYIIYQQVk2SGSGE2aSlpaEoCsePHzd3KEIIKybJjBCiwmJjY3n44YerrD+tVsvFixcJDw8HIDk5GUVR+PPPP6tsHUKI6k8eNCmEMBtbW1t8fX2rvF9VVdHr9dSoIV9xQvwbyMiMEBZAVVUK8vVmmW73WbP16tVj0aJFJmWRkZFMnz7dOK8oCsuXL6dHjx44OjoSFBTEhx9+aKwveZgpLS2Njh07AlCzZk0URSE2Nta4febMmUNwcDCOjo5ERESwefNmYz/FIzpffvklLVu2RKPR8M0339zW+xJCWB/52SKEBSjUGVj53G6zrPupxe2x09jesf6nTZvGm2++yeLFi3n//fd5/PHHCQ8PJywszKSdVqvlo48+om/fvpw+fRo3NzccHR0BmDp1Klu2bGH58uWEhoayZ88ennjiCby9vWnfvr2xjxdffJF58+YRHByMh4fHHXtPQgjLIsmMEOKO6tevH8OHDwfg9ddfZ8eOHbz11lu8/fbbJu1sbW3x9PQEoHbt2sZkJCcnhwULFvD111/Tpk0bAIKDg9m7dy8rVqwwSWZee+01unTpchfelRDCkkgyI4QFqGFvw1OL29+64R1a951UnICUnK/M1UsnT54kLy+vVJKi0+lo3ry5SVnLli1vO04hhPWSZEYIC6Aoyh091HMn2NjYlDrfpqCgoELLKopS4fUYDAYAPv/8c+rUqWNSp9FoTOadnZ0r3K8QovqQZEYIcVu8vb25ePGicT47O5vU1NRS7Q4ePMiQIUNM5m8cUSlmb28PgF6vN5Y1btwYjUZDenq6ySElIYQoJsmMEOK2PPjgg6xZs4aYmBhq1qzJtGnTsLUtPbr04Ycf0rJlS9q1a0dSUhKHDh1i9erVZfYZGBiIoih89tln9OzZE0dHR1xdXZk4cSLPP/88BoOBdu3akZ2dzf79+3FxcWHo0KF3+q0KISycJDNCiAozGAzGe7dMnjyZs2fP0qtXL9zd3Xn99dfLHJmZMWMGGzduZNSoUfj6+pKUlETjxo3L7L9OnTrMmDGDl156iSeffJIhQ4awZs0aXn/9dWrXrk18fDxnz57Fw8ODFi1aMGXKlDv6foUQ1kFRb/cmE1YiOzsbd3d3srKycHNzM3c4QpCXl0dqaipBQUE4ODiYO5xK6d69OyEhISxdurRC7RVFYevWrVV612BxZ1jz51JUT5XZf8tN84QQt3TlyhU+//xzkpOT6dy5s7nDEUIIE3KYSQhxS3FxcRw+fJgJEybw0EMPmTscIYQwIcmMEOKWtm7delvLVfOj2EIICyGHmYQQQghh1SSZEUIIIYRVk2RGCCGEEFZNkhkhhBBCWDWzJjPx8fHce++9uLq6Urt2bR5++GFOnz5t0iY2NhZFUUym++67z0wRCyGEEMLSmDWZ2b17N88++ywHDx5kx44dFBYW0rVrV3Jyckzade/enYsXLxqnL774wkwRCyGEEMLSmPXS7O3bt5vMJyYmUrt2bY4cOcIDDzxgLNdoNPj6+t7t8IQQQghhBSzqnJmsrCwAPD09TcqTk5OpXbs2DRo0YMSIEWRmZpbbR35+PtnZ2SaTEKLqZGZm8vTTTxMQEGD8odGtWzcOHDhQoeXXrFmDh4fHnQ1SCPGvYjE3zVNVlfHjx9OuXTvCw8ON5T169KBfv34EBgaSmprKtGnTePDBBzly5AgajaZUP/Hx8cyYMeNuhi7Ev0rfvn0pKChg7dq1BAcH89tvv/Hf//6XP/74467HUlBQgJ2d3V1frxDCwqgWYtSoUWpgYKB67ty5m7a7cOGCamdnp3700Udl1ufl5alZWVnG6dy5cyqgZmVl3Ymwhai03Nxc9eTJk2pubq6xzGAwqLrcXLNMBoOhwrFfuXJFBdTk5ORy28yfP18NDw9XnZyc1Lp166rPPPOM+tdff6mqqqq7du1SAZPp1VdfVVVVVQF169atJn25u7uriYmJqqqqampqqgqomzZtUtu3b69qNBo1ISFBHTp0qPrQQw+pc+fOVX19fVVPT0911KhRqk6nq/D7EmV/LoUwp6ysrArvvy1iZGbMmDF8+umn7Nmzh7p16960rZ+fH4GBgfz8889l1ms0mjJHbISwZIX5+SwZ+qhZ1j127WbsKviUZBcXF1xcXPj444+57777yvxbs7GxYcmSJdSrV4/U1FRGjRrFiy++yNtvv01UVBSLFi3ilVdeMV656OLiUql4J02axPz580lMTESj0bB792527dqFn58fu3bt4pdffmHAgAFERkYyYsSISvUthLBOZj1nRlVVRo8ezZYtW/j6668JCgq65TKXL1/m3Llz+Pn53YUIhRAl1ahRgzVr1rB27Vo8PDxo27YtU6ZM4YcffjC2GTduHB07diQoKIgHH3yQ119/nQ8++AAAe3t73N3dURQFX19ffH19K53MjBs3jj59+hAUFIS/vz8ANWvWZOnSpTRq1IhevXoRHR3Nf//736p740IIi2bWkZlnn32W9evX88knn+Dq6kpGRgYA7u7uODo6cvXqVaZPn07fvn3x8/MjLS2NKVOmUKtWLR555BFzhi5Elaqh0TB27Wazrbsy+vbtS3R0NN988w0HDhxg+/btzJkzh3fffZfY2Fh27drFrFmzOHnyJNnZ2RQWFpKXl0dOTg7Ozs7/ON6WLVuWKmvSpAm2trbGeT8/P1JSUv7xuoQQ1sGsyczy5csB6NChg0l5YmIisbGx2NrakpKSwnvvvceff/6Jn58fHTt2ZNOmTbi6upohYiHuDEVRKnyoxxI4ODjQpUsXunTpwiuvvMLw4cN59dVX6dixIz179mTkyJG8/vrreHp6snfvXoYNG0ZBQcFN+1QUpdRTtstapqyE6MaTgBVFwWAw3MY7E0JYI7MmMzd+cd3I0dGRL7/88i5FI4S4XY0bN+bjjz/mu+++o7CwkPnz52NjU3QUu/gQUzF7e3v0en2pPry9vbl48aJx/ueff+batWt3NnAhRLVgEScACyGsw+XLl+nXrx9xcXE0a9YMV1dXvvvuO+bMmcNDDz1E/fr1KSws5K233iImJoZ9+/bxzjvvmPRRr149rl69yn//+18iIiJwcnLCycmJBx98kKVLl3LfffdhMBiYNGmSXHYthKgQi7ppnhDCsrm4uNC6dWsWLlzIAw88QHh4ONOmTWPEiBEsXbqUyMhIFixYwOzZswkPDycpKYn4+HiTPqKiohg5ciQDBgzA29ubOXPmADB//ny0Wi0PPPAAAwcOZOLEiTg5OZnjbQohrIyi3upYj5XLzs7G3d2drKws3NzczB2OEOTl5ZGamkpQUBAOVnSejKje5HMpLE1l9t8yMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQogK69ChA+PGjStV/vHHH6Moyt0PSAghkGRGCHGHFRQUmDsEIUQ1J8mMEBZAVVUMOr1Zpqp+1uz06dOJjIwkISGB4OBgNBoNqqqyfft22rVrh4eHB15eXvTq1YszZ84Yl0tLS0NRFLZs2ULHjh1xcnIiIiKCAwcOlOq7pEWLFlGvXj3jfHJyMq1atcLZ2RkPDw/atm3Lr7/+WqXvUQhhWWqYOwAhBKgFBi68st8s6/Z/LQrF3rZK+/zll1/44IMP+Oijj7C1Leo7JyeH8ePH07RpU3JycnjllVd45JFHOH78ODY2f/+uevnll5k3bx6hoaG8/PLLPP744/zyyy/UqHHrr6vCwkIefvhhRowYwYYNG9DpdBw6dEgOgQlRzUkyI4Socjqdjvfffx9vb29jWd++fU3arF69mtq1a3Py5EnCw8ON5RMnTiQ6OhqAGTNm0KRJE3755RcaNWp0y/VmZ2eTlZVFr169qF+/PgBhYWFV8ZaEEBZMkhkhLIBiZ4P/a1FmW3dVCwwMNElkAM6cOcO0adM4ePAgly5dwmAwAJCenm6SzDRr1sz42s/PD4DMzMwKJTOenp7ExsbSrVs3unTpQufOnenfv7+xHyFE9STnzAhhARRFwcbe1ixTZQ7BuLm5kZWVVar8zz//xM3NzTjv7Oxcqk1MTAyXL19m1apVfPvtt3z77bdA0ShOSXZ2dibbBTAmPjY2NqXO8bnxBOPExEQOHDhAVFQUmzZtokGDBhw8eLDC71EIYX0kmRFCVFijRo347rvvSpUfPnyYhg0blrvc5cuXOXXqFFOnTqVTp06EhYVx5cqVSq/f29ubjIwMk4Tm+PHjpdo1b96cyZMns3//fsLDw1m/fn2l1yWEsB6SzAghKmzUqFGcOXOGZ599lu+//56ffvqJZcuWsXr1al544YVyl6tZsyZeXl6sXLmSX375ha+//prx48dXev0dOnTg999/Z86cOZw5c4Zly5axbds2Y31qaiqTJ0/mwIED/Prrr3z11Vf89NNPct6MENWcJDNCiAqrV68e33zzDWfOnKFr167ce++9rFmzhjVr1tCvX79yl7OxsWHjxo0cOXKE8PBwnn/+eebOnVvp9YeFhfH222+zbNkyIiIiOHToEBMnTjTWOzk58b///Y++ffvSoEEDnnrqKUaPHs3TTz99W+9XCGEdFLWqbzJhYbKzs3F3dycrK8vkmL4Q5pKXl0dqaipBQUE4ODiYOxwhAPlcCstTmf23jMwIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqZk1m4uPjuffee3F1daV27do8/PDDnD592qSNqqpMnz4df39/HB0d6dChAydOnDBTxEIIIYSwNGZNZnbv3s2zzz7LwYMH2bFjB4WFhXTt2pWcnBxjmzlz5rBgwQKWLl3K4cOH8fX1pUuXLvz1119mjFwIIYQQlsKsycz27duJjY2lSZMmREREkJiYSHp6OkeOHAGKRmUWLVrEyy+/TJ8+fQgPD2ft2rVcu3ZNnoIrhBnExsaiKApvvvmmSfnHH3+MoihVui5FUfj444+rtE8hRPVU458sfPLkSdLT09HpdCblvXv3vq3+srKyAPD09ASKnoCbkZFB165djW00Gg3t27dn//79ZT48Lj8/n/z8fON8dnb2bcUihCibg4MDs2fP5umnn6ZmzZrmDkcIIW5vZObs2bNEREQQHh5OdHQ0Dz/8MA8//DCPPPIIjzzyyG0Foqoq48ePp127doSHhwOQkZEBgI+Pj0lbHx8fY92N4uPjcXd3N05arfa24hHiblJVFZ1OZ5apss+a7dy5M76+vsTHx5fbZv/+/TzwwAM4Ojqi1WoZO3asyeHjevXq8frrrzNw4EBcXFzw9/fnrbfeMqkHeOSRR1AUxTgfGxvLww8/bLKucePG0aFDB+N8hw4dGDt2LC+++CKenp74+voyffp0k2WysrJ46qmnqF27Nm5ubjz44IN8//33ldoOQgjLcVsjM8899xxBQUHs3LmT4OBgDh06xOXLl5kwYQLz5s27rUBGjx7NDz/8wN69e0vV3Th8rapquUPakydPZvz48cb57OxsSWiExSsoKGDWrFlmWfeUKVOwt7evcHtbW1tmzZrFwIEDGTt2LHXr1jWpT0lJoVu3brz++uusXr2a33//ndGjRzN69GgSExON7ebOncuUKVOYPn06X375Jc8//zyNGjWiS5cuHD58mNq1a5OYmEj37t2xtbWt1Htau3Yt48eP59tvv+XAgQPExsbStm1bunTpgqqqREdH4+npyRdffIG7uzsrVqygU6dO/PTTT8aRYSGE9bitkZkDBw7w2muv4e3tjY2NDTY2NrRr1474+HjGjh1b6f7GjBnDp59+yq5du0y+GH19fQFKjcJkZmaWGq0pptFocHNzM5mEEFXrkUceITIykldffbVU3dy5cxk4cCDjxo0jNDSUqKgolixZwnvvvUdeXp6xXdu2bXnppZdo0KABY8aM4dFHH2XhwoUAeHt7A+Dh4YGvr69xvqKaNWvGq6++SmhoKEOGDKFly5b897//BWDXrl2kpKTw4Ycf0rJlS0JDQ5k3bx4eHh5s3rz5djeJEMKMbmtkRq/X4+LiAkCtWrW4cOECDRs2JDAwsNSl1Tejqipjxoxh69atJCcnExQUZFIfFBSEr68vO3bsoHnz5gDodDp2797N7Nmzbyd0ISySnZ0dU6ZMMdu6b8fs2bN58MEHmTBhgkn5kSNH+OWXX0hKSjKWqaqKwWAgNTWVsLAwANq0aWOyXJs2bVi0aNFtxXKjZs2amcz7+fmRmZlpjO/q1at4eXmZtMnNzeXMmTNVsn4hxN11W8lMeHg4P/zwA8HBwbRu3Zo5c+Zgb2/PypUrCQ4OrnA/zz77LOvXr+eTTz7B1dXVOALj7u6Oo6MjiqIwbtw4Zs2aRWhoKKGhocyaNQsnJycGDhx4O6ELYZEURanUoR5L8MADD9CtWzemTJlCbGyssdxgMPD000+XOUobEBBw0z5vdUWUjY2N8Ryf4v8XX4BQ8tyfGjVqmMwrioLBYEBVVfR6PX5+fuzatQtuOF3Iw8MD1VC5c4jKVMnzkMzqeqhqoQHVoKL/S0dhfsWvTCuzZUULK7qaCrcrv6FJVfHMjc2VG17cUF+02O0ti3Lrz7e4fbeVzEydOtV4Mt/MmTPp1asX999/P15eXmzatKnC/SxfvhzA5OQ9gMTEROOX44svvkhubi6jRo3iypUrtG7dmq+++gpXV9fbCb3K5J35k7zTf/z9ZVjWd9eNX2hltCl18uWNbcr6UizVppy+bre/Eu3UCrSpaF+oZTS5ne12Q32pMEoWlFd34ypKvtFytu+t4lNvUmdcToUCRxV9RA0KLl3D1k5PKZXYD5ZqWu6yt+i0gus0XCvAkFeI7sJVAF4f/wr3do2ivl89AHT/d5XIRs1IOfoDAY6+pTv4XYcOHehV9ifvRTfs7/tF7d/1DQ3qhaA7X1RmZ2dHfuZV4zyAp4M7Kb9+b1J2/PBR7GrYUfB/RTGp+XoMOQXGeQBDbiEGu6KyZgFhZGRkoF7Kp5420DQ+HRRcuMq/UUGhDn1WPr9//AM1/rKiZMya3SwhUm4sLyow5kMVaFuUQJVRd7P2t+q7ZIclqp1b++H6gOn5c3fTbSUz3bp1M74ODg7m5MmT/PHHH9SsWbNSmWdFrqJQFIXp06eXuhrB3HTn/uLqnv8zdxjCChW6Kqh6Z9RCFRWDucOpHPX6dH30IrxhYx5/pD9vJ664Xq8ycdQ47u/dibFTnmfYwFicnJz4388/8d9vvmbR639fIHDgu2+Z9/YieneL5r97dvHRZx/z8ZoPjfWBdQP4em8ybVreh8benpoeNekQ9QAL3lnMus3raX1PKzZs2cSJ06eIbGJ6WOlmOt3fkftatKLf8IG8MXkGDeqHcvG3DLZ//SW9u/XinogWVbGlrI9yfbJRoHLnW/+tItl1RfKkf0suVe4PpzJ/5ZVbU5lV3CmGa4V3aU1l+0f3mSnp33YFgH1dV1weqHN9Tinrfzdk3bcamiyeV242W/YwalmZehl9K2UFVlbuWcl2pm+znPdZumGF2ikVGdItr+56LKa/NMpY/sb6G+pK/xuUs1wZ6y4VoqKQr9dxreASNWpqqKFx4KZu+tugnMrbHcmuwHKKYw2UfFtq+DgZy15/8w02f7YVADsfZ1r4tib5611MfWUqDz7aHVVVqR9cn/79+2Pn61y0kK3C+PHj+f5kCm8sehNXV1fmzZlH9GMPGWOZv2A+E16YSMKGtdSpU4fUM2eJfvwhpp6eypQ3XyUvL48nY59k8JDB/JjyI3Z+RX0r9rbYONsZ5wFsHGyxcaqBnX/RuX5f7NjOy1Nf5ulJo/n999/x9fXlgfsfoE54kLFNVW+723E3D0sY8vKocc0B3+fDcHC4xefSAtxyVLvMupuNvBYfvryNZf/RaPENL9Qbi0oUlLU+k/ZlrKecdatqZdvfvK2tmwZzUtTK3mTCymRnZ+Pu7k5WVpZc2SQsQl5eHqmpqQQFBVnFTuNOqFevHuPGjWPcuHHmDkVcJ59LYWkqs/+Wp2YLIYQQwqpJMiOEEEIIq1Zl58wIIURFpaWlmTsEIUQ1IiMzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqkswIIYQQwqpJMiOEqBLJyckoisKff/5p7lCEEP8ykswIISosNjYWRVFQFAU7OzuCg4OZOHEiOTk5t93fww8/XLVBCiH+deSmeUKISunevTuJiYkUFBTwzTffMHz4cHJychgwYIC5QxNC/EvJyIwQFkBVVfT6a2aZKvusWY1Gg6+vL1qtloEDBzJo0CA+/vjjUu2mT59OZGSkSdmiRYuoV6+esX7t2rV88sknxtGe5ORkAFJSUnjwwQdxdHTEy8uLp556iqtXr97GlhVC/BvIyIwQFsBgyCV5d1OzrLtD+xRsbZ1ue3lHR0cKCgoqvdzEiRM5deoU2dnZJCYmAuDp6cm1a9fo3r079913H4cPHyYzM5Phw4czevRo1qxZc9txCiGqL0lmhBC37dChQ6xfv55OnTpVelkXFxccHR3Jz8/H19fXWL527Vpyc3N57733cHZ2BmDp0qXExMQwe/ZsfHx8qix+IUT1IMmMEBbAxsaRDu1TzLbuyvjss89wcXGhsLCQgoICHnroId566y1OnjxZJfGcOnWKiIgIYyID0LZtWwwGA6dPn5ZkRghRiiQzQlgARVH+0aGeu6ljx44sX74cOzs7/P39sbOzAyiVzNjY2JQ6H6cih6NUVUVRlDLryisXQvy7yQnAQohKcXZ2JiQkhMDAQGMiUxZvb28yMjJMEprjx4+btLG3t0ev15uUNW7cmOPHj5tc7r1v3z5sbGxo0KBB1bwJIUS1IsmMEOKO6NChA7///jtz5szhzJkzLFu2jG3btpm0qVevHj/88AOnT5/m0qVLFBQUMGjQIBwcHBg6dCg//vgju3btYsyYMQwePFgOMQkhyiTJjBDijggLC+Ptt99m2bJlREREcOjQISZOnGjSZsSIETRs2JCWLVvi7e3Nvn37cHJy4ssvv+SPP/7g3nvv5dFHH6VTp04sXbrUTO9ECGHpFLWyN5mwMtnZ2bi7u5OVlYWbm5u5wxGCvLw8UlNTCQoKwsHBwdzhCAHI51JYnsrsv2VkRgghhBBWTZIZIYQQQlg1SWaEEEIIYdUkmRFCCCGEVZNkRgghhBBWTZIZIYQQQlg1SWaEEEIIYdUkmRFCCCGEVZNkRgghhBBWzazJzJ49e4iJicHf3x9FUfj4449N6mNjY1EUxWS67777zBOsEMJi1KtXj0WLFpk7jApJS0tDUZRSD9kUQlQdsyYzOTk5RERE3PSZK927d+fixYvG6YsvvriLEQohSir+gTFy5MhSdaNGjUJRFGJjY+94HIcPH+app56qsv6KE47iyd7enpCQEGbOnImlPfElPz+fMWPGUKtWLZydnenduzfnz583afPGG28QFRWFk5MTHh4e5glUiLuohjlX3qNHD3r06HHTNhqNBl9f37sUkRDiVrRaLRs3bmThwoU4OjoCRc/12bBhAwEBAXclBm9v7zvS786dO2nSpAn5+fns3buX4cOH4+fnx7Bhw+7I+m7HuHHj+M9//sPGjRvx8vJiwoQJ9OrViyNHjmBrawuATqejX79+tGnThtWrV5s5YiHuPIs/ZyY5OZnatWvToEEDRowYQWZm5k3b5+fnk52dbTIJYelUVSVHrzfLVNmRhxYtWhAQEMCWLVuMZVu2bEGr1dK8eXNj2fbt22nXrh0eHh54eXnRq1cvzpw5Y6zX6XSMHj0aPz8/HBwcqFevHvHx8cb66dOnExAQgEajwd/fn7FjxxrrSh5mevzxx3nsscdMYiwoKKBWrVokJiYat++cOXMIDg7G0dGRiIgINm/eXOq9eXl54evrS2BgIIMGDSIqKoqjR48a6w0GA6+99hp169ZFo9EQGRnJ9u3bTfo4dOgQzZs3x8HBgZYtW3Ls2DFjnaqqhISEMG/ePJNlfvzxR2xsbEy2T1mysrJYvXo18+fPp3PnzjRv3px169aRkpLCzp07je1mzJjB888/T9OmTW/anxDVhVlHZm6lR48e9OvXj8DAQFJTU5k2bRoPPvggR44cQaPRlLlMfHw8M2bMuMuRCvHPXDMYqL8nxSzrPvNAU5yv/6KvqCeffJLExEQGDRoEQEJCAnFxcSQnJxvb5OTkMH78eJo2bUpOTg6vvPIKjzzyCMePH8fGxoYlS5bw6aef8sEHHxAQEMC5c+c4d+4cAJs3b2bhwoVs3LiRJk2akJGRwffff19mLIMGDaJ///5cvXoVFxcXAL788ktycnLo27cvAFOnTmXLli0sX76c0NBQ9uzZwxNPPIG3tzft27cvs9/vvvuOo0ePMnToUGPZ4sWLmT9/PitWrKB58+YkJCTQu3dvTpw4QWhoKDk5OfTq1YsHH3yQdevWkZqaynPPPWdcXlEU4uLiSExMZOLEicbyhIQE7r//furXr3/T7X7kyBEKCgro2rWrsczf35/w8HD2799Pt27dbrq8ENWVRSczAwYMML4ODw+nZcuWBAYG8vnnn9OnT58yl5k8eTLjx483zmdnZ6PVau94rEL8mwwePJjJkycbzzXZt28fGzduNElmihOJYqtXr6Z27dqcPHmS8PBw0tPTCQ0NpV27diiKQmBgoLFteno6vr6+dO7cGTs7OwICAmjVqlWZsXTr1g1nZ2e2bt3K4MGDAVi/fj0xMTG4ubmRk5PDggUL+Prrr2nTpg0AwcHB7N27lxUrVpgkM1FRUdjY2KDT6SgoKOCpp55iyJAhxvp58+YxadIk40jQ7Nmz2bVrF4sWLWLZsmUkJSWh1+tJSEjAycmJJk2acP78eZ555hljH08++SSvvPIKhw4dolWrVhQUFLBu3Trmzp17y+2ekZGBvb09NWvWNCn38fEhIyPjlssLUV1ZdDJzIz8/PwIDA/n555/LbaPRaModtRHCUjnZ2HDmAfMcEnCyqfzR5lq1ahEdHc3atWtRVZXo6Ghq1apl0ubMmTNMmzaNgwcPcunSJQwGA1CUqISHhxMbG0uXLl1o2LAh3bt3p1evXsYRh379+rFo0SKCg4Pp3r07PXv2JCYmhho1Sn9l2dnZ0a9fP5KSkhg8eDA5OTl88sknrF+/HoCTJ0+Sl5dHly5dTJbT6XQmh8UANm3aRFhYGAUFBaSkpDB27Fhq1qzJm2++SXZ2NhcuXKBt27Ymy7Rt29Y4anTq1CkiIiJwcnIy1hcnUMX8/PyIjo4mISGBVq1a8dlnn5GXl0e/fv0qvP1vpKoqiqLc9vJCWDurSmYuX77MuXPn8PPzM3coQlQpRVEqfajH3OLi4hg9ejQAy5YtK1UfExODVqtl1apV+Pv7YzAYCA8PR6fTAUXn3qSmprJt2zZ27txJ//796dy5M5s3b0ar1XL69Gl27NjBzp07GTVqFHPnzmX37t3Y2dmVWtegQYNo3749mZmZ7NixAwcHB+PFBcVJ1Oeff06dOnVMlrvxh49WqyUkJASAsLAwzp49y7Rp05g+fbqxzY1JQ8lEoqLnHw0fPpzBgwezcOFCEhMTGTBggEkCVB5fX190Oh1XrlwxGZ3JzMwkKiqqQusWojoy6wnAV69e5fjx48b7L6SmpnL8+HHS09O5evUqEydO5MCBA6SlpZGcnExMTAy1atXikUceMWfYQgiKbpug0+nQ6XSlztW4fPkyp06dYurUqXTq1ImwsDCuXLlSqg83NzcGDBjAqlWr2LRpEx999BF//PEHAI6OjvTu3ZslS5aQnJzMgQMHSEkp+7yiqKgotFotmzZtIikpiX79+mFvbw9A48aN0Wg0pKenExISYjLd6hC0ra0thYWF6HQ63Nzc8Pf3Z+/evSZt9u/fT1hYmHFd33//Pbm5ucb6gwcPluq3Z8+eODs7s3z5crZt20ZcXNxN4yh2zz33YGdnx44dO4xlFy9e5Mcff5RkRvyrmXVk5rvvvqNjx47G+eJzXYYOHcry5ctJSUnhvffe488//8TPz4+OHTuyadMmXF1dzRWyEOI6W1tbTp06ZXxdUs2aNfHy8mLlypX4+fmRnp7OSy+9ZNJm4cKF+Pn5ERkZiY2NDR9++CG+vr54eHiwZs0a9Ho9rVu3xsnJiffffx9HR0eT82pKUhSFgQMH8s477/DTTz+xa9cuY52rqysTJ07k+eefx2Aw0K5dO7Kzs9m/fz8uLi4mJ/hevnyZjIwMCgsLSUlJYfHixXTs2BE3NzcAXnjhBV599VXq169PZGQkiYmJHD9+nKSkJAAGDhzIyy+/zLBhw5g6dSppaWmlrlwq3l6xsbFMnjyZkJCQUoeiyuPu7s6wYcOYMGECXl5eeHp6MnHiRJo2bUrnzp2N7dLT0/njjz9IT09Hr9cbfzCGhIQYT5IWolpRq7msrCwVULOysswdihCqqqpqbm6uevLkSTU3N9fcoVTa0KFD1Yceeqjc+oceekgdOnSoqqqqumPHDjUsLEzVaDRqs2bN1OTkZBVQt27dqqqqqq5cuVKNjIxUnZ2dVTc3N7VTp07q0aNHVVVV1a1bt6qtW7dW3dzcVGdnZ/W+++5Td+7caVxPYGCgunDhQpN1nzhxQgXUwMBA1WAwmNQZDAZ18eLFasOGDVU7OzvV29tb7datm7p7925VVVU1NTVVBYyTra2tWrduXXXEiBFqZmamsR+9Xq/OmDFDrVOnjmpnZ6dGRESo27ZtM1nXgQMH1IiICNXe3l6NjIxUP/roIxVQjx07ZtLuzJkzKqDOmTPnVpvdRG5urjp69GjV09NTdXR0VHv16qWmp6ebtBk6dKjJ+ymedu3addN+rfVzKaqnyuy/FVW1sNtbVrHs7Gzc3d3Jysoy/roSwpzy8vJITU0lKCgIBwcHc4cjzGTfvn106NCB8+fP4+PjY+5w5HMpLE5l9t9WdQKwEEJYu/z8fM6dO8e0adPo37+/RSQyQlg7i78DsBBCVCcbNmygYcOGZGVlMWfOHJO6pKQkXFxcypyaNGlipoiFsHwyMiOEEHdRbGxsuQ/j7N27N61bty6zrqxL0oUQRSSZEUIIC+Hq6ipXawpxG+QwkxBCCCGsmiQzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgirU69ePRYtWmTuMCokLS0NRVGMz0cSQlQ9SWaEEBUWGxuLoiiMHDmyVN2oUaNQFKXce6hUpcOHD/PUU09VWX/FCUfxZG9vT0hICDNnzsTSnviSn5/PmDFjqFWrFs7OzvTu3Zvz588b69PS0hg2bBhBQUE4OjpSv359Xn31VXQ6nRmjFuLOkmRGCFEpWq2WjRs3kpubayzLy8tjw4YNBAQE3JUYvL29cXJyqvJ+d+7cycWLF/n555+ZMWMGb7zxBgkJCVW+nn9i3LhxbN26lY0bN7J3716uXr1Kr1690Ov1APzvf//DYDCwYsUKTpw4wcKFC3nnnXeYMmWKmSMX4s6RZEYIC6CqKtd0hWaZKjvy0KJFCwICAtiyZYuxbMuWLWi1Wpo3b24s2759O+3atcPDwwMvLy969erFmTNnjPU6nY7Ro0fj5+eHg4MD9erVIz4+3lg/ffp0AgIC0Gg0+Pv7M3bsWGNdycNMjz/+OI899phJjAUFBdSqVYvExETj9p0zZw7BwcE4OjoSERHB5s2bS703Ly8vfH19CQwMZNCgQURFRXH06FFjvcFg4LXXXqNu3bpoNBoiIyPZvn27SR+HDh2iefPmODg40LJlS44dO2asU1WVkJAQ5s2bZ7LMjz/+iI2Njcn2KUtWVharV69m/vz5dO7cmebNm7Nu3TpSUlLYuXMnAN27dycxMZGuXbsSHBxM7969mThxosm/lxDVjdwBWAgLkFugp/ErX5pl3Sdf64aTfeW+Cp588kkSExMZNGgQAAkJCcTFxZGcnGxsk5OTw/jx42natCk5OTm88sorPPLIIxw/fhwbGxuWLFnCp59+ygcffEBAQADnzp3j3LlzAGzevJmFCxeyceNGmjRpQkZGBt9//32ZsQwaNIj+/ftz9epVXFxcAPjyyy/Jycmhb9++AEydOpUtW7awfPlyQkND2bNnD0888QTe3t60b9++zH6/++47jh49ytChQ41lixcvZv78+axYsYLmzZuTkJBA7969OXHiBKGhoeTk5NCrVy8efPBB1q1bR2pqKs8995xxeUVRiIuLIzExkYkTJxrLExISuP/++6lfv/5Nt/uRI0coKCiga9euxjJ/f3/Cw8PZv38/3bp1K3O5rKwsPD09b9q3ENZMkhkhRKUNHjyYyZMnG8812bdvHxs3bjRJZooTiWKrV6+mdu3anDx5kvDwcNLT0wkNDaVdu3YoikJgYKCxbXp6Or6+vnTu3Bk7OzsCAgJo1apVmbF069YNZ2dntm7dyuDBgwFYv349MTExuLm5kZOTw4IFC/j6669p06YNAMHBwezdu5cVK1aYJDNRUVHY2Nig0+koKCjgqaeeYsiQIcb6efPmMWnSJONI0OzZs9m1axeLFi1i2bJlJCUlodfrSUhIwMnJiSZNmnD+/HmeeeYZYx9PPvkkr7zyCocOHaJVq1YUFBSwbt065s6de8vtnpGRgb29PTVr1jQp9/HxISMjo8xlzpw5w1tvvcX8+fNv2b8Q1kqSGSEsgKOdLSdfK/tX9d1Yd2XVqlWL6Oho1q5di6qqREdHU6tWLZM2Z86cYdq0aRw8eJBLly5hMBiAokQlPDyc2NhYunTpQsOGDenevTu9evUyjjj069ePRYsWERwcTPfu3enZsycxMTHUqFH6K8vOzo5+/fqRlJTE4MGDycnJ4ZNPPmH9+vUAnDx5kry8PLp06WKynE6nMzksBrBp0ybCwsIoKCggJSWFsWPHUrNmTd58802ys7O5cOECbdu2NVmmbdu2xlGjU6dOERERYXI+T3ECVczPz4/o6GgSEhJo1aoVn332GXl5efTr16/C2/9GqqqiKEqp8gsXLtC9e3f69evH8OHDb7t/ISydJDNCWABFUSp9qMfc4uLiGD16NADLli0rVR8TE4NWq2XVqlX4+/tjMBgIDw83XlXTokULUlNT2bZtGzt37qR///507tyZzZs3o9VqOX36NDt27GDnzp2MGjWKuXPnsnv37jKfHj1o0CDat29PZmYmO3bswMHBgR49egAYk6jPP/+cOnXqmCyn0WhM5rVaLSEhIQCEhYVx9uxZpk2bxvTp041tbkwaSiYSFT3/aPjw4QwePJiFCxeSmJjIgAEDKnRCs6+vLzqdjitXrpiMzmRmZhIVFWXS9sKFC3Ts2JE2bdqwcuXKCsUlhLWSE4CFELele/fu6HQ6dDpdqXM1Ll++zKlTp5g6dSqdOnUiLCyMK1eulOrDzc2NAQMGsGrVKjZt2sRHH33EH3/8AYCjoyO9e/dmyZIlJCcnc+DAAVJSUsqMJSoqCq1Wy6ZNm0hKSqJfv37Y29sD0LhxYzQaDenp6YSEhJhMWq32pu/R1taWwsJCdDodbm5u+Pv7s3fvXpM2+/fvJywszLiu77//3uRKr4MHD5bqt2fPnjg7O7N8+XK2bdtGXFzcTeMods8992BnZ8eOHTuMZRcvXuTHH380SWb+7//+jw4dOtCiRQsSExOxsZGvelG9WddPQSGExbC1teXUqVPG1yXVrFkTLy8vVq5ciZ+fH+np6bz00ksmbRYuXIifnx+RkZHY2Njw4Ycf4uvri4eHB2vWrEGv19O6dWucnJx4//33cXR0NDmvpiRFURg4cCDvvPMOP/30E7t27TLWubq6MnHiRJ5//nkMBgPt2rUjOzub/fv34+LiYnKC7+XLl8nIyKCwsJCUlBQWL15Mx44dcXNzA+CFF17g1VdfpX79+kRGRpKYmMjx48dJSkoCYODAgbz88ssMGzaMqVOnkpaWVurKpeLtFRsby+TJkwkJCSl1KKo87u7uDBs2jAkTJuDl5YWnpycTJ06kadOmdO7cGSgakenQoQMBAQHMmzeP33//3bi8r69vhdYjhNVRq7msrCwVULOysswdihCqqqpqbm6uevLkSTU3N9fcoVTa0KFD1Yceeqjc+oceekgdOnSoqqqqumPHDjUsLEzVaDRqs2bN1OTkZBVQt27dqqqqqq5cuVKNjIxUnZ2dVTc3N7VTp07q0aNHVVVV1a1bt6qtW7dW3dzcVGdnZ/W+++5Td+7caVxPYGCgunDhQpN1nzhxQgXUwMBA1WAwmNQZDAZ18eLFasOGDVU7OzvV29tb7datm7p7925VVVU1NTVVBYyTra2tWrduXXXEiBFqZmamsR+9Xq/OmDFDrVOnjmpnZ6dGRESo27ZtM1nXgQMH1IiICNXe3l6NjIxUP/roIxVQjx07ZtLuzJkzKqDOmTPnVpvdRG5urjp69GjV09NTdXR0VHv16qWmp6cb6xMTE03eS8npVv1a6+dSVE+V2X8rqmpht7esYtnZ2bi7u5OVlWX8dSWEOeXl5ZGamkpQUBAODg7mDkeYyb59++jQoQPnz5/Hx8fH3OHI51JYnMrsv+UwkxBC3EX5+fmcO3eOadOm0b9/f4tIZISwdnJWmBBC3EUbNmygYcOGZGVlMWfOHJO6pKQkXFxcypyaNGlipoiFsHwyMiOEEHdRbGxsuQ/j7N27N61bty6zrqxL0oUQRSSZEUIIC+Hq6oqrq6u5wxDC6shhJiGEEEJYNUlmhBBCCGHVJJkRQgghhFWTZEYIIYQQVk2SGSGEEEJYNUlmhBBWp169eixatMjcYVRIWloaiqJw/Phxc4ciRLUlyYwQosJiY2NRFIWRI0eWqhs1ahSKopR7D5WqdPjwYZ566qkq66844Sie7O3tCQkJYebMmVjaE1/y8/MZM2YMtWrVwtnZmd69e3P+/HmTNr179yYgIAAHBwf8/PwYPHgwFy5cMFPEQtx5Zk1m9uzZQ0xMDP7+/iiKwscff2xSr6oq06dPx9/fH0dHRzp06MCJEyfME6wQAgCtVsvGjRvJzc01luXl5bFhwwYCAgLuSgze3t44OTlVeb87d+7k4sWL/Pzzz8yYMYM33niDhISEKl/PPzFu3Di2bt3Kxo0b2bt3L1evXqVXr17o9Xpjm44dO/LBBx9w+vRpPvroI86cOcOjjz5qxqiFuLPMmszk5OQQERHB0qVLy6yfM2cOCxYsYOnSpRw+fBhfX1+6dOnCX3/9dZcjFeIOU1XQ5ZhnquTIQ4sWLQgICGDLli3Gsi1btqDVamnevLmxbPv27bRr1w4PDw+8vLzo1asXZ86cMdbrdDpGjx6Nn58fDg4O1KtXj/j4eGP99OnTCQgIQKPR4O/vz9ixY411JQ8zPf744zz22GMmMRYUFFCrVi0SExOvb16VOXPmEBwcjKOjIxEREWzevLnUe/Py8sLX15fAwEAGDRpEVFQUR48eNdYbDAZee+016tati0ajITIyku3bt5v0cejQIZo3b46DgwMtW7bk2LFjxjpVVQkJCWHevHkmy/z444/Y2NiYbJ+yZGVlsXr1aubPn0/nzp1p3rw569atIyUlhZ07dxrbPf/889x3330EBgYSFRXFSy+9xMGDBykoKLhp/0JYK7PeAbhHjx706NGjzDpVVVm0aBEvv/wyffr0AWDt2rX4+Piwfv16nn766TKXy8/PJz8/3zifnZ1d9YELUdUKrsEsf/Ose8oFsHeu1CJPPvkkiYmJDBo0CICEhATi4uJITk42tsnJyWH8+PE0bdqUnJwcXnnlFR555BGOHz+OjY0NS5Ys4dNPP+WDDz4gICCAc+fOce7cOQA2b97MwoUL2bhxI02aNCEjI4Pvv/++zFgGDRpE//79uXr1Ki4uLgB8+eWX5OTk0LdvXwCmTp3Kli1bWL58OaGhoezZs4cnnngCb29v2rdvX2a/3333HUePHmXo0KHGssWLFzN//nxWrFhB8+bNSUhIoHfv3pw4cYLQ0FBycnLo1asXDz74IOvWrSM1NZXnnnvOuLyiKMTFxZGYmMjEiRON5QkJCdx///3Ur1//ptv9yJEjFBQU0LVrV2OZv78/4eHh7N+/n27dupVa5o8//iApKYmoqCh5JIKotiz2nJnU1FQyMjJM/mg1Gg3t27dn//795S4XHx+Pu7u7cdJqtXcjXCH+VQYPHszevXtJS0vj119/Zd++fTzxxBMmbfr27UufPn0IDQ0lMjKS1atXk5KSwsmTJwFIT08nNDSUdu3aERgYSLt27Xj88ceNdb6+vnTu3JmAgABatWrFiBEjyoylW7duODs7s3XrVmPZ+vXriYmJwc3NjZycHBYsWEBCQgLdunUjODiY2NhYnnjiCVasWGHSV1RUFC4uLtjb23PvvffSv39/hgwZYqyfN28ekyZN4rHHHqNhw4bMnj2byMhI4yhRUlISer2ehIQEmjRpQq9evXjhhRdM1vHkk09y+vRpDh06BBSNIq1bt464uLhbbveMjAzs7e2pWbOmSbmPjw8ZGRkmZZMmTcLZ2RkvLy/S09P55JNPbtm/ENbKYp/NVPyH6ePjY1Lu4+PDr7/+Wu5ykydPZvz48cb57OxsSWiE5bNzKhohMde6K6lWrVpER0ezdu1aVFUlOjqaWrVqmbQ5c+YM06ZN4+DBg1y6dAmDwQAUJSrh4eHExsbSpUsXGjZsSPfu3enVq5fxx0u/fv1YtGgRwcHBdO/enZ49exITE0ONGqW/suzs7OjXrx9JSUkMHjyYnJwcPvnkE9avXw/AyZMnycvLo0uXLibL6XQ6k8NiAJs2bSIsLIyCggJSUlIYO3YsNWvW5M033yQ7O5sLFy7Qtm1bk2Xatm1rHDU6deoUERERJufztGnTxqS9n58f0dHRJCQk0KpVKz777DPy8vLo169fhbf/jVRVRVEUk7IXXniBYcOG8euvvzJjxgyGDBnCZ599VqqdENWBxSYzxW78wyvrj7YkjUaDRqO502EJUbUUpdKHeswtLi6O0aNHA7Bs2bJS9TExMWi1WlatWoW/vz8Gg4Hw8HB0Oh1QdO5Namoq27ZtY+fOnfTv35/OnTuzefNmtFotp0+fZseOHezcuZNRo0Yxd+5cdu/eXeahkkGDBtG+fXsyMzPZsWMHDg4OxkPYxUnU559/Tp06dUyWu/G7QqvVEhISAkBYWBhnz55l2rRpTJ8+3djmZt9JFb3yafjw4QwePJiFCxeSmJjIgAEDKnRCs6+vLzqdjitXrpiMzmRmZhIVFWXStlatWtSqVYsGDRoQFhaGVqvl4MGDpZIrIaoDi01mfH19gaIRGj8/P2N5ZmZmqdEaczj25Wcc+uT6CYQlv8BKvDb5WrteXt6XnVqJPkq3KVleRn8l6m/5XXuTBiq3WPhm1bdY8S37/gcU/sEv0X/yK7aMRRUUHGt60uzRJ7hkZ4udrW3VdF2hMP/5L/K8q3+Rfy2HzLSztGjUgLzrVzQ1bxhKZtpZ8q/lkFfDlv8dO8KpU6d4c/qrNK0fBMC3330HwJ+ZGWT+etbYZ8f7WtHxvlZ0vr8djw2N5fT3x/D08ACgTURT2kQ05fGHHyKqU2f2fLWdZuHhGAoLufrHZX7/NRWA0Dp+1PHzY/Xyt/lvcjK9uncj6+L/AVDbxQmNvT0/Hj1C43o3XG1lKOT39DQu/1/Rpc1/XLzA754exurc7CwKCwu5cOZnXF1d8fXx4cvP/kNYiX72JO+ieUQEv6enofWtzXtrj5H+0/9wdHAAYMcXn5fq+97wxjg6ODBv1hts27aNTz/YxKX0NMr9N7peXM/HGzs7Oz7asJ6HY3oBkPFbJj/++CMvvzCBS+fKHrW+fP2y7Mxzv3KpbtnnZhUUFnL1ymW2fvAeur/KOdewnL+Hm36yyvsbKrevm/RWblc3i6B0XZnNyygsM5Z/tGwFY7lh2VLvr9TsrddVej03dnKLdZahZJuGbe4nvGOXm7S+syw2mQkKCsLX15cdO3YYh4J1Oh27d+9m9uzZZo4OdLm5XL18ydxhCCtkp9OBqqIaDKi3mSyZ684nqqqiqioGvR4F2LN9G1D0tWjQ6431bi4u1KzpwXvr1+Pt5cX/XbzAG3OLruBRDQYMhXpWJCTiU9ubJmFh2NjY8Olnn1Hb2xtXJyfWb/oAvV5Pi8gIHB0c2PTRR0X3TPHxQV9YiErRiIu+sNAY28MxvVizLomzaWlsXveesc7RwYGRw4cx7fXX0RcW0Oqelvx19SrfHT2Ks7MT/fv0wVBYdJXP5Uu/c/HCBQoL9Zz66TSrEhJpe999ODk4oC8o4Jnhw5i3eAkBdevQJCyMTZs/4seTp1g6fx76ggIe7tmTWXPmMW7iCzw3ahTn/+88b69YCYChsAB9iauJ+vd5hJmz51AvMIDmzZpSWIErjZwcHHi836O8+vpM3F1d8HD34LU33ySsYQPatmpFoU7Hse+/59j3P9Cq5T24u7uTnn6OOYsXUy8ggMjwphReHxm7UaFej6FQz58ZF7n2h3y3icqpHXTzk9fvNLMmM1evXuWXX34xzqempnL8+HE8PT0JCAhg3LhxzJo1i9DQUEJDQ5k1axZOTk4MHDjQjFEXCe/QmXrNShxvL7FTKi+jNZaXbGvaoJw+lLJfltPGWGyyHpMFy4yvYtW3Wrbyv6oq0eA2/IPd/j9ZtMyRqKIyXUEBmVeyqOlXB4frv94rtdI7lsmU7vjGEo2zM/YFhXjVLRqZ8Lqh3t7RCY2TM94B9Vi/LonnJ0ygY3QvGjRowML58+nctSuuXt541Q3Ap05d3lm5kl9++QVbW1ta3nMP//nPf/AOqEedoGDmzpvLjDdno9frCQ9vwsdbthDatBkAtra2OHt44FX37/Ph4p56iiXL3yEwIIAevR82+SzOmTefwOD6LFu5kokvT8PDw4PmkZFMevFFPOtoyS4sOhTVf0issX8/X196REfz2vTpeHp7A/DilJfR29jw+uy5ZP7+O2Fhjfho84e0bNsOAE/g461bGT12DF0fepiwsEa8+WY8Ax4fiHttHzzr1DVu2GdGj2HJ8neIixuGp3/dCv8rvbV0GS9NmcLIcePJzc2lY4cOrFn7Ht516wIqPn9cYceSt5i/dBk5OTn4+vrStUsXJk+ahF+dv0dlbvy3zc/P58/cfLo9M44aNqX/Fss9hHaTz2O5I67l9nWzkeHyKio5mlxW0U3+Xm+5qjIKy15vBWO5sfCG5Up3U1act+jjVrGV+Z5uvoyXNrDUMneToprx9pbJycl07NixVPnQoUNZs2YNqqoyY8YMVqxYwZUrV2jdujXLli0jPDy8wuvIzs7G3d2drKws3NzcqjJ8IW5LXl4eqampBAUFlZHMiH+Lffv20aFDB86fP28Rh87lcyksTWX232ZNZu4GSWaEpZGdxr9bfn4+586d46mnnsLPz4+kpCRzhwTI51JYnsrsvy32PjNCCFEdbdiwgYYNG5KVlcWcOXNM6pKSknBxcSlzatKkiZkiFsLyWewJwEIIUR3FxsaW+zDO3r1707p16zLr5O69QpRPkhkhhLAQrq6uuLq6mjsMIayOHGYSQgghhFWTZEYIIYQQVk2SGSGEEEJYNUlmhBBCCGHVJJkRQgghhFWTZEYIYXXq1avHokWLzB1GhaSlpaEoCsePHzd3KEJUW5LMCCEqLDY2FkVRGDlyZKm6UaNGoShKufdQqUqHDx/mqaeeqrL+ihOO4sne3p6QkBBmzpxZ/jOJzCQ/P58xY8ZQq1YtnJ2d6d27N+fPny+3bWRkpCRTotqTZEYIUSlarZaNGzeSm5trLMvLy2PDhg0EBATclRi8vb1xcnKq8n537tzJxYsX+fnnn5kxYwZvvPEGCQkJVb6ef2LcuHFs3bqVjRs3snfvXq5evUqvXr3Q6/Wl2r744ov4+/uX0YsQ1YskM0JYAFVVuVZwzSxTZUceWrRoQUBAAFu2bDGWbdmyBa1WS/Pmfz9Jfvv27bRr1w4PDw+8vLzo1asXZ86cMdbrdDpGjx6Nn58fDg4O1KtXj/j4eGP99OnTCQgIQKPR4O/vz9ixY411JQ8zPf744zz22GMmMRYUFFCrVi0SExON23fOnDkEBwfj6OhIREQEmzdvLvXevLy88PX1JTAwkEGDBhEVFcXRo0eN9QaDgddee426deui0WiIjIxk+/btJn0cOnSI5s2b4+DgQMuWLTl27JixTlVVQkJCmDdvnskyP/74IzY2NibbpyxZWVmsXr2a+fPn07lzZ5o3b866detISUlh586dJm23bdvGV199VWpdQlRHcgdgISxAbmEurdeXfRv7O+3bgd/iZFe5UY4nn3ySxMREBg0aBEBCQgJxcXEkJycb2+Tk5DB+/HiaNm1KTk4Or7zyCo888gjHjx/HxsaGJUuW8Omnn/LBBx8QEBDAuXPnOHfuHACbN29m4cKFbNy4kSZNmpCRkcH3339fZiyDBg2if//+XL16FRcXFwC+/PJLcnJy6Nu3LwBTp05ly5YtLF++nNDQUPbs2cMTTzyBt7c37du3L7Pf7777jqNHjzJ06FBj2eLFi5k/fz4rVqygefPmJCQk0Lt3b06cOEFoaCg5OTn06tWLBx98kHXr1pGamspzzz1nXF5RFOLi4khMTGTixInG8oSEBO6//37q169/0+1+5MgRCgoK6Nq1q7HM39+f8PBw9u/fT7du3QD47bffGDFiBB9//PEdGcESwtJIMiOEqLTBgwczefJk47km+/btY+PGjSbJTHEiUWz16tXUrl2bkydPEh4eTnp6OqGhobRr1w5FUQgMDDS2TU9Px9fXl86dO2NnZ0dAQACtWrUqM5Zu3brh7OzM1q1bGTx4MADr168nJiYGNzc3cnJyWLBgAV9//TVt2rQBIDg4mL1797JixQqTZCYqKgobGxt0Oh0FBQU89dRTDBkyxFg/b948Jk2aZBwJmj17Nrt27WLRokUsW7aMpKQk9Ho9CQkJODk50aRJE86fP88zzzxj7OPJJ5/klVde4dChQ7Rq1YqCggLWrVvH3Llzb7ndMzIysLe3p2bNmiblPj4+ZGRkAEWjP7GxsYwcOZKWLVuSlpZ2y36FsHaSzAhhARxrOPLtwG/Ntu7KqlWrFtHR0axduxZVVYmOjqZWrVombc6cOcO0adM4ePAgly5dwmAwAEWJSnh4OLGxsXTp0oWGDRvSvXt3evXqZRxx6NevH4sWLSI4OJju3bvTs2dPYmJiqFGj9FeWnZ0d/fr1IykpicGDB5OTk8Mnn3zC+vXrATh58iR5eXl06dLFZDmdTmdyWAxg06ZNhIWFUVBQQEpKCmPHjqVmzZq8+eabZGdnc+HCBdq2bWuyTNu2bY2jRqdOnSIiIsJkNKQ4gSrm5+dHdHQ0CQkJtGrVis8++4y8vDz69etX4e1/I1VVURQFgLfeeovs7GwmT5582/0JYW0kmRHCAiiKUulDPeYWFxfH6NGjAVi2bFmp+piYGLRaLatWrcLf3x+DwUB4eDg6nQ4oOvcmNTWVbdu2sXPnTvr370/nzp3ZvHkzWq2W06dPs2PHDnbu3MmoUaOYO3cuu3fvLvPp0YMGDaJ9+/ZkZmayY8cOHBwc6NGjB4Axifr888+pU6eOyXIajcZkXqvVEhISAkBYWBhnz55l2rRpTJ8+3dimOGkoVjKRqOj5R8OHD2fw4MEsXLiQxMREBgwYUKHDQb6+vuh0Oq5cuWIyOpOZmUlUVBQAX3/9NQcPHiz13lq2bMmgQYNYu3ZthWIUwprICcBCiNvSvXt3dDodOp3OeK5GscuXL3Pq1CmmTp1Kp06dCAsL48qVK6X6cHNzY8CAAaxatYpNmzbx0Ucf8ccffwDg6OhI7969WbJkCcnJyRw4cICUlJQyY4mKikKr1bJp0yaSkpLo168f9vb2ADRu3BiNRkN6ejohISEmk1arvel7tLW1pbCwEJ1Oh5ubG/7+/uzdu9ekzf79+wkLCzOu6/vvvze50uvgwYOl+u3ZsyfOzs4sX76cbdu2ERcXd9M4it1zzz3Y2dmxY8cOY9nFixf58ccfjcnMkiVL+P777zl+/DjHjx/niy++AIpGnd54440KrUcIayMjM0KI22Jra8upU6eMr0uqWbMmXl5erFy5Ej8/P9LT03nppZdM2ixcuBA/Pz8iIyOxsbHhww8/xNfXFw8PD9asWYNer6d169Y4OTnx/vvv4+joaHJeTUmKojBw4EDeeecdfvrpJ3bt2mWsc3V1ZeLEiTz//PMYDAbatWtHdnY2+/fvx8XFxeQE38uXL5ORkUFhYSEpKSksXryYjh074ubmBsALL7zAq6++Sv369YmMjCQxMZHjx4+TlJQEwMCBA3n55ZcZNmwYU6dOJS0trcyriWxtbYmNjWXy5MmEhISUOhRVHnd3d4YNG8aECRPw8vLC09OTiRMn0rRpUzp37gxQ6vL44pOi69evT926dSu0HiGsjSQzQojbVryTv5GNjQ0bN25k7NixhIeH07BhQ5YsWUKHDh2MbVxcXJg9ezY///wztra23HvvvXzxxRfY2Njg4eHBm2++yfjx49Hr9TRt2pT//Oc/eHl5lRvLoEGDmDVrFoGBgaXOa3n99depXbs28fHxnD17Fg8PD1q0aMGUKVNM2hUnBLa2tvj5+dGzZ0+T0YyxY8eSnZ3NhAkTyMzMpHHjxnz66aeEhoYa39N//vMfRo4cSfPmzWncuDGzZ88udTI0wLBhw5g1a1aFR2WKLVy4kBo1atC/f39yc3Pp1KkTa9asKZVQCvFvoqiWdnvLKpadnY27uztZWVnlfvEKcTfl5eWRmppKUFAQDg4O5g5HmMm+ffvo0KED58+fx8fHx9zhyOdSWJzK7L9lZEYIIe6i/Px8zp07x7Rp0+jfv79FJDJCWDs5AVgIIe6iDRs20LBhQ7KyspgzZ45JXVJSEi4uLmVOTZo0MVPEQlg+GZkRQoi7KDY2ttyHcfbu3ZvWrcu+E3RZl6QLIYpIMiOEEBbC1dUVV1dXc4chhNWRw0xCCCGEsGqSzAghhBDCqkkyI4QQQgirJsmMEEIIIayaJDNCCCGEsGqSzAghrE69evVYtGiRucOokLS0NBRF4fjx4+YORYhqS5IZIUSFxcbGoigKI0eOLFU3atQoFEUp9x4qVenw4cM89dRTVdZfccJRPNnb2xMSEsLMmTOxtCe+5OfnM2bMGGrVqoWzszO9e/fm/PnzJm3q1atn8n4URSn1oE8hqhNJZoQQlaLVatm4cSO5ubnGsry8PDZs2FDqic13ire3N05OTlXe786dO7l48SI///wzM2bM4I033iAhIaHK1/NPjBs3jq1bt7Jx40b27t3L1atX6dWrF3q93qTda6+9xsWLF43T1KlTzRSxEHeeRd80b/r06cyYMcOkzMfHh4yMDDNF9LdfT1zm7PHfzR2GsEI29gbcggq5eiUPnb1l/eq/lYJ8Pc3CI0j7NY2ktRvp3+8xAD7cvIk6/nWpF1iPgnw92Zdz2fnfr5i7YDanTp3ExtaWVi1b8easeQQHBQOg0+mYMm0Sn/7nY/7M+hOf2j7EDh3GhHEvABA/eybr1r9H5u+ZeNb05KHejzAnfj4ATZs34pmnRzNq5GjiRgxFVVUS333v7zgLCmjQJJjXp7/BEwOHoKoqi99aSMLad/nttwxC6ofywoSXeLj3IwD8dSUPAI2tM0527jjZQUz3PrS6910O7j9Ev4cHAmAwGJg7/03WvJfApcuXaNigIdOnvU7nTl2N6z5y9DDPTRjDTz+dJqxRYyaOfxGAq3/mkXXpGs1bNSVu6HCeGz3OuMzJUydo80Arjh1KMW6fsmRlZ7F69WpWvr2a1s2Lngz+zpJVhEU04NMtn9P5wS4AqAYVOxsHnO3cjcuq+fBXfm6Z/Rb9e+STl1PAof+cxaCT37nmYl3fCH8LbOxFcHNvs63fopMZgCZNmrBz507jvKU85v7Sub84+c0Fc4chrJCDuw1N/dzJv1Zo3Gmoqgr5eeYJSOOAoigVaqovMKAvVOnfZyDvrVtL7x59AVj7/hoG9B3E/oN70RcYyLtawJ9/ZPPUk6MIa9iEa9dymL1wFgMHD+DrL/ZiY2PD2yvf4ottn7Ny6Rrq+NflwsX/4/8u/B95Vwv4zxcfs+ydpaxYkkDDBo3I/D2TE6dSyLtaAIBqKEqs8q4W8HCvR3nq2Vgu/3YFZ2cXAL7673au5eTQrWM0eVcLmDX3Nb7Y/h9mvzafoKD6HPx2P0+NjMPNyYOo+9qRn1PUr+5aoXEdx384yvffH6Pfw48Zy955dxlvvb2EuW8spGmTCDZ88D6PPdGPPV99S3BQfXKu5dDv8b60a/MAS+evJP3cr7w8bZKx7/ycQh57dBDrkt7jqdhnjds1cc0a7rs3Cj9vLbnX11WWbw8epqCggKhWDxjbebh606hBY/bt3U/bVh2M22fh4vnMmRePv19dYqIf4tmnnsPe3r7cvgsKCynI0/Pzd3+Ql2Wo0OdBiGIOznaSzNxMjRo18PX1rXD7/Px88vPzjfPZ2dl3Iiz8QzxoFRN0R/oW1UO5+YGtnhpOuTi62aOx1wBgyM3lfPcOdy22kuomH8DG0aFCbWvY21DDzoahTw5h1twZXMq6iKIoHD7yLe+/l8ShI/upYW+Ds4eGAY/3N1l2RcOVBIVoSc84S5PGTfjt0kVCQkJ4sEsHFEWhUXiose3vf2Tg6+NDj17dsLOzoyEh3N8hyliv2IDGsQbOHhp69e6J04vO/Peb7Tz+2CAAPt32ET16ROOr9SYnJ4cVq5fx+afbad3qPgCaNGvE0R8OsX7zWrp074RTVtFOvtejXbGxsUGn01FQUMCTscN4cliscb3vrF7K+HETeGJI0Xoi7pnNgcP7SEhawcJ5i9n08fuoBgOrVr6Lk5MT97SK5HJWJuPGj8HR1R5nDw3DhsUxZ8EsTv3yAy3vuZeCggI++uQD3nhtFs7umptu/+yrf2Bvb0+dANPvRF8/H/7IumRc/tlRo4mIiKSme02+O3qYV197hYu/nWfZknfK7TtfB/ZZNWjWsS7oLeNHo7AefvXdb93oDrL4ZObnn3/G398fjUZD69atmTVrFsHB5Q/DxsfHlzo0dSf4hXjgF+Jxx9cjqp+8vDxSU1NxcrXHweF6MmOnv8VSd46zmz02TjffiRarYW+LrZ0NgcF1iI6O5sOtG1FVlejoaAKD62BrZ0MNe1uc3TWcOXOGadOmcfDgQS5duoTBUPRr/9KVDJzdWzDi6WF06dKFFq2a0b17d3r16kXXrkWHawYNeZy3VyylafMwunfvTs+ePYmJiaFGjaKvLMVGwd6xxvWdt4b+/fuxeesHDH86jpycHD7/4jPWr1+Ps7uGkz/9QF5eHr0fiTZ5LzqdjubNm+PsrsHJrej9b9q0ibCwMAoKCkhJSWHs2LHU9qnFm2++SXZ2NhcvXqBjp/YmScf9D7Tj+++/x9ldw9m0X4iIjMDbr6axvsOD9wMUJTPuGuq71yM6OpoNH66jfad2bN36Bfn5eTwxdCBOt/h30DgXPWzS2cO0nY2tgr1DDWP5pCkvGOtat2uJb53aPProo8xfOA8vL68y+7bNU9E41qBRh7o4OFQsuRXCUlh0MtO6dWvee+89GjRowG+//cbMmTOJiorixIkT5f5BTp48mfHjxxvns7Oz0Wq1dytkIW6L4uhIw6NHzLbu2xEXF8fo0aMBWLZsWan6mJgYtFotq1atwt/fH4PBQHh4ODqdDoAWLVqQmprKtm3b2LlzJ/3796dz585s3rwZrVbL6dOn2bFjBzt37mTUqFHMnTuX3bt3l/n06EGDBtG+fXsyMzPZsWMHDg4O9OjRA8CYRH3++efUqVPHZDmNxjQp0Gq1hISEABAWFsbZs2eZNm0a06dP/3t73TDkpqqqsayiVz4NHz6cwYMHs3DhQhITExkwYECFTmj29fVFp9Nx5coVatb8O2HKzMwkKiqq3OXuu69oROqXX34p97tTCGtm0clM8ZcRQNOmTWnTpg3169dn7dq1JglLSRqNptQXlBCWTlEUlDtwdc6d1L17d2Ni0q1bN5O6y5cvc+rUKVasWMH99xeNTOzdu7dUH25ubgwYMIABAwbw6KOP0r17d/744w88PT1xdHSkd+/e9O7dm2effZZGjRqRkpJCixYtSvUTFRWFVqtl06ZNbNu2jX79+hnPD2ncuDEajYb09HTat29fqfdoa2tLYWEhOp0ONzc3/P392bt3Lw888ICxzf79+2nVqpVxXe+//z65ubk4Xk8SDx48WKrfnj174uzszPLly9m2bRt79uypUDz33HMPdnZ27Nixg/79iw7jXbx4kR9//JE5c+aUu9yxY8cA8PPzq9gbF8LKWHQycyNnZ2eaNm3Kzz//bO5QhPjXs7W15dSpU8bXJdWsWRMvLy9WrlyJn58f6enppe5zsnDhQvz8/IiMjMTGxoYPP/wQX19fPDw8WLNmDXq9ntatW+Pk5MT777+Po6MjgYGBZcaiKAoDBw7knXfe4aeffmLXrl3GOldXVyZOnMjzzz+PwWCgXbt2ZGdns3//flxcXBg6dKix7eXLl8nIyKCwsJCUlBQWL15Mx44dcXNzA+CFF17g1VdfpX79+kRGRpKYmMjx48dJSkoCYODAgbz88ssMGzaMqVOnkpaWxrx588rcdrGxsUyePJmQkBDatGlToW3u7u7OsGHDmDBhAl5eXnh6ejJx4kSaNm1K586dAThw4AAHDx6kY8eOuLu7c/jwYZ5//nl69+591y6dF+KuU61IXl6eWqdOHXXGjBkVXiYrK0sF1KysrDsYmRAVl5ubq548eVLNzc01dyiVNnToUPWhhx4qt/6hhx5Shw4dqqqqqu7YsUMNCwtTNRqN2qxZMzU5OVkF1K1bt6qqqqorV65UIyMjVWdnZ9XNzU3t1KmTevToUVVVVXXr1q1q69atVTc3N9XZ2Vm977771J07dxrXExgYqC5cuNBk3SdOnFABNTAwUDUYDCZ1BoNBXbx4sdqwYUPVzs5O9fb2Vrt166bu3r1bVVVVTU1NVSm6KlYFVFtbW7Vu3brqiBEj1MzMTGM/er1enTFjhlqnTh3Vzs5OjYiIULdt22ayrgMHDqgRERGqvb29GhkZqX700UcqoB47dsyk3ZkzZ1RAnTNnzq02u4nc3Fx19OjRqqenp+ro6Kj26tVLTU9PN9YfOXJEbd26teru7q46ODioDRs2VF999VU1Jyfnlv1a6+dSVE+V2X8rqmpht7csYeLEicTExBAQEEBmZiYzZ85k9+7dpKSklPsL7UbZ2dm4u7uTlZVl/HUlhDkVnwAcFBQkJ1r+i+3bt48OHTpw/vx5fHx8zB2OfC6FxanM/tuiDzOdP3+exx9/nEuXLuHt7c19993HwYMHK5zICCGEpcnPz+fcuXNMmzaN/v37W0QiI4S1s+hkZuPGjeYOQQghqtSGDRsYNmwYkZGRvP/++yZ1SUlJPP3002UuFxgYyIkTJ+5GiEJYHYtOZoQQorqJjY0t92GcvXv3pnXr1mXWlXVJuhCiiCQzQghhIVxdXXF1dTV3GEJYHXmamBBmYsHn3ot/Ifk8CmsmyYwQd1nx4YJr166ZORIh/lb8eZTDWcIayWEmIe4yW1tbPDw8yMzMBMDJyanCT60Woqqpqsq1a9fIzMzEw8Oj1A0QhbAGkswIYQbFT4IvTmiEMDcPDw/j51IIayPJjBBmoCgKfn5+1K5dm4KCAnOHI/7l7OzsZERGWDVJZoQwI1tbW9mJCCHEPyQnAAshhBDCqkkyI4QQQgirJsmMEEIIIayaJDNCCCGEsGqSzAghhBDCqkkyI4QQQgirJsmMEEIIIayaJDNCCCGEsGqSzAghhBDCqkkyI4QQQgirJsmMEEIIIayaJDNCCCGEsGqSzAghhBDCqslTs2/TH++9z6UVK8wdhhBCCGF2noMHU2vk02ZbvyQzt8mQl4f+8mVzhyGEEEKYneHaNbOuX5KZ22Tr6YlT69bmDkMIAaCYO4B/AdXcAQhLVqN2bfOu36xrt2L6P/7g2rffmjsMIYQQwuwcIyLMun5JZm7T783q8n+jeps7DCGEEMLsCiPqYc6xGUlmbtMau7/YVM/b3GEIIYQQZtfLIYelZly/JDO3KdfWnzzXBuYOQwghhDC7nBrZZl2/JDO3yeWCG057fjV3GEIIIYTZubWuC/eYb/2SzNymuvY1MMg9B4UQQgjq2Js3nZBk5jbVr3WCIfftMHcYQgghhNmF1b4faGq29Usyc5v2ff8hWzll7jCEEEIIs8u7/H/E3DvYbOu3imTm7bffZu7cuVy8eJEmTZqwaNEi7r//frPGVMcugOjfJZkRQggh6ngEmnX9iqqqFn1fx02bNjF48GDefvtt2rZty4oVK3j33Xc5efIkAQEBt1w+Ozsbd3d3srKycHNzuwsRCyGEEOKfqsz+2+KTmdatW9OiRQuWL19uLAsLC+Phhx8mPj6+VPv8/Hzy8/ON89nZ2Wi1WklmhBBCCCtSmWTGoi/H0el0HDlyhK5du5qUd+3alf3795e5THx8PO7u7sZJq9XejVCFEEIIYSYWncxcunQJvV6Pj4+PSbmPjw8ZGRllLjN58mSysrKM07lz5+5GqEIIIYQwE6s4AVhRTB+Jq6pqqbJiGo0GjUZzN8ISQgghhAWw6JGZWrVqYWtrW2oUJjMzs9RojRBCCCH+nSw6mbG3t+eee+5hxw7Tm9Pt2LGDqKgoM0UlhBBCCEti8YeZxo8fz+DBg2nZsiVt2rRh5cqVpKenM3LkSHOHJoQQQggLYPHJzIABA7h8+TKvvfYaFy9eJDw8nC+++ILAQPPeoEcIIYQQlsHi7zPzT8lN84QQQgjrU23uMyOEEEIIcSuSzAghhBDCqkkyI4QQQgirJsmMEEIIIayaJDNCCCGEsGoWf2n2P1V8sVZ2draZIxFCCCFERRXvtyty0XW1T2b++usvAHl6thBCCGGF/vrrL9zd3W/aptrfZ8ZgMHDhwgVcXV3LfTjl7crOzkar1XLu3Dm5h80dJNv57pDtfHfIdr47ZDvfHXdyO6uqyl9//YW/vz82Njc/K6baj8zY2NhQt27dO7oONzc3+WO5C2Q73x2yne8O2c53h2znu+NObedbjcgUkxOAhRBCCGHVJJkRQgghhFWTZOYf0Gg0vPrqq2g0GnOHUq3Jdr47ZDvfHbKd7w7ZzneHpWznan8CsBBCCCGqNxmZEUIIIYRVk2RGCCGEEFZNkhkhhBBCWDVJZoQQQghh1SSZuU1vv/02QUFBODg4cM899/DNN9+YO6RqJT4+nnvvvRdXV1dq167Nww8/zOnTp80dVrUXHx+PoiiMGzfO3KFUS//3f//HE088gZeXF05OTkRGRnLkyBFzh1WtFBYWMnXqVIKCgnB0dCQ4OJjXXnsNg8Fg7tCs2p49e4iJicHf3x9FUfj4449N6lVVZfr06fj7++Po6EiHDh04ceLEXYtPkpnbsGnTJsaNG8fLL7/MsWPHuP/+++nRowfp6enmDq3a2L17N88++ywHDx5kx44dFBYW0rVrV3JycswdWrV1+PBhVq5cSbNmzcwdSrV05coV2rZti52dHdu2bePkyZPMnz8fDw8Pc4dWrcyePZt33nmHpUuXcurUKebMmcPcuXN56623zB2aVcvJySEiIoKlS5eWWT9nzhwWLFjA0qVLOXz4ML6+vnTp0sX4fMQ7ThWV1qpVK3XkyJEmZY0aNVJfeuklM0VU/WVmZqqAunv3bnOHUi399ddfamhoqLpjxw61ffv26nPPPWfukKqdSZMmqe3atTN3GNVedHS0GhcXZ1LWp08f9YknnjBTRNUPoG7dutU4bzAYVF9fX/XNN980luXl5anu7u7qO++8c1dikpGZStLpdBw5coSuXbualHft2pX9+/ebKarqLysrCwBPT08zR1I9Pfvss0RHR9O5c2dzh1Jtffrpp7Rs2ZJ+/fpRu3ZtmjdvzqpVq8wdVrXTrl07/vvf//LTTz8B8P3337N371569uxp5siqr9TUVDIyMkz2ixqNhvbt29+1/WK1f9BkVbt06RJ6vR4fHx+Tch8fHzIyMswUVfWmqirjx4+nXbt2hIeHmzucamfjxo0cPXqUw4cPmzuUau3s2bMsX76c8ePHM2XKFA4dOsTYsWPRaDQMGTLE3OFVG5MmTSIrK4tGjRpha2uLXq/njTfe4PHHHzd3aNVW8b6vrP3ir7/+eldikGTmNimKYjKvqmqpMlE1Ro8ezQ8//MDevXvNHUq1c+7cOZ577jm++uorHBwczB1OtWYwGGjZsiWzZs0CoHnz5pw4cYLly5dLMlOFNm3axLp161i/fj1NmjTh+PHjjBs3Dn9/f4YOHWru8Ko1c+4XJZmppFq1amFra1tqFCYzM7NUVir+uTFjxvDpp5+yZ88e6tata+5wqp0jR46QmZnJPffcYyzT6/Xs2bOHpUuXkp+fj62trRkjrD78/Pxo3LixSVlYWBgfffSRmSKqnl544QVeeuklHnvsMQCaNm3Kr7/+Snx8vCQzd4ivry9QNELj5+dnLL+b+0U5Z6aS7O3tueeee9ixY4dJ+Y4dO4iKijJTVNWPqqqMHj2aLVu28PXXXxMUFGTukKqlTp06kZKSwvHjx41Ty5YtGTRoEMePH5dEpgq1bdu21O0FfvrpJwIDA80UUfV07do1bGxMd222trZyafYdFBQUhK+vr8l+UafTsXv37ru2X5SRmdswfvx4Bg8eTMuWLWnTpg0rV64kPT2dkSNHmju0auPZZ59l/fr1fPLJJ7i6uhpHwtzd3XF0dDRzdNWHq6trqfOQnJ2d8fLykvOTqtjzzz9PVFQUs2bNon///hw6dIiVK1eycuVKc4dWrcTExPDGG28QEBBAkyZNOHbsGAsWLCAuLs7coVm1q1ev8ssvvxjnU1NTOX78OJ6engQEBDBu3DhmzZpFaGgooaGhzJo1CycnJwYOHHh3Arwr10xVQ8uWLVMDAwNVe3t7tUWLFnLJcBUDypwSExPNHVq1J5dm3zn/+c9/1PDwcFWj0aiNGjVSV65cae6Qqp3s7Gz1ueeeUwMCAlQHBwc1ODhYffnll9X8/Hxzh2bVdu3aVeZ38tChQ1VVLbo8+9VXX1V9fX1VjUajPvDAA2pKSspdi09RVVW9O2mTEEIIIUTVk3NmhBBCCGHVJJkRQgghhFWTZEYIIYQQVk2SGSGEEEJYNUlmhBBCCGHVJJkRQgghhFWTZEYIIYQQVk2SGSGEEEJYNUlmhBAWKTk5GUVR+PPPP80dihDCwskdgIUQFqFDhw5ERkayaNEioOhBdX/88Qc+Pj4oimLe4IQQFk0eNCmEsEj29vb4+vqaOwwhhBWQw0xCCLOLjY1l9+7dLF68GEVRUBSFNWvWmBxmWrNmDR4eHnz22Wc0bNgQJycnHn30UXJycli7di316tWjZs2ajBkzBr1eb+xbp9Px4osvUqdOHZydnWndujXJycnmeaNCiDtCRmaEEGa3ePFifvrpJ8LDw3nttdcAOHHiRKl2165dY8mSJWzcuJG//vqLPn360KdPHzw8PPjiiy84e/Ysffv2pV27dgwYMACAJ598krS0NDZu3Ii/vz9bt26le/fupKSkEBoaelffpxDizpBkRghhdu7u7tjb2+Pk5GQ8tPS///2vVLuCggKWL19O/fr1AXj00Ud5//33+e2333BxcaFx48Z07NiRXbt2MWDAAM6cOcOGDRs4f/48/v7+AEycOJHt27eTmJjIrFmz7t6bFELcMZLMCCGshpOTkzGRAfDx8aFevXq4uLiYlGVmZgJw9OhRVFWlQYMGJv3k5+fj5eV1d4IWQtxxkswIIayGnZ2dybyiKGWWGQwGAAwGA7a2thw5cgRbW1uTdiUTICGEdZNkRghhEezt7U1O3K0KzZs3R6/Xk5mZyf3331+lfQshLIdczSSEsAj16tXj22+/JS0tjUuXLhlHV/6JBg0aMGjQIIYMGcKWLVtITU3l8OHDzJ49my+++KIKohZCWAJJZoQQFmHixInY2trSuHFjvL29SU9Pr5J+ExMTGTJkCBMmTKBhw4b07t2bb7/9Fq1WWyX9CyHMT+4ALIQQQgirJiMzQgghhLBqkswIIYQQwqpJMiOEEEIIqybJjBBCCCGsmiQzQgghhLBqkswIIcT/t1sHJAAAAACC/r9uR6ArBNZkBgBYkxkAYE1mAIA1mQEA1mQGAFgLiLFkSit7X7AAAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "sim.data.where(sim.data['particle_type'] == 'Massive Body',drop=True)['a'].plot(hue=\"name\")" - ] - }, { "cell_type": "code", "execution_count": null, @@ -137,9 +45,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python (My debug_env Kernel)", + "display_name": "swiftest", "language": "python", - "name": "debug_env" + "name": "swiftest" }, "language_info": { "codemirror_mode": { From 168d2de5a206d4ac7576f57efa8074070d1a2606 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 19:47:51 -0500 Subject: [PATCH 07/13] Fixed some inconsistencies with the new API and the id vs name dimensioning --- examples/helio_gr_test/init_cond.py | 48 ++++ examples/whm_gr_test/init_cond.py | 226 ++++++++++++++++++ .../whm_gr_test/swiftest_relativity.ipynb | 193 +++++++++++++++ python/swiftest/swiftest/io.py | 24 +- python/swiftest/swiftest/simulation_class.py | 77 ++++-- 5 files changed, 539 insertions(+), 29 deletions(-) create mode 100755 examples/helio_gr_test/init_cond.py create mode 100644 examples/whm_gr_test/init_cond.py create mode 100644 examples/whm_gr_test/swiftest_relativity.ipynb diff --git a/examples/helio_gr_test/init_cond.py b/examples/helio_gr_test/init_cond.py new file mode 100755 index 000000000..9e89dc358 --- /dev/null +++ b/examples/helio_gr_test/init_cond.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python3 +import swiftest + +sim = swiftest.Simulation() +sim.param['PL_IN'] = "pl.swiftest.in" +sim.param['TP_IN'] = "tp.swiftest.in" +sim.param['CB_IN'] = "cb.swiftest.in" +sim.param['BIN_OUT'] = "bin.swiftest.nc" + +sim.param['MU2KG'] = swiftest.MSun +sim.param['TU2S'] = swiftest.YR2S +sim.param['DU2M'] = swiftest.AU2M +sim.param['T0'] = 0.0 +sim.param['DT'] = 0.25 * swiftest.JD2S / swiftest.YR2S +sim.param['TSTOP'] = 1000.0 +sim.param['ISTEP_OUT'] = 1461 +sim.param['ISTEP_DUMP'] = 1461 +sim.param['CHK_QMIN_COORD'] = "HELIO" +sim.param['CHK_QMIN'] = swiftest.RSun / swiftest.AU2M +sim.param['CHK_QMIN_RANGE'] = f"{swiftest.RSun / swiftest.AU2M} 1000.0" +sim.param['CHK_RMIN'] = swiftest.RSun / swiftest.AU2M +sim.param['CHK_RMAX'] = 1000.0 +sim.param['CHK_EJECT'] = 1000.0 +sim.param['OUT_STAT'] = "UNKNOWN" +sim.param['IN_FORM'] = "EL" +sim.param['OUT_FORM'] = "XVEL" +sim.param['OUT_TYPE'] = "NETCDF_DOUBLE" +sim.param['RHILL_PRESENT'] = "YES" +sim.param['GR'] = 'YES' + +bodyid = { + "Sun": 0, + "Mercury": 1, + "Venus": 2, + "Earth": 3, + "Mars": 4, + "Jupiter": 5, + "Saturn": 6, + "Uranus": 7, + "Neptune": 8, +} + +for name, id in bodyid.items(): + sim.add(name, idval=id, date="2027-04-30") + +sim.save("param.swiftest.in") + + diff --git a/examples/whm_gr_test/init_cond.py b/examples/whm_gr_test/init_cond.py new file mode 100644 index 000000000..7904eb100 --- /dev/null +++ b/examples/whm_gr_test/init_cond.py @@ -0,0 +1,226 @@ +import numpy as np +import sys +from astroquery.jplhorizons import Horizons +import astropy.constants as const + +#Values from JPL Horizons +AU2M = const.au.value +GMSunSI = const.GM_sun.value +Rsun = const.R_sun.value +GC = const.G.value +JD = 86400 +year = 365.25 * JD +c = 299792458.0 +MSun_over_Mpl = [6023600.0, + 408523.71, + 328900.56, + 3098708., + 1047.3486, + 3497.898, + 22902.98, + 19412.24, + 1.35e8] + +MU2KG = GMSunSI / GC #Conversion from mass unit to kg +DU2M = AU2M #Conversion from radius unit to centimeters +TU2S = year #Conversion from time unit to seconds +GU = GC / (DU2M**3 / (MU2KG * TU2S**2)) + +GMSun = GMSunSI / (DU2M**3 / TU2S**2) + +t_print = 10.e0 * year / TU2S #output interval to print results +deltaT = 0.25 * JD / TU2S #timestep simulation +end_sim = 1.0e3 * year / TU2S + t_print #end time + +# Solar oblatenes values: From Mecheri et al. (2004), using Corbard (b) 2002 values (Table II) +J2 = 2.198e-7 * (Rsun / DU2M)**2 +J4 = -4.805e-9 * (Rsun / DU2M)**4 + +tstart = '2021-01-28' +tend = '2021-01-29' +tstep = '1d' +planetid = { + 'mercury' : '1', + 'venus' : '2', + 'earthmoon' : '3', + 'mars' : '4', + 'jupiter' : '5', + 'saturn' : '6', + 'uranus' : '7', + 'neptune' : '8', + 'plutocharon' : '9' +} +npl = 9 + +#Planet Msun/M ratio +MSun_over_Mpl = { + 'mercury' : 6023600.0, + 'venus' : 408523.71, + 'earthmoon' : 328900.56, + 'mars' : 3098708., + 'jupiter' : 1047.3486, + 'saturn' : 3497.898, + 'uranus' : 22902.98, + 'neptune' : 19412.24, + 'plutocharon' : 1.35e8 +} + +#Planet radii in meters +Rpl = { + 'mercury' : 2439.4e3, + 'venus' : 6051.8e3, + 'earthmoon' : 6371.0084e3, # Earth only for radius + 'mars' : 3389.50e3, + 'jupiter' : 69911e3, + 'saturn' : 58232.0e3, + 'uranus' : 25362.e3, + 'neptune' : 24622.e3, + 'plutocharon' : 1188.3e3 +} + +pdata = {} +plvec = {} +Rhill = {} + +for key,val in planetid.items(): + pdata[key] = Horizons(id=val, id_type='majorbody',location='@sun', + epochs={'start': tstart, 'stop': tend, + 'step': tstep}) + plvec[key] = np.array([pdata[key].vectors()['x'][0], + pdata[key].vectors()['y'][0], + pdata[key].vectors()['z'][0], + pdata[key].vectors()['vx'][0], + pdata[key].vectors()['vy'][0], + pdata[key].vectors()['vz'][0] + ]) + Rhill[key] = pdata[key].elements()['a'][0] * (3 * MSun_over_Mpl[key])**(-1.0 / 3.0) + + +if __name__ == '__main__': + # Convert from AU-day to AU-year just because I find it easier to keep track of the sim progress + for plid in plvec: + plvec[plid][3:] *= year / JD + + # Names of all output files + swifter_input = "param.swifter.in" + swifter_pl = "pl.swifter.in" + swifter_tp = "tp.swifter.in" + swifter_bin = "bin.swifter.dat" + swifter_enc = "enc.swifter.dat" + + swiftest_input = "param.swiftest.in" + swiftest_pl = "pl.swiftest.in" + swiftest_tp = "tp.swiftest.in" + swiftest_cb = "cb.swiftest.in" + swiftest_bin = "bin.swiftest.dat" + swiftest_enc = "enc.swiftest.dat" + + # Simulation start, stop, and output cadence times + t_0 = 0 # simulation start time + end_sim = 1000.0e0 * year / TU2S # simulation end time + deltaT = 0.25 * JD / TU2S # simulation step size + t_print = 1.0 * year / TU2S #output interval to print results + + iout = int(np.ceil(t_print / deltaT)) + rmin = Rsun / DU2M + rmax = 1000.0 + + #Make Swifter files + plfile = open(swifter_pl, 'w') + print(f'{npl+1} ! Planet input file generated using init_cond.py using JPL Horizons data for the major planets (and Pluto) for epoch {tstart}' ,file=plfile) + print(f'1 {GMSun}',file=plfile) + print(f'0.0 0.0 0.0',file=plfile) + print(f'0.0 0.0 0.0',file=plfile) + for i, plid in enumerate(plvec): + print(f'{i + 2} {GMSun / MSun_over_Mpl[plid]} {Rhill[plid]}', file=plfile) + print(f'{Rpl[plid] / DU2M}', file=plfile) + print(f'{plvec[plid][0]} {plvec[plid][1]} {plvec[plid][2]}', file=plfile) + print(f'{plvec[plid][3]} {plvec[plid][4]} {plvec[plid][5]}', file=plfile) + plfile.close() + + tpfile = open(swifter_tp, 'w') + print(0,file=tpfile) + tpfile.close() + + sys.stdout = open(swifter_input, "w") + print(f'! Swifter input file generated using init_cond.py') + print(f'T0 {t_0} ') + print(f'TSTOP {end_sim}') + print(f'DT {deltaT}') + print(f'PL_IN {swifter_pl}') + print(f'TP_IN {swifter_tp}') + print(f'IN_TYPE ASCII') + print(f'ISTEP_OUT {iout:d}') + print(f'ISTEP_DUMP {iout:d}') + print(f'BIN_OUT {swifter_bin}') + print(f'OUT_TYPE REAL8') + print(f'OUT_FORM EL') + print(f'OUT_STAT NEW') + print(f'J2 {J2}') + print(f'J4 {J4}') + print(f'CHK_CLOSE yes') + print(f'CHK_RMIN {rmin}') + print(f'CHK_RMAX {rmax}') + print(f'CHK_EJECT {rmax}') + print(f'CHK_QMIN {rmin}') + print(f'CHK_QMIN_COORD HELIO') + print(f'CHK_QMIN_RANGE {rmin} {rmax}') + print(f'ENC_OUT {swifter_enc}') + print(f'EXTRA_FORCE no') + print(f'BIG_DISCARD no') + print(f'RHILL_PRESENT yes') + print(f'C {c / (DU2M / TU2S)}') + + #Now make Swiftest files + cbfile = open(swiftest_cb, 'w') + print(f'{1.0}',file=cbfile) + print(f'{rmin}',file=cbfile) + print(f'{J2}',file=cbfile) + print(f'{J4}',file=cbfile) + + plfile = open(swiftest_pl, 'w') + print(npl,file=plfile) + + for i, plid in enumerate(plvec): + print(f'{i + 2} {1.0 / MSun_over_Mpl[plid]}', file=plfile) + print(f'{Rpl[plid] / DU2M}', file=plfile) + print(f'{plvec[plid][0]} {plvec[plid][1]} {plvec[plid][2]}', file=plfile) + print(f'{plvec[plid][3]} {plvec[plid][4]} {plvec[plid][5]}', file=plfile) + plfile.close() + tpfile = open(swiftest_tp, 'w') + print(0,file=tpfile) + tpfile.close() + + sys.stdout = open(swiftest_input, "w") + print(f'! Swiftest input file generated using init_cond.py') + print(f'T0 {t_0} ') + print(f'TSTOP {end_sim}') + print(f'DT {deltaT}') + print(f'CB_IN {swiftest_cb}') + print(f'PL_IN {swiftest_pl}') + print(f'TP_IN {swiftest_tp}') + print(f'IN_TYPE ASCII') + print(f'ISTEP_OUT {iout:d}') + print(f'ISTEP_DUMP {iout:d}') + print(f'BIN_OUT {swiftest_bin}') + print(f'OUT_TYPE REAL8') + print(f'OUT_FORM EL') + print(f'OUT_STAT REPLACE') + print(f'CHK_CLOSE yes') + print(f'CHK_RMIN {rmin}') + print(f'CHK_RMAX {rmax}') + print(f'CHK_EJECT {rmax}') + print(f'CHK_QMIN {rmin}') + print(f'CHK_QMIN_COORD HELIO') + print(f'CHK_QMIN_RANGE {rmin} {rmax}') + print(f'ENC_OUT {swiftest_enc}') + print(f'EXTRA_FORCE no') + print(f'BIG_DISCARD no') + print(f'ROTATION no') + print(f'GR yes') + print(f'MU2KG {MU2KG}') + print(f'DU2M {DU2M}') + print(f'TU2S {TU2S}') + + + sys.stdout = sys.__stdout__ diff --git a/examples/whm_gr_test/swiftest_relativity.ipynb b/examples/whm_gr_test/swiftest_relativity.ipynb new file mode 100644 index 000000000..ae586907c --- /dev/null +++ b/examples/whm_gr_test/swiftest_relativity.ipynb @@ -0,0 +1,193 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "import swiftest\n", + "from astroquery.jplhorizons import Horizons" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swifter file param.swifter.in\n", + "Reading in time 1.000e+03\n", + "Creating Dataset\n", + "Successfully converted 1001 output frames.\n", + "Swifter simulation data stored as xarray DataSet .ds\n" + ] + } + ], + "source": [ + "swiftersim = swiftest.Simulation(param_file=\"param.swifter.in\", codename=\"Swifter\")\n", + "swiftersim.bin2xr()\n", + "swifterdat = swiftersim.ds" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swiftest file param.swiftest.in\n", + "Reading in time 1.000e+03\n", + "Creating Dataset\n", + "Successfully converted 1001 output frames.\n", + "Swiftest simulation data stored as xarray DataSet .ds\n" + ] + } + ], + "source": [ + "swiftestsim = swiftest.Simulation(param_file=\"param.swiftest.in\")\n", + "swiftestsim.bin2xr()\n", + "swiftestdat = swiftestsim.ds" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "swifterdat['varpi'] = swifterdat['omega'] + swifterdat['capom']\n", + "swiftestdat['varpi'] = swiftestdat['omega'] + swiftestdat['capom']" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "obj = Horizons(id='1', id_type='majorbody',location='@sun',\n", + " epochs={'start':'2021-01-28', 'stop':'3021-02-05',\n", + " 'step':'1y'})\n", + "el = obj.elements()\n", + "t = (el['datetime_jd']-el['datetime_jd'][0]) / 365.25\n", + "varpi_obs = el['w'] + el['Omega']" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "varpiswiftest = swiftestdat['varpi'].sel(id=2) * 180.0 / np.pi\n", + "varpiswifter = swifterdat['varpi'].sel(id=2) * 180.0 / np.pi\n", + "tsim = swiftestdat['time']" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "dvarpi_swiftest = np.diff(varpiswiftest) * 3600 * 100 \n", + "dvarpi_swifter = np.diff(varpiswifter) * 3600 * 100 \n", + "dvarpi_obs = np.diff(varpi_obs) / np.diff(t) * 3600 * 100 " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Mean precession rate for Mercury long. peri. (arcsec/100 y)\n", + "JPL Horizons : 571.3210506300043\n", + "Swifter GR : 571.1981012667947\n", + "Swiftest GR : 1.5844780122245083\n", + "Obs - Swifter : 0.12294936320971743\n", + "Obs - Swiftest : 569.7365726177798\n", + "Swiftest - Swifter: -569.61362325457\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAArP0lEQVR4nO3deXyV5bnv/8+VmTlzyACEeRQjRHCoioh2a1Grp93V7T6V9rSW396n7a4vh7pb3W5ObbU/W92ttW51O/xOu6mtY4/WrW1tHfqzIiCiDDEEApkISSAMCSHTdf5YixgwQBYkayV5vu/XK6+s536m616Bda3nvp/nvs3dERGR4IqLdQAiIhJbSgQiIgGnRCAiEnBKBCIiAadEICIScAmxDuBkZGZmemFhYazDEBEZVNasWVPv7llHlw/KRFBYWMjq1atjHYaIyKBiZtt7KlfTkIhIwCkRiIgEnBKBiEjAKRGIiAScEoGISMApEYiIBJwSgYhIwEXtOQIzmw481a1oEnAH8CfgIWAkUA5c5+77ohWXiMhA5e7sbqijpnI7+yvWE1/zHvmX3UJ+wfg+PU/UEoG7lwBFAGYWD1QBzwFPAze5++tm9mXgZuD2aMUlIhJrB/btprZ8M3u2vUfb7h3E793B6KZtjGrfTb7XkhHers3j2fTRuYM3ERzlIqDM3beHrxTeCJf/HngFJQIRGUIOtTSzb0899RWl7K/aSEfDVuIO7iZnzxpGdu4nk0ZGhrftdKPe0qhPyqd2xAyqcz7PiKwJjBl/GlnjZzB3ZFqfxxerRHANsDL8+kPgCuAF4PPAuJ52MLMbgBsAxo/v22woInKqOtrb2bOrksa6Khp2bCSu6l06Ww+Svncjee2VZFkLhwf5afc4mi2F8pSZ1CbPoTR9Mgmp+WTPPJdRqVlkZ+eRHcXYLdpTVZpZElANzHb3WjObAfwEyAB+C3zD3TOOd4zi4mLXWEMiEnXu7GmopXTT+yTu205HXSnJez5iZHMlaW07SeVA16YtnshBS2FnUiH7RxbC2NOIH57O2Jlnk5FbSMrwkcc+Tz8xszXuXnx0eSyuCC4F1rp7LYC7bwYuATCzacBnYhCTiAgAba2HqK/eRmPNVg7sLKN9TwUJe7eRdqCM7I4a0mhmQXjbDjeqLIe6pALq06bTObaIhFFZZBXOInfyaaQlD6fvG3L6XiwSwbV83CyEmWW7+y4ziwO+S+gOIhGRfuOdndTXVrK7agt7y9fR2bCVpH3lpB2sYHzHdnLNye22/S7S2ZkykY1pp0NaIakFM0nInEz2+KmMHzWawd5YHdVEYGbDgYuBr3UrvtbM/jH8+lng8WjGJCJDk3d2UvbhX2naXcPBqg0k7NtByp6PSG7fz6SObWSZd7XZt3oCtXHZ7EnOZ1XGhcSlT2RYViFp+VPIzp9IdsrwqLbZR1vU+wj6gvoIROSw8k2raawqpbm2lNb9DWTuXkdqWy05nbtItI6u7do8nq2JU7CEZPakzYUxBSSmjyNz4lzGTZ6DxQ3952sHUh+BiEivdR5qpn7HJnbt2klT9WasvoQR+8rIailndOdeCq3tiO3L4idTN3wSlamXEJ8xmZH5M0nLm0R23kSmx8fHqBYDmxKBiMRcc/MBqkrWsr++ggM1pdBQxujWXWQe2kFeRw3Z5l1NM82eTGX8OCqGz6RtRD6WOYXUSfPILZzF6NGjmZw0IqZ1GYyUCEQkKvY3NlC7o4R9FRtg/05aa0sYdaCckW31jPU6plr7x9sygrq4LHamTGbrmL/BsqeTmpFL5oQZZBdMYZq+2fcpJQIR6TON9bU072ugeuNfONRYQ8fuHeTs/5DMthrG+D6mWGfXtnsZQU3iBPaMnknD6CXYuAWMyBpP/sRZjErNYpQZk2JYlyBRIhCRiFRv20zzvgYaKzbQtrOEzsYKRrdUkd1eTQ67SQXywtu2kER5wkS2pJ1H5/Asksedzpj8mSRnjCMvZyxj4od+B+1goEQgIkfwzk721Newq3wDB3aVc6huKykNmxjRUkta+y7yaOjattONhrh0DiaMoTq1mO3p0/Bh6WRMXUD2hOmMTs1ihlkMayO9oUQgEkCdHR3sqtlOw47NHNq1hc6GrSTsLSftwBZGde4lnX2kd9u+2nLYk5jNjjHz2ZY9h6T08aSOm0lqVgFZOQUAg/6hqiBTIhAZojra29m5o5TdFZto3lmK795G8v7tpLZUkttRw1hrY2x42zaPp8ay2D2skNrk0+jMmklS5kTSx81gTGYeedn5Xc09MvQoEYgMYt7Zya7qbVRvfJvWvTvprN9CclMVww41UNC2lXwOkh/etsUT2RmfS2NKAbtGn4ulT2L42KmMzJ1GzrgpjE9J1rf6gFIiEBngvLOTloNNVJW8y77ydbTWlTG6fh1pbbWkeyM51kZOeNsWT6Q+LpP9CelsSl8C+cWMzJtO1vgZZOaOpzBOt13KJykRiMSYu7OvsZ62g03UV3xE45a/wr5qUvaXM7y1gcK2rQyzdqaEt2/1BHYkTKByzBlsH5YFo3NJn3EeozJyySmYQoHusZcIKRGIRElHexu1O0rZU13KnqotJNRtYPiBcsa01jKhsxKAzPC2rR5PTVwuBxIzWJvz3/CUVBLz55I5eT75hVOZkqD/utJ39K9JpA913Xq57UP2V20OTVyybyvpB3eQ21FDnrV3dbo2ezLVCQU0JudRnX05ccPTGJY9iYwJs8ifNJsJMa2JBIkSgUiEDjbtp7GmjJ1bP6Szci1JByoY1lRJ2qFqRngT6dbWdetlqydQHZ9LQ8oEqscsIj5rKsNzJpOVP5mc8dOZEoARL2XgUyIQ6YF3dtBQs4Oq0rUc3LGOlF3vMfrQTsa01pLG/q6JSzrc2EkGdQm51I0+C0seBWmFDMudQeaEWeSMm0phYiKFsa6QyHEoEUhgtbe1sWPrJurK3oN91cRXryax7QBjDlUxtmMnmdbW1WZfRQ47Ewuoy5iJjc6lM3USo7LGMb5oEdkpw8jXUAkyiCkRyJDW2LiHnVvX07j5zdCH/cE6Upu2kdNeRaK3M8lauwY2a2AM++LS2JMyjppRnwrNUpUzhbwZC8jLzidfzTgyRCkRyKC3t6GWmqrt7N+xHna+D/t2En+wgbFt28mjntTwdoc8kd2Wyt7ETDalX4wnDicpdxYZE4tITs8jO28iGbr1UgIoaonAzKYDT3UrmgTcAfyZ0IT1KUA78A/uvipacckg4M7e+ip2bV1P465K2nd9ROLebYxoriKzvYYsGhkT3rTV46mPy+BQwmh2pZ5BddpUEnNnkjVtIbnjp5IbF3fEpOQiEsVE4O4lQBGAmcUDVcBzwCPAv7r7y2Z2GfBDYFG04pKBoaOjk/rGRho2/JlDtSX47nIS91cw6mAlOR07GcOhrg/7Djd2xmWzN2ksO0adw9aM6QzPHMeYcbPInlxEXsqwmNZFZLCJVdPQRUCZu283MwdGh8vHANUxikn6WUd7Owf276Hmo7Xs31lGR80H0NFG+p73KWwrI8c6uoZKaPJkauLGUp+cR1X6WVhaIcNyppKRN4Gc8VPJH5HWNYaOiJyaWCWCa4CV4df/BLxiZvcCccA5Pe1gZjcANwCMH6+hsQaqto5OKnfuoqZ8E4l1m+ioeZ/4lt2Maqogv62cMXbwiGacTuLYmjyDtXnXkDwyncTc2aRNO5fMnDymJKoLSyQazN2je0KzJELf+me7e62Z/QR43d2fMbO/BW5w9yXHO0ZxcbGvXr06GuFKD9raO6ip2k59xWaaakrpbNhG4r7tjGmpYmxHDRm2r2vbQ57IXhtFXVI+B0ZOpDN9KsmZE0gtmMGEGfOJi4vDdDeOSFSY2Rp3Lz66PBZfuS4F1rp7bXj5euCb4de/AR6NQUxylNZDLdRWlLK7ooTm2lK8YRspB3aQ2lLF2M5axtuhriGLO9yoj8+iMTmfylEXUpU+kRFjpxKXNZXx0+eTnZBAdkxrIyLHE4tEcC0fNwtB6OrgAkJ3Dy0GSmMQUyC1t7WyfdNqDjXvY9/WNdjuMlKaKhjZWseE9u2Ms07Ghbc96EnUJoxl3/ACNo4+h7j0SYzMnUrmuBmk5U8mJyG5q31fRAaXqCYCMxsOXAx8rVvxV4F/M7MEoIVwP4D0jfb2DrZ88Ff2NdTQuXsrvqeCEXs2knGogqzOBiZbe9e2zZ5MTXwezYlpvDv2AuIzpzAqbyrZ42eSnjOOQjXhiAxJUU0E7t4MZBxV9hYwP5pxDDVtrYeo3LKexsoSDtV+BE31jNizifxDW0jwdmZY88fbejzbEwrZNWwyFaOXkJB/OgnDx5A/8yzSswuYrOGNRQJH/+sHiX2NDTRUb+VAQzVNFR8QV7eJUQe2kdNWwShvYqJ1dG17eOKSsjFn0ZkwnKRx80kfN50RuVMZNiqdKaNSY1cRERlwlAgGkKb9jdRs3cC+6o9ord9GfMNHDDtYQ/ahHWSzu+thC4A9jKYmcTxb0s6nY1gGibmzGZZeQMGshSSnDGPKsBExq4eIDC5KBFHk7uzZ3UBDxSYat66lfU8F8furGNVUTlpbHWOp65qOEEKDoNUl5LF9TDFlGTNIyphA8ugscqfNJyOngLSY1UREhhIlgj7mnZ007KqkqbGe+q3v01qzgYSD9WTv+4DU1p2kc6Br0pJONxoslfrEPKpHzmRb+lUk585idN40ssZPJyMt88gOFRGRfqBEcBK8s5Nd1dvYu6uC/TvLaCt/h8TmWka21DC6fTe51JEJTADaPY4DDKcyYRyVaUtIyphAQuZEsqacSfa4KWQlp5AV6wqJSKApERzD4blnq0tWc3B3FS27q0mo30Tewc2M6DxADo1d9823egK1cdnsS8ykZuQstucWE5cyivTJxWQXziQ1NbNrKGQRkYEm8ImgsX4nO7d+wP7qj2iv30LS3nJGN+8gp6OadJq7mnEA9tgYKobNoC5hFFty55GUmkfahDlk5E1iXKoacURkcApUIlj3h5W0lPyR5KZqhrXtIaetkjT2dX1b73CjNi6b+uQCNqXNxdMnMTx3JmnjppOelU/amHR10IrIkBOoRNDy0Z84bdf/oS4+i6aEdErTL6AzYyrDcqeTPm4mOeOnk5ecQl6sAxURiaJAJYL5X/kpCQkPaagEEZFuApUIEpOSYx2CiMiAo6/GIiIBp0QgIhJwSgQiIgGnRCAiEnBKBCIiAadEICIScEoEIiIBF7XnCMxsOvBUt6JJwB3A2cD0cFkq0OjuRdGKS0Qk6KKWCNy9BCgCMLN4oAp4zt3vP7yNmf0I2ButmEREJHZPFl8ElLn79sMFZmbA3wKLYxSTiEggxaqP4Bpg5VFl5wG17l4ag3hERAIr6onAzJKAK4DfHLXqWj6ZHLrvd4OZrTaz1XV1df0ZoohIoMTiiuBSYK271x4uMLME4GqO7Ew+grs/7O7F7l6claXJHUVE+kosEkFP3/yXAJvdvTIG8YiIBFrEicDMRoTv+omYmQ0HLgaePWpVT30GIiISBSe8a8jM4gh9UF8HnAkcApLNrA74HfBwbzt43b0Z+MTkvu6+LIKYRUSkD/XmiuBPwGTgNmCsu49z92xCd/n8FbjbzP6+H2MUEZF+1JvnCJa4e9vRhe6+G3gGeMbMEvs8MhERiYoTXhH0lAROZhsRERmYev1ksZnd2EPxXmCNu6/rs4hERCSqIrlrqBhYDuSHf24AFgGPmNktfR+aiIhEQyRjDWUA89z9AICZ/QvwNHA+sAb4Yd+HJyIi/S2SK4LxQGu35TZggrsfJHRLqYiIDEKRXBH8J/BXM3shvHw5sNLMRgAb+zwyERGJil4nAnf/X2b2O+BTgAHL3X11ePV1/RGciIj0v143DYXnC5gJjAlPJrPLzBb0V2AiIhIdkfQRPEhoWslrw8v7gZ/1eUQiIhJVkfQRLHT3eWb2HoC77wnPLSAiIoNYJFcEbeFRRx3AzLKAzn6JSkREoiaSK4KfAM8BOWZ2F/A54Lv9EpWIDGptbW1UVlbS0tIS61ACKSUlhYKCAhITezcMXCR3Df3SzNYQmnge4LPuvukkYhSRIa6yspJRo0ZRWFhI6D4TiRZ3p6GhgcrKSiZOnNirfXozH0FPYwwBXGpml7r7jyMJUkSGvpaWFiWBGDEzMjIyiGRu995cEYwK/55OaGKa34aXLwfeiChCEQkMJYHYifS9P2EicPd/DR/4VUJjDe0PL98J/CbyEEVEZCA5lbGGWoHC3u5sZtPNbF23n31m9k/hdV83sxIz22BmGrxORPrEyJEjKS8vZ9iwYRQVFTFr1iyWL19OZ2cn5eXlzJkz57j733nnndx7771HlBUWFlJfXx9RHJdddhmNjY2Rhh81kdw19L+BVWb2HKFbSK8Cnuztzu5eAhQBhG9DrQKeM7MLgSuBue5+yMyyI4hJROSEJk+ezLp162hvb2fx4sU8//zzzJs3r9/P6+64O7/73e/6/VynotdXBO5+F/AlYA/QCHzJ3X9wkue9CChz9+3A/wPc7e6HwufZdZLHFBE5roSEBM455xy2bNnSJ8f78Y9/zJw5c5gzZw73338/AOXl5cycOZN/+Id/YN68eVRUVHRdRTz00EMUFRVRVFTExIkTufDCCwFYuXIlp512GnPmzOHWW2/tOv7IkSP5zne+w+mnn85ZZ51FbW0tAL/5zW+YM2cOp59+Oueff/4p16M3dw2ZuzuAu68F1h5vm166BlgZfj0NOC/8bEILcJO7vxvBsURkAPvX/7OBjdX7+vSYs/JG8y+Xz454v+bmZv74xz+yYsWKXu9z33338Ytf/KJrubq6GoA1a9bw+OOP88477+DuLFy4kAsuuIC0tDRKSkp4/PHHefDBB4841vLly1m+fDltbW0sXryYG2+8kerqam699VbWrFlDWloal1xyCc8//zyf/exnaWpq4qyzzuKuu+7illtu4ZFHHuG73/0uK1as4JVXXiE/P79Pmpx6c0Xwp3Ab/vjuhWaWZGaLzexJ4PrenjA8LMUVfNzRnACkAWcBNwO/th66vM3sBjNbbWarI7ktSkSkrKyMoqIizj33XD7zmc9w6aWX9nrfb33rW6xbt67rJy8vD4C33nqLq666ihEjRjBy5Eiuvvpq3nzzTQAmTJjAWWeddcxjfvOb32Tx4sVcfvnlvPvuuyxatIisrCwSEhK47rrreOON0A2ZSUlJLF26FID58+dTXl4OwLnnnsuyZct45JFH6OjoOJm35Ai96SP4G+DLhOYemEioWWgYoSTyKnBfhHMWXwqsdffa8HIl8Gz4imKVmXUCmcARn/bu/jDwMEBxcXEkVx8iEkMn8829rx3uI+hLx2sEGTFixDHXPfHEE2zfvp0HHnjghMdJTEzsuhU0Pj6e9vZ2AB566CHeeecdXnrpJYqKili3bh0ZGRknUw2gF1cE7t7i7g+6+7nABELt+2e4+wR3/+pJTFx/LR83CwE8DywGMLNpQBIQWZe8iEiUnX/++Tz//PM0NzfT1NTEc889x3nnnXfcfdasWcO9997LL37xC+LiQh+/Cxcu5PXXX6e+vp6Ojg5WrlzJBRdccNzjlJWVsXDhQlasWEFmZiYVFRWnVJdI7hrC3duAmpM9mZkNBy4Gvtat+DHgMTP7kNAtqddH2N8gIvIJ7e3tJCcnH3ebkpISCgoKupbvu+8+Pv/5z/fq+PPmzWPZsmUsWBCaluUrX/kKZ5xxRlfzTU8eeOABdu/e3dVJXFxczKOPPsoPfvADLrzwQtydyy67jCuvvPK457755pspLS3F3bnooos4/fTTexXzsdhg/MwtLi721atXn3hDEYmJTZs2MXPmzJjG8P777/PVr36VVatWxTSOWOnpb2Bma9y9+OhtI3mgTERkUHjooYe49tpr+d73vhfrUAaFSKaq/B/9GYiISF9Zvnw5Gzdu5JJLLol1KINCJH0EPzKz64B2YBWw0t039E9YIiISLZE0DTUA3wN+DBwgdL//146/i4iIDHSRXBHsdffXwq//y8z+DXgH+Pe+D0tERKIlottHAczsVkIPfI0B9vd5RCIiElUnc9fQM8AWoAD4ft+GIyLSN+666y5mz57N3LlzKSoq4p133unVfnfccQd/+MMfAHjzzTeZPXs2RUVFvP32230yimhtbS1/93d/x6RJk5g/fz5nn302zz33HAB//vOfGTNmDGeccQYzZszgpptuOuXz9UYkiSDNzMa5+xZ3/3dgKUoEIjIAvf3227z44ousXbuW9evX84c//IFx48b1at8VK1awZMkSAH75y19y0003sW7dOkpKSiJOBIeHhDjM3fnsZz/L+eefz9atW1mzZg2/+tWvqKys7NrmvPPO47333uO9997jxRdf5C9/+UtE5zwZkTQNjQb+bGb1wEYgFTj10Y5ERPpYTU0NmZmZXU8WZ2ZmArBq1Sruvvtunn32WV544QWuueYa9u7dS2dnJ7NmzWLr1q0sW7aMpUuX0tjYyK9//WteeeUVXn31Vf7yl79w8OBB3nrrLW677TaWLl3K17/+dT744APa29u58847ufLKK3niiSd46aWXaGlpoampiddee60rrtdee42kpCSWL1/eVTZhwgS+/vWvf6IOhyfTqaqq6ud3K7JEcCHwIbCQ0PzFDrzUH0GJyBDy8rdh5wd9e8yxp8Gldx9z9SWXXMKKFSuYNm0aS5Ys4Qtf+AIXXHAB8+bN47333gNCzT5z5szh3Xffpb29nYULFx5xjK985Su89dZbLF26lM997nM88cQTrF69umuwuH/+539m8eLFPPbYYzQ2NrJgwYKuK4m3336b9evXk56efsQxN2zY0OsJcfbs2UNpaWmfzDdwIpFMTLPe3Tvd/W13f8Ldn3R3DQ4nIgPOyJEjWbNmDQ8//DBZWVl84Qtf4IknniAhIYEpU6awadMmVq1axY033sgbb7zBm2++ecIB44726quvcvfdd1NUVMSiRYtoaWlhx44dAFx88cWfSAI9+cd//EdOP/10zjzzzK6yN998k7lz5zJ27FiWLl3K2LFjI6v8SYj4riERkYgc55t7f4qPj2fRokUsWrSI0047jSeffJJly5Zx3nnn8fLLL5OYmMiSJUtYtmwZHR0dn5ib+ETcnWeeeYbp06cfUf7OO+8ccxjq2bNn88wzz3Qt/+xnP6O+vp7i4o+H/znvvPN48cUX+eijj/jUpz7FVVddRVFRUUSxRUpjDYnIkFNSUkJpaWnX8rp165gwYQIQGj76/vvv5+yzzyYrK4uGhgY2b97M7NnHnzdh1KhR7N//8R3zn/70p/npT3/aNZ/A4San41m8eDEtLS38/Oc/7yprbm7ucdtp06Zx2223cc8995zwuKcqkrGG/qeZpfVnMCIifeHAgQNcf/31zJo1i7lz57Jx40buvPNOIDT+f21tbVfb+9y5c5k7dy49TIx4hAsvvJCNGzdSVFTEU089xe23305bWxtz585lzpw53H777SeMy8x4/vnnef3115k4cSILFizg+uuvP+aH/fLly3njjTfYtm1bZG9AhHo9DLWZfY/QXMNrCc0h8Eqs5g3QMNQiA9tAGIY66PplGGp3/y4wFfgPYBlQambfN7PJpxauiIjEUkR9BOErgJ3hn3ZCk84/bWY/7IfYREQkCnp915CZfQO4ntB8wo8CN7t7m5nFAaXALf0TooiI9KdeJQIL9aKcDlzt7tu7r3P3TjNb2h/BiYhI/+tVInB3N7Mzjk4C3dZvOtExzGw68FS3oknAHYSGqvgqUBcu/2d3P/WRnUREpFci6SN428zOPPFmPXP3EncvcvciYD7QDDwXXn3f4XVKAiIi0RVJIrgQ+KuZlZnZejP7wMzWn+R5LwLKjnWFISJyqgbSMNSNjY08+OCDx1wf66GpI0kElxJqzlkMXE5oGOrLT/K81wAruy3/z3ByeexYD62Z2Q1mttrMVtfV1fW0iYgIMHCGoT7seIlgIAxNHUkiuP4YPxExsyTgCuA34aKfA5OBIqAG+FFP+7n7w+5e7O7FWVlZkZ5WRAKkp2Go8/LyWLVqFVdffTUAL7zwAsOGDaO1tZWWlhYmTZoEwLJly3j66ad59NFH+fWvf82KFSu49tprueOOO3jqqae6nixuamriy1/+MmeeeSZnnHEGL7zwAhAaYXTBggUUFRUxd+5cSktL+fa3v01ZWRlFRUXcfPPNR8Q6EIamjmTQuaZur1MIXRGcsJO4B5cCa929FuDwbwAzewR48SSOKSID1D2r7mHz7s19eswZ6TO4dcGtx1wfy2GoH3roIb75zW9y3XXX0draSkdHB3fffTcffvgh69at+0SsA2Fo6kieLP5Rt5+7gEVA/kmc81q6NQuZWW63dVcRmvNAROSkxXIY6rPPPpvvf//73HPPPWzfvp1hw4ZFdNxYDE19KsNQDyfUZ9BrZjYcuBj4WrfiH5pZEaGJbsqPWicig9zxvrn3p1gNQz1z5kwWLlzISy+9xKc//WkeffTRrmanngyEoakjGX30g3CH7noz2wCUAP8WycncvdndM9x9b7ey/+7up7n7XHe/wt1rIjmmiMjRYjkM9datW5k0aRLf+MY3uOKKK1i/fv0n9u1uIAxNHUln8eG7hC4HLgHy3P2BPo1GRKQPxHIY6qeeeoo5c+ZQVFTE5s2b+eIXv0hGRgbnnnsuc+bM+URn8UAYmrrXw1APJBqGWmRg0zDUsdcvw1Cb2ZNmltptOc3MHjuVQEVEJPYiaRqa6+6NhxfcfQ9wRp9HJCIiURVJIojr/tSvmaVzancdicgQNhibnYeKSN/7SD7IfwT8/2b2NKFbPf8WuCuis4lIIKSkpNDQ0EBGRsYJO2Glb7k7DQ0NpKSk9HqfSOYj+BOwmtBYQ0ZoboKNJxOoiAxtBQUFVFZWonHBYiMlJYWCgoJebx/JfATPu/t8QB/+InJciYmJTJw4MdZhSC9F0kfw11OZj0BERAamSPoILgSWm1k5oQHojNDFwtz+CExERKIjkkRwab9FISIiMRNJ09AO4Dzg+vDMYg7k9EtUIiISNZEkggeBswkNIw2wH/hZn0ckIiJRFUnT0EJ3n2dm70HoyeLwbGMiIjKIRXJF0GZm8YSahDCzLKCzX6ISEZGoiSQR/AR4Dsgxs7uAt4Af9EtUIiISNb1uGnL3X5rZGuCicNGV7t63E5GKiEjUnTARmNlvjy4K//60meHuV/R9WCIiEi29uSI4G6ggNOH8O3ycCCJiZtOBp7oVTQLucPf7w+tvAv5fIMvd60/mHCIiErneJIKxhCacvxb4O+AlYKW7b4jkRO5eAhQBhDudqwj1OWBm48Ln2BHJMUVE5NSdsLPY3Tvc/b/c/XrgLGAL8Gcz+/opnPcioCz8YBrAfcAthO9IEhGR6OntMNTJwGcIXRUUErqD6NlTOO81hJqaMLMrgCp3f/9445ab2Q3ADQDjx48/hVOLiEh3J5y83syeBOYALwO/cvcPT+mEoYfQqoHZhJ5O/hNwibvvDQ9oV3yiPgJNXi8iErljTV7fmyuC/05otNFpwDe6fWs/PPro6AhjuRRY6+61ZnYaMBE4fDVQAKw1swXuvjPC44qIyEk4YSJw90geOuuNawk3C7n7B0D24RW9vSIQEZG+09cf8sdlZsMJ3R10Kv0LIiLShyIZdO6UuXszkHGc9YXRi0ZERCDKVwQiIjLwKBGIiAScEoGISMApEYiIBJwSgYhIwCkRiIgEnBKBiEjAKRGIiAScEoGISMApEYiIBJwSgYhIwCkRiIgEnBKBiEjAKRGIiAScEoGISMApEYiIBJwSgYhIwCkRiIgEXNSmqjSz6cBT3YomAXcQmrrySqAT2AUsc/fqaMUlIhJ0UUsE7l4CFAGYWTxQBTwH7HH328Pl3yCUHJZHKy4RkaCL6uT13VwElLn79qPKRwAeg3hERAIrVongGmDl4QUzuwv4IrAXuLCnHczsBuAGgPHjx0chRBGRYDD36H4BN7MkoBqY7e61R627DUhx93853jGKi4t99erV/RiliMjQY2Zr3L346PJY3DV0KbD26CQQ9p/Af4tyPCIigRaLRHAtRzYLTe227gpgc9QjEhEJsKj2EZjZcOBi4Gvdiu8O31raCWxHdwyJiERVVBOBuzcTem6ge5magkREYkhPFouIBJwSgYhIwCkRiIgEnBKBiEjAKRGIiAScEoGISMApEYiIBJwSgYhIwCkRiIgEnBKBiEjAKRGIiAScEoGISMApEYiIBJwSgYhIwCkRiIgEnBKBiEjAKRGIiAScEoGISMBFbarK8LzET3UrmgTcAeQDlwOtQBnwJXdvjFZcIiJBF7UrAncvcfcidy8C5gPNwHPA74E57j4X+Ai4LVoxiYhI7JqGLgLK3H27u7/q7u3h8r8CBTGKSUQkkGKVCK4BVvZQ/mXg5Z52MLMbzGy1ma2uq6vr1+BERIIk6onAzJKAK4DfHFX+HaAd+GVP+7n7w+5e7O7FWVlZ/R+oiEhARK2zuJtLgbXuXnu4wMyuB5YCF7m7xyAmEZHAikUiuJZuzUJm9jfArcAF7t4cg3hERAItqk1DZjYcuBh4tlvxA8Ao4Pdmts7MHopmTCIiQRfVK4LwN/6Mo8qmROv8D73/EL/b9ju6tz453rV8+LXTbf1Ry0c77rpjtHIdb5/jr4r8eMdraTtuHH18vD5/n6JlQIQQ+yAGRAwDoNV4ILwP9194P+fkndOnx4xF01DMZA/PZmrqVMwMwwBCv42Pl8PrDMPMuvY9vL4n3bf7xLpj7Hcy+xzPyR7vZOvV18c7mfcpWk7m7zEU6W8xMGIYO2Jsnx8zUIng6qlXc/XUq2MdhojIgKKxhkREAk6JQEQk4JQIREQCTolARCTglAhERAJOiUBEJOCUCEREAk6JQEQk4GwgPLYdKTOrA7af5O6ZQH0fhjMYqM7BoDoHw6nUeYK7f2Ic/0GZCE6Fma129+JYxxFNqnMwqM7B0B91VtOQiEjAKRGIiARcEBPBw7EOIAZU52BQnYOhz+scuD4CERE5UhCvCEREpBslAhGRgAtUIjCzvzGzEjPbYmbfjnU8fcHMxpnZn8xsk5ltMLNvhsvTzez3ZlYa/p3WbZ/bwu9BiZl9OnbRnxozizez98zsxfDykK6zmaWa2dNmtjn89z47AHX+Vvjf9YdmttLMUoZanc3sMTPbZWYfdiuLuI5mNt/MPgiv+4lFMqWcuwfiB4gHyoBJQBLwPjAr1nH1Qb1ygXnh16OAj4BZwA+Bb4fLvw3cE349K1z3ZGBi+D2Jj3U9TrLuNwL/CbwYXh7SdQaeBL4Sfp0EpA7lOgP5wDZgWHj518CyoVZn4HxgHvBht7KI6wisAs4GDHgZuLS3MQTpimABsMXdt7p7K/Ar4MoYx3TK3L3G3deGX+8HNhH6D3QloQ8Owr8/G359JfArdz/k7tuALYTem0HFzAqAzwCPdisesnU2s9GEPjD+A8DdW929kSFc57AEYJiZJQDDgWqGWJ3d/Q1g91HFEdXRzHKB0e7+toeywv/XbZ8TClIiyAcqui1XhsuGDDMrBM4A3gFy3L0GQskCyA5vNlTeh/uBW4DObmVDuc6TgDrg8XBz2KNmNoIhXGd3rwLuBXYANcBed3+VIVznbiKtY3749dHlvRKkRNBTe9mQuXfWzEYCzwD/5O77jrdpD2WD6n0ws6XALndf09tdeigbVHUm9M14HvBzdz8DaCLUZHAsg77O4XbxKwk1geQBI8zs74+3Sw9lg6rOvXCsOp5S3YOUCCqBcd2WCwhdZg56ZpZIKAn80t2fDRfXhi8XCf/eFS4fCu/DucAVZlZOqIlvsZn9gqFd50qg0t3fCS8/TSgxDOU6LwG2uXudu7cBzwLnMLTrfFikdawMvz66vFeClAjeBaaa2UQzSwKuAX4b45hOWfjOgP8ANrn7j7ut+i1wffj19cAL3cqvMbNkM5sITCXUyTRouPtt7l7g7oWE/o6vufvfM7TrvBOoMLPp4aKLgI0M4ToTahI6y8yGh/+dX0SoD2wo1/mwiOoYbj7ab2Znhd+rL3bb58Ri3WMe5d75ywjdVVMGfCfW8fRRnT5F6BJwPbAu/HMZkAH8ESgN/07vts93wu9BCRHcWTAQf4BFfHzX0JCuM1AErA7/rZ8H0gJQ538FNgMfAv+b0N0yQ6rOwEpCfSBthL7Z/4+TqSNQHH6fyoAHCI8c0ZsfDTEhIhJwQWoaEhGRHigRiIgEnBKBiEjAKRGIiAScEoGISMApEUigmVmGma0L/+w0s6rw6wNm9mA/nfOfzOyLJ9jmV2Y2tT/OL3I03T4qEmZmdwIH3P3efjxHArCW0Iix7cfZ7gLg7939q/0Vi8hhuiIQ6YGZLeo2z8GdZvakmb1qZuVmdrWZ/TA89vt/hYf4ODwe/OtmtsbMXjk8RMBRFgNr3b3dzCab2dpu55xqZofHT3oTWBJOHCL9SolApHcmExr2+krgF8Cf3P004CDwmXAy+CnwOXefDzwG3NXDcc4F1gC4exmw18yKwuu+BDwRXtdJaIjh0/upPiJd9G1DpHdedvc2M/uA0CRH/xUu/wAoBKYDc4DfhyeGiic0bMDRcgmNl3PYo8CXzOxG4AscOX7+LkKjbvZ2lFWRk6JEINI7hyD0Td3M2vzjzrVOQv+PDNjg7mef4DgHgZRuy88A/wK8Bqxx94Zu61LC24v0KzUNifSNEiDLzM6G0NDgZja7h+02AVMOL7h7C/AK8HPg8aO2nQZs6J9wRT6mRCDSBzw0/enngHvM7H1Co8Ce08OmLxOacrK7XxIaQfbVwwVmlgMc9PAsVSL9SbePikSZmT0H3OLupeHlm4Ax7n57t22+Bexz9/+IUZgSIOojEIm+bxPqNC4NJ4XJhG4r7a6R0Pj7Iv1OVwQiIgGnPgIRkYBTIhARCTglAhGRgFMiEBEJOCUCEZGA+7+0PkMTQghxhAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "\n", + "ax.plot(t, varpi_obs, label=\"JPL Horizons\")\n", + "ax.plot(tsim, varpiswifter, label=\"Swifter GR\")\n", + "ax.plot(tsim, varpiswiftest, label=\"Swiftest GR\")\n", + "ax.set_xlabel('Time (y)')\n", + "ax.set_ylabel('Mercury $\\\\varpi$ (deg)')\n", + "ax.legend()\n", + "print('Mean precession rate for Mercury long. peri. (arcsec/100 y)')\n", + "print(f'JPL Horizons : {np.mean(dvarpi_obs)}')\n", + "print(f'Swifter GR : {np.mean(dvarpi_swifter)}')\n", + "print(f'Swiftest GR : {np.mean(dvarpi_swiftest)}')\n", + "print(f'Obs - Swifter : {np.mean(dvarpi_obs - dvarpi_swifter)}')\n", + "print(f'Obs - Swiftest : {np.mean(dvarpi_obs - dvarpi_swiftest)}')\n", + "print(f'Swiftest - Swifter: {np.mean(dvarpi_swiftest - dvarpi_swifter)}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "swiftestOOF", + "language": "python", + "name": "swiftestoof" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.10" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/python/swiftest/swiftest/io.py b/python/swiftest/swiftest/io.py index 2b8c47cb3..4fd797b76 100644 --- a/python/swiftest/swiftest/io.py +++ b/python/swiftest/swiftest/io.py @@ -1048,13 +1048,22 @@ def select_active_from_frame(ds, param, framenum=-1): frame = ds.isel(time=[framenum]) iframe = frame.isel(time=0) + if "name" in ds.dims: + count_dim = "name" + elif "id" in ds.dims: + count_dim = "id" + # Select only the active particles at this time step # Remove the inactive particles if param['OUT_FORM'] == 'XV' or param['OUT_FORM'] == 'XVEL': - iactive = iframe['id'].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['xhx'])), drop=True).id + iactive = iframe[count_dim].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['xhx'])), drop=True)[count_dim] else: - iactive = iframe['id'].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['a'])), drop = True).id - frame = frame.sel(id=iactive.values) + iactive = iframe[count_dim].where((~np.isnan(iframe['Gmass'])) | (~np.isnan(iframe['a'])), drop = True)[count_dim] + if count_dim == "id": + frame = frame.sel(id=iactive.values) + elif count_dim == "name": + frame = frame.sel(name=iactive.values) + return frame @@ -1079,6 +1088,10 @@ def swiftest_xr2infile(ds, param, in_type="NETCDF_DOUBLE", infile_name=None,fram param_tmp['OUT_FORM'] = param['IN_FORM'] frame = select_active_from_frame(ds, param_tmp, framenum) + if "name" in frame.dims: + frame = frame.swap_dims({"name" : "id"}) + frame = frame.reset_coords("name") + if in_type == "NETCDF_DOUBLE" or in_type == "NETCDF_FLOAT": # Convert strings back to byte form and save the NetCDF file # Note: xarray will call the character array dimension string32. The Fortran code @@ -1186,8 +1199,13 @@ def swifter_xr2infile(ds, param, framenum=-1): ------- A set of input files for a Swifter run """ + frame = ds.isel(time=framenum) + if "name" in frame.dims: + frame = frame.swap_dims({"name" : "id"}) + frame = frame.reset_coords("name") + cb = frame.where(frame.id == 0, drop=True) pl = frame.where(frame.id > 0, drop=True) pl = pl.where(np.invert(np.isnan(pl['Gmass'])), drop=True).drop_vars(['j2rp2', 'j4rp4']) diff --git a/python/swiftest/swiftest/simulation_class.py b/python/swiftest/swiftest/simulation_class.py index 7c17c3d32..c6fe5e8da 100644 --- a/python/swiftest/swiftest/simulation_class.py +++ b/python/swiftest/swiftest/simulation_class.py @@ -68,6 +68,8 @@ def __init__(self,read_param: bool = True, **kwargs: Any): integrator : {"symba","rmvs","whm","helio"}, default "symba" Name of the n-body integrator that will be used when executing a run. Parameter input file equivalent: None + read_param : bool, default False + Read the parameter file given by `param_file`. param_file : str, path-like, or file-lke, default "param.in" Name of the parameter input file that will be passed to the integrator. Parameter input file equivalent: None @@ -291,7 +293,7 @@ def __init__(self,read_param: bool = True, **kwargs: Any): # Set the location of the parameter input file param_file = kwargs.pop("param_file",self.param_file) - read_param = kwargs.pop("read_param",True) + read_param = kwargs.pop("read_param",False) self.set_parameter(verbose=False,param_file=param_file) #----------------------------------------------------------------- @@ -301,9 +303,8 @@ def __init__(self,read_param: bool = True, **kwargs: Any): # If the user asks to read in an old parameter file, override any default parameters with values from the file # If the file doesn't exist, flag it for now so we know to create it if read_param: - if os.path.exists(self.param_file): - self.read_param(self.param_file, codename=self.codename, verbose=self.verbose) - param_file_found = True + #good_param is self.read_param() + if self.read_param(): # We will add the parameter file to the kwarg list. This will keep the set_parameter method from # overriding everything with defaults when there are no arguments passed to Simulation() kwargs['param_file'] = self.param_file @@ -743,7 +744,7 @@ def get_parameter(self, **kwargs): return param_dict def set_integrator(self, - codename: Literal["swiftest", "swifter", "swift"] | None = None, + codename: Literal["Swiftest", "Swifter", "Swift"] | None = None, integrator: Literal["symba","rmvs","whm","helio"] | None = None, mtiny: float | None = None, gmtiny: float | None = None, @@ -2316,58 +2317,82 @@ def _combine_and_fix_dsnew(self,dsnew): Updated Dataset with ntp, npl values and types fixed. """ + if "id" not in self.data.dims: + if len(np.unique(dsnew['name'])) == len(dsnew['name']): + dsnew = dsnew.swap_dims({"id" : "name"}) + dsnew = dsnew.reset_coords("id") + else: + warnings.warn("Non-unique names detected for bodies. The Dataset will be dimensioned by integer id instead of name.") + warnings.warn("Consider using unique names instead.") + + if self.param['OUT_TYPE'] == "NETCDF_DOUBLE": + dsnew = io.fix_types(dsnew, ftype=np.float64) + elif self.param['OUT_TYPE'] == "NETCDF_FLOAT": + dsnew = io.fix_types(dsnew, ftype=np.float32) self.data = xr.combine_by_coords([self.data, dsnew]) def get_nvals(ds): + if "name" in ds.dims: + count_dim = "name" + elif "id" in ds.dims: + count_dim = "id" if "Gmass" in ds: - ds['ntp'] = ds['id'].where(np.isnan(ds['Gmass'])).count(dim="id") - ds['npl'] = ds['id'].where(~(np.isnan(ds['Gmass']))).count(dim="id") - 1 + ds['ntp'] = ds[count_dim].where(np.isnan(ds['Gmass'])).count(dim=count_dim) + ds['npl'] = ds[count_dim].where(~(np.isnan(ds['Gmass']))).count(dim=count_dim) - 1 else: - ds['ntp'] = ds['id'].count(dim="id") + ds['ntp'] = ds[count_dim].count(dim=count_dim) ds['npl'] = xr.full_like(ds['ntp'],0) return ds dsnew = get_nvals(dsnew) self.data = get_nvals(self.data) - if self.param['OUT_TYPE'] == "NETCDF_DOUBLE": - dsnew = io.fix_types(dsnew, ftype=np.float64) - self.data = io.fix_types(self.data, ftype=np.float64) - elif self.param['OUT_TYPE'] == "NETCDF_FLOAT": - dsnew = io.fix_types(dsnew, ftype=np.float32) - self.data = io.fix_types(self.data, ftype=np.float32) - return dsnew - def read_param(self, param_file, codename="Swiftest", verbose=True): + def read_param(self, + param_file : os.PathLike | str | None = None, + codename: Literal["Swiftest", "Swifter", "Swift"] | None = None, + verbose: bool | None = None): """ Reads in an input parameter file and stores the values in the param dictionary. Parameters ---------- - param_file : string + param_file : str or path-like, default is the value of the Simulation object's internal `param_file`. File name of the input parameter file - codename : string + codename : {"Swiftest", "Swifter", "Swift"}, default is the value of the Simulation object's internal`codename` Type of parameter file, either "Swift", "Swifter", or "Swiftest" + verbose : bool, default is the value of the Simulation object's internal `verbose` + If set to True, then more information is printed by Simulation methods as they are executed. Setting to + False suppresses most messages other than errors. Returns ------- - + True if the parameter file exists and is read correctly. False otherwise. """ + if param_file is None: + param_file = self.param_file + + if coename is None: + codename = self.codename + + if verbose is None: + verbose = self.verbose + + if not os.path.exists(param_file): + return False + if codename == "Swiftest": - param_old = self.param.copy() - self.param = io.read_swiftest_param(param_file, param_old, verbose=verbose) - self.codename = "Swiftest" + self.param = io.read_swiftest_param(param_file, param, verbose=verbose) elif codename == "Swifter": self.param = io.read_swifter_param(param_file, verbose=verbose) - self.codename = "Swifter" elif codename == "Swift": self.param = io.read_swift_param(param_file, verbose=verbose) - self.codename = "Swift" else: warnings.warn(f'{codename} is not a recognized code name. Valid options are "Swiftest", "Swifter", or "Swift".') - self.codename = "Unknown" - return + return False + + return True def write_param(self, codename: Literal["Swiftest", "Swifter", "Swift"] | None = None, From be6893487d307a92c876af8e3ed2ffd85d6289a9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 20:26:51 -0500 Subject: [PATCH 08/13] Added the GR/no GR comparison test Notebook to the whm_gr_test example --- .../whm_gr_test/swiftest_relativity.ipynb | 997 ++++++++++++++++-- 1 file changed, 927 insertions(+), 70 deletions(-) diff --git a/examples/whm_gr_test/swiftest_relativity.ipynb b/examples/whm_gr_test/swiftest_relativity.ipynb index ae586907c..0cd73753f 100644 --- a/examples/whm_gr_test/swiftest_relativity.ipynb +++ b/examples/whm_gr_test/swiftest_relativity.ipynb @@ -6,11 +6,11 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "import pandas as pd\n", "import swiftest\n", - "from astroquery.jplhorizons import Horizons" + "from astroquery.jplhorizons import Horizons\n", + "import datetime\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" ] }, { @@ -22,18 +22,442 @@ "name": "stdout", "output_type": "stream", "text": [ - "Reading Swifter file param.swifter.in\n", - "Reading in time 1.000e+03\n", - "Creating Dataset\n", - "Successfully converted 1001 output frames.\n", - "Swifter simulation data stored as xarray DataSet .ds\n" + "Creating the Sun as a central body\n", + "Fetching ephemerides data for Mercury from JPL/Horizons\n", + "Fetching ephemerides data for Venus from JPL/Horizons\n", + "Fetching ephemerides data for Earth from JPL/Horizons\n", + "Fetching ephemerides data for Mars from JPL/Horizons\n", + "Fetching ephemerides data for Jupiter from JPL/Horizons\n", + "Fetching ephemerides data for Saturn from JPL/Horizons\n", + "Fetching ephemerides data for Uranus from JPL/Horizons\n", + "Fetching ephemerides data for Neptune from JPL/Horizons\n", + "Writing initial conditions to file init_cond.nc\n", + "Writing parameter inputs to file /home/daminton/git_debug/swiftest/examples/whm_gr_test/param.gr.in\n" ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:        (name: 9, time: 1)\n",
+       "Coordinates:\n",
+       "  * name           (name) <U32 'Sun' 'Mercury' 'Venus' ... 'Uranus' 'Neptune'\n",
+       "  * time           (time) float64 0.0\n",
+       "Data variables: (12/14)\n",
+       "    particle_type  (name) <U32 'Central Body' 'Massive Body' ... 'Massive Body'\n",
+       "    id             (name) int64 0 1 2 3 4 5 6 7 8\n",
+       "    Gmass          (time, name) float64 39.48 6.554e-06 ... 0.001724 0.002034\n",
+       "    radius         (time, name) float64 0.00465 1.631e-05 ... 0.0001646\n",
+       "    j2rp2          (time, name) float64 4.754e-12 nan nan nan ... nan nan nan\n",
+       "    j4rp4          (time, name) float64 -2.247e-18 nan nan nan ... nan nan nan\n",
+       "    ...             ...\n",
+       "    xhz            (time, name) float64 nan -0.002529 ... -0.0407 -0.7287\n",
+       "    vhx            (time, name) float64 nan -9.241 3.454 ... -1.315 -0.08755\n",
+       "    vhy            (time, name) float64 nan 7.755 6.477 ... 1.932 0.5372 1.149\n",
+       "    vhz            (time, name) float64 nan 1.481 -0.1103 ... 0.01905 -0.02168\n",
+       "    ntp            (time) int64 0\n",
+       "    npl            (time) int64 8
" + ], + "text/plain": [ + "\n", + "Dimensions: (name: 9, time: 1)\n", + "Coordinates:\n", + " * name (name) \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:        (name: 9, time: 1)\n",
+       "Coordinates:\n",
+       "  * name           (name) <U32 'Sun' 'Mercury' 'Venus' ... 'Uranus' 'Neptune'\n",
+       "  * time           (time) float64 0.0\n",
+       "Data variables: (12/14)\n",
+       "    particle_type  (name) <U32 'Central Body' 'Massive Body' ... 'Massive Body'\n",
+       "    id             (name) int64 0 1 2 3 4 5 6 7 8\n",
+       "    Gmass          (time, name) float64 39.48 6.554e-06 ... 0.001724 0.002034\n",
+       "    radius         (time, name) float64 0.00465 1.631e-05 ... 0.0001646\n",
+       "    j2rp2          (time, name) float64 4.754e-12 nan nan nan ... nan nan nan\n",
+       "    j4rp4          (time, name) float64 -2.247e-18 nan nan nan ... nan nan nan\n",
+       "    ...             ...\n",
+       "    xhz            (time, name) float64 nan -0.002529 ... -0.0407 -0.7287\n",
+       "    vhx            (time, name) float64 nan -9.241 3.454 ... -1.315 -0.08755\n",
+       "    vhy            (time, name) float64 nan 7.755 6.477 ... 1.932 0.5372 1.149\n",
+       "    vhz            (time, name) float64 nan 1.481 -0.1103 ... 0.01905 -0.02168\n",
+       "    ntp            (time) int64 0\n",
+       "    npl            (time) int64 8
" + ], + "text/plain": [ + "\n", + "Dimensions: (name: 9, time: 1)\n", + "Coordinates:\n", + " * name (name) " - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], + "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "\n", "ax.plot(t, varpi_obs, label=\"JPL Horizons\")\n", - "ax.plot(tsim, varpiswifter, label=\"Swifter GR\")\n", - "ax.plot(tsim, varpiswiftest, label=\"Swiftest GR\")\n", + "ax.plot(tsim, varpisim_gr, label=\"Swiftest WHM GR\")\n", + "ax.plot(tsim, varpisim_nogr, label=\"Swiftest WHM No GR\")\n", "ax.set_xlabel('Time (y)')\n", "ax.set_ylabel('Mercury $\\\\varpi$ (deg)')\n", "ax.legend()\n", "print('Mean precession rate for Mercury long. peri. (arcsec/100 y)')\n", - "print(f'JPL Horizons : {np.mean(dvarpi_obs)}')\n", - "print(f'Swifter GR : {np.mean(dvarpi_swifter)}')\n", - "print(f'Swiftest GR : {np.mean(dvarpi_swiftest)}')\n", - "print(f'Obs - Swifter : {np.mean(dvarpi_obs - dvarpi_swifter)}')\n", - "print(f'Obs - Swiftest : {np.mean(dvarpi_obs - dvarpi_swiftest)}')\n", - "print(f'Swiftest - Swifter: {np.mean(dvarpi_swiftest - dvarpi_swifter)}')" + "print(f'JPL Horizons : {np.mean(dvarpi_obs)}')\n", + "print(f'Swiftest No GR : {np.mean(dvarpi_nogr)}')\n", + "print(f'Swiftest GR : {np.mean(dvarpi_gr)}')\n", + "print(f'Obs - Swiftest GR : {np.mean(dvarpi_obs - dvarpi_gr)}')\n", + "print(f'Obs - Swiftest No GR : {np.mean(dvarpi_obs - dvarpi_nogr)}')" ] }, { @@ -171,9 +1028,9 @@ ], "metadata": { "kernelspec": { - "display_name": "swiftestOOF", + "display_name": "swiftest", "language": "python", - "name": "swiftestoof" + "name": "swiftest" }, "language_info": { "codemirror_mode": { @@ -185,7 +1042,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.10" + "version": "3.8.5" } }, "nbformat": 4, From 9eb46d7652683b2732a523a412a49ba3560094ef Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 22:50:00 -0500 Subject: [PATCH 09/13] Fixed problem where orbital element input files were being read in as radians instead of degrees. Converted from degrees to radians on input. --- src/netcdf/netcdf.f90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index 16333dec8..aa599f94e 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -626,24 +626,28 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr) if (ntp > 0) tp%e(:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%inc_varid, rtemp, start=[1, tslot]), "netcdf_read_frame_system nf90_getvar inc_varid" ) + rtemp = rtemp * DEG2RAD if (.not.allocated(pl%inc)) allocate(pl%inc(npl)) if (.not.allocated(tp%inc)) allocate(tp%inc(ntp)) if (npl > 0) pl%inc(:) = pack(rtemp, plmask) if (ntp > 0) tp%inc(:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%capom_varid, rtemp, start=[1, tslot]), "netcdf_read_frame_system nf90_getvar capom_varid" ) + rtemp = rtemp * DEG2RAD if (.not.allocated(pl%capom)) allocate(pl%capom(npl)) if (.not.allocated(tp%capom)) allocate(tp%capom(ntp)) if (npl > 0) pl%capom(:) = pack(rtemp, plmask) if (ntp > 0) tp%capom(:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%omega_varid, rtemp, start=[1, tslot]), "netcdf_read_frame_system nf90_getvar omega_varid" ) + rtemp = rtemp * DEG2RAD if (.not.allocated(pl%omega)) allocate(pl%omega(npl)) if (.not.allocated(tp%omega)) allocate(tp%omega(ntp)) if (npl > 0) pl%omega(:) = pack(rtemp, plmask) if (ntp > 0) tp%omega(:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%capm_varid, rtemp, start=[1, tslot]), "netcdf_read_frame_system nf90_getvar capm_varid" ) + rtemp = rtemp * DEG2RAD if (.not.allocated(pl%capm)) allocate(pl%capm(npl)) if (.not.allocated(tp%capm)) allocate(tp%capm(ntp)) if (npl > 0) pl%capm(:) = pack(rtemp, plmask) From ecbaa9301be95b7d27de6b617a1c476a8ff2f94b Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 22:50:27 -0500 Subject: [PATCH 10/13] Updated the GR test Notebook --- .../whm_gr_test/swiftest_relativity.ipynb | 941 +----------------- 1 file changed, 34 insertions(+), 907 deletions(-) diff --git a/examples/whm_gr_test/swiftest_relativity.ipynb b/examples/whm_gr_test/swiftest_relativity.ipynb index 0cd73753f..6abb9524a 100644 --- a/examples/whm_gr_test/swiftest_relativity.ipynb +++ b/examples/whm_gr_test/swiftest_relativity.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -15,913 +15,48 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Creating the Sun as a central body\n", - "Fetching ephemerides data for Mercury from JPL/Horizons\n", - "Fetching ephemerides data for Venus from JPL/Horizons\n", - "Fetching ephemerides data for Earth from JPL/Horizons\n", - "Fetching ephemerides data for Mars from JPL/Horizons\n", - "Fetching ephemerides data for Jupiter from JPL/Horizons\n", - "Fetching ephemerides data for Saturn from JPL/Horizons\n", - "Fetching ephemerides data for Uranus from JPL/Horizons\n", - "Fetching ephemerides data for Neptune from JPL/Horizons\n", - "Writing initial conditions to file init_cond.nc\n", - "Writing parameter inputs to file /home/daminton/git_debug/swiftest/examples/whm_gr_test/param.gr.in\n" - ] - }, - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "
<xarray.Dataset>\n",
-       "Dimensions:        (name: 9, time: 1)\n",
-       "Coordinates:\n",
-       "  * name           (name) <U32 'Sun' 'Mercury' 'Venus' ... 'Uranus' 'Neptune'\n",
-       "  * time           (time) float64 0.0\n",
-       "Data variables: (12/14)\n",
-       "    particle_type  (name) <U32 'Central Body' 'Massive Body' ... 'Massive Body'\n",
-       "    id             (name) int64 0 1 2 3 4 5 6 7 8\n",
-       "    Gmass          (time, name) float64 39.48 6.554e-06 ... 0.001724 0.002034\n",
-       "    radius         (time, name) float64 0.00465 1.631e-05 ... 0.0001646\n",
-       "    j2rp2          (time, name) float64 4.754e-12 nan nan nan ... nan nan nan\n",
-       "    j4rp4          (time, name) float64 -2.247e-18 nan nan nan ... nan nan nan\n",
-       "    ...             ...\n",
-       "    xhz            (time, name) float64 nan -0.002529 ... -0.0407 -0.7287\n",
-       "    vhx            (time, name) float64 nan -9.241 3.454 ... -1.315 -0.08755\n",
-       "    vhy            (time, name) float64 nan 7.755 6.477 ... 1.932 0.5372 1.149\n",
-       "    vhz            (time, name) float64 nan 1.481 -0.1103 ... 0.01905 -0.02168\n",
-       "    ntp            (time) int64 0\n",
-       "    npl            (time) int64 8
" - ], - "text/plain": [ - "\n", - "Dimensions: (name: 9, time: 1)\n", - "Coordinates:\n", - " * name (name) \n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "
<xarray.Dataset>\n",
-       "Dimensions:        (name: 9, time: 1)\n",
-       "Coordinates:\n",
-       "  * name           (name) <U32 'Sun' 'Mercury' 'Venus' ... 'Uranus' 'Neptune'\n",
-       "  * time           (time) float64 0.0\n",
-       "Data variables: (12/14)\n",
-       "    particle_type  (name) <U32 'Central Body' 'Massive Body' ... 'Massive Body'\n",
-       "    id             (name) int64 0 1 2 3 4 5 6 7 8\n",
-       "    Gmass          (time, name) float64 39.48 6.554e-06 ... 0.001724 0.002034\n",
-       "    radius         (time, name) float64 0.00465 1.631e-05 ... 0.0001646\n",
-       "    j2rp2          (time, name) float64 4.754e-12 nan nan nan ... nan nan nan\n",
-       "    j4rp4          (time, name) float64 -2.247e-18 nan nan nan ... nan nan nan\n",
-       "    ...             ...\n",
-       "    xhz            (time, name) float64 nan -0.002529 ... -0.0407 -0.7287\n",
-       "    vhx            (time, name) float64 nan -9.241 3.454 ... -1.315 -0.08755\n",
-       "    vhy            (time, name) float64 nan 7.755 6.477 ... 1.932 0.5372 1.149\n",
-       "    vhz            (time, name) float64 nan 1.481 -0.1103 ... 0.01905 -0.02168\n",
-       "    ntp            (time) int64 0\n",
-       "    npl            (time) int64 8
" - ], - "text/plain": [ - "\n", - "Dimensions: (name: 9, time: 1)\n", - "Coordinates:\n", - " * name (name) Date: Thu, 17 Nov 2022 22:55:46 -0500 Subject: [PATCH 11/13] Converted the swiftest_relativity Notebook to a Python script. --- examples/whm_gr_test/init_cond.py | 226 -------------------- examples/whm_gr_test/swiftest_relativity.py | 60 ++++++ 2 files changed, 60 insertions(+), 226 deletions(-) delete mode 100644 examples/whm_gr_test/init_cond.py create mode 100644 examples/whm_gr_test/swiftest_relativity.py diff --git a/examples/whm_gr_test/init_cond.py b/examples/whm_gr_test/init_cond.py deleted file mode 100644 index 7904eb100..000000000 --- a/examples/whm_gr_test/init_cond.py +++ /dev/null @@ -1,226 +0,0 @@ -import numpy as np -import sys -from astroquery.jplhorizons import Horizons -import astropy.constants as const - -#Values from JPL Horizons -AU2M = const.au.value -GMSunSI = const.GM_sun.value -Rsun = const.R_sun.value -GC = const.G.value -JD = 86400 -year = 365.25 * JD -c = 299792458.0 -MSun_over_Mpl = [6023600.0, - 408523.71, - 328900.56, - 3098708., - 1047.3486, - 3497.898, - 22902.98, - 19412.24, - 1.35e8] - -MU2KG = GMSunSI / GC #Conversion from mass unit to kg -DU2M = AU2M #Conversion from radius unit to centimeters -TU2S = year #Conversion from time unit to seconds -GU = GC / (DU2M**3 / (MU2KG * TU2S**2)) - -GMSun = GMSunSI / (DU2M**3 / TU2S**2) - -t_print = 10.e0 * year / TU2S #output interval to print results -deltaT = 0.25 * JD / TU2S #timestep simulation -end_sim = 1.0e3 * year / TU2S + t_print #end time - -# Solar oblatenes values: From Mecheri et al. (2004), using Corbard (b) 2002 values (Table II) -J2 = 2.198e-7 * (Rsun / DU2M)**2 -J4 = -4.805e-9 * (Rsun / DU2M)**4 - -tstart = '2021-01-28' -tend = '2021-01-29' -tstep = '1d' -planetid = { - 'mercury' : '1', - 'venus' : '2', - 'earthmoon' : '3', - 'mars' : '4', - 'jupiter' : '5', - 'saturn' : '6', - 'uranus' : '7', - 'neptune' : '8', - 'plutocharon' : '9' -} -npl = 9 - -#Planet Msun/M ratio -MSun_over_Mpl = { - 'mercury' : 6023600.0, - 'venus' : 408523.71, - 'earthmoon' : 328900.56, - 'mars' : 3098708., - 'jupiter' : 1047.3486, - 'saturn' : 3497.898, - 'uranus' : 22902.98, - 'neptune' : 19412.24, - 'plutocharon' : 1.35e8 -} - -#Planet radii in meters -Rpl = { - 'mercury' : 2439.4e3, - 'venus' : 6051.8e3, - 'earthmoon' : 6371.0084e3, # Earth only for radius - 'mars' : 3389.50e3, - 'jupiter' : 69911e3, - 'saturn' : 58232.0e3, - 'uranus' : 25362.e3, - 'neptune' : 24622.e3, - 'plutocharon' : 1188.3e3 -} - -pdata = {} -plvec = {} -Rhill = {} - -for key,val in planetid.items(): - pdata[key] = Horizons(id=val, id_type='majorbody',location='@sun', - epochs={'start': tstart, 'stop': tend, - 'step': tstep}) - plvec[key] = np.array([pdata[key].vectors()['x'][0], - pdata[key].vectors()['y'][0], - pdata[key].vectors()['z'][0], - pdata[key].vectors()['vx'][0], - pdata[key].vectors()['vy'][0], - pdata[key].vectors()['vz'][0] - ]) - Rhill[key] = pdata[key].elements()['a'][0] * (3 * MSun_over_Mpl[key])**(-1.0 / 3.0) - - -if __name__ == '__main__': - # Convert from AU-day to AU-year just because I find it easier to keep track of the sim progress - for plid in plvec: - plvec[plid][3:] *= year / JD - - # Names of all output files - swifter_input = "param.swifter.in" - swifter_pl = "pl.swifter.in" - swifter_tp = "tp.swifter.in" - swifter_bin = "bin.swifter.dat" - swifter_enc = "enc.swifter.dat" - - swiftest_input = "param.swiftest.in" - swiftest_pl = "pl.swiftest.in" - swiftest_tp = "tp.swiftest.in" - swiftest_cb = "cb.swiftest.in" - swiftest_bin = "bin.swiftest.dat" - swiftest_enc = "enc.swiftest.dat" - - # Simulation start, stop, and output cadence times - t_0 = 0 # simulation start time - end_sim = 1000.0e0 * year / TU2S # simulation end time - deltaT = 0.25 * JD / TU2S # simulation step size - t_print = 1.0 * year / TU2S #output interval to print results - - iout = int(np.ceil(t_print / deltaT)) - rmin = Rsun / DU2M - rmax = 1000.0 - - #Make Swifter files - plfile = open(swifter_pl, 'w') - print(f'{npl+1} ! Planet input file generated using init_cond.py using JPL Horizons data for the major planets (and Pluto) for epoch {tstart}' ,file=plfile) - print(f'1 {GMSun}',file=plfile) - print(f'0.0 0.0 0.0',file=plfile) - print(f'0.0 0.0 0.0',file=plfile) - for i, plid in enumerate(plvec): - print(f'{i + 2} {GMSun / MSun_over_Mpl[plid]} {Rhill[plid]}', file=plfile) - print(f'{Rpl[plid] / DU2M}', file=plfile) - print(f'{plvec[plid][0]} {plvec[plid][1]} {plvec[plid][2]}', file=plfile) - print(f'{plvec[plid][3]} {plvec[plid][4]} {plvec[plid][5]}', file=plfile) - plfile.close() - - tpfile = open(swifter_tp, 'w') - print(0,file=tpfile) - tpfile.close() - - sys.stdout = open(swifter_input, "w") - print(f'! Swifter input file generated using init_cond.py') - print(f'T0 {t_0} ') - print(f'TSTOP {end_sim}') - print(f'DT {deltaT}') - print(f'PL_IN {swifter_pl}') - print(f'TP_IN {swifter_tp}') - print(f'IN_TYPE ASCII') - print(f'ISTEP_OUT {iout:d}') - print(f'ISTEP_DUMP {iout:d}') - print(f'BIN_OUT {swifter_bin}') - print(f'OUT_TYPE REAL8') - print(f'OUT_FORM EL') - print(f'OUT_STAT NEW') - print(f'J2 {J2}') - print(f'J4 {J4}') - print(f'CHK_CLOSE yes') - print(f'CHK_RMIN {rmin}') - print(f'CHK_RMAX {rmax}') - print(f'CHK_EJECT {rmax}') - print(f'CHK_QMIN {rmin}') - print(f'CHK_QMIN_COORD HELIO') - print(f'CHK_QMIN_RANGE {rmin} {rmax}') - print(f'ENC_OUT {swifter_enc}') - print(f'EXTRA_FORCE no') - print(f'BIG_DISCARD no') - print(f'RHILL_PRESENT yes') - print(f'C {c / (DU2M / TU2S)}') - - #Now make Swiftest files - cbfile = open(swiftest_cb, 'w') - print(f'{1.0}',file=cbfile) - print(f'{rmin}',file=cbfile) - print(f'{J2}',file=cbfile) - print(f'{J4}',file=cbfile) - - plfile = open(swiftest_pl, 'w') - print(npl,file=plfile) - - for i, plid in enumerate(plvec): - print(f'{i + 2} {1.0 / MSun_over_Mpl[plid]}', file=plfile) - print(f'{Rpl[plid] / DU2M}', file=plfile) - print(f'{plvec[plid][0]} {plvec[plid][1]} {plvec[plid][2]}', file=plfile) - print(f'{plvec[plid][3]} {plvec[plid][4]} {plvec[plid][5]}', file=plfile) - plfile.close() - tpfile = open(swiftest_tp, 'w') - print(0,file=tpfile) - tpfile.close() - - sys.stdout = open(swiftest_input, "w") - print(f'! Swiftest input file generated using init_cond.py') - print(f'T0 {t_0} ') - print(f'TSTOP {end_sim}') - print(f'DT {deltaT}') - print(f'CB_IN {swiftest_cb}') - print(f'PL_IN {swiftest_pl}') - print(f'TP_IN {swiftest_tp}') - print(f'IN_TYPE ASCII') - print(f'ISTEP_OUT {iout:d}') - print(f'ISTEP_DUMP {iout:d}') - print(f'BIN_OUT {swiftest_bin}') - print(f'OUT_TYPE REAL8') - print(f'OUT_FORM EL') - print(f'OUT_STAT REPLACE') - print(f'CHK_CLOSE yes') - print(f'CHK_RMIN {rmin}') - print(f'CHK_RMAX {rmax}') - print(f'CHK_EJECT {rmax}') - print(f'CHK_QMIN {rmin}') - print(f'CHK_QMIN_COORD HELIO') - print(f'CHK_QMIN_RANGE {rmin} {rmax}') - print(f'ENC_OUT {swiftest_enc}') - print(f'EXTRA_FORCE no') - print(f'BIG_DISCARD no') - print(f'ROTATION no') - print(f'GR yes') - print(f'MU2KG {MU2KG}') - print(f'DU2M {DU2M}') - print(f'TU2S {TU2S}') - - - sys.stdout = sys.__stdout__ diff --git a/examples/whm_gr_test/swiftest_relativity.py b/examples/whm_gr_test/swiftest_relativity.py new file mode 100644 index 000000000..a4ea53c3b --- /dev/null +++ b/examples/whm_gr_test/swiftest_relativity.py @@ -0,0 +1,60 @@ +#!/usr/bin/env python +import swiftest +from astroquery.jplhorizons import Horizons +import datetime +import numpy as np +import matplotlib.pyplot as plt + +sim_gr = swiftest.Simulation(param_file="param.gr.in", output_file_name="bin.gr.nc") +sim_gr.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"]) + +sim_nogr = swiftest.Simulation(param_file="param.nogr.in", output_file_name="bin.nogr.nc") +sim_nogr.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"]) + +tstep_out = 10.0 +sim_gr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator="whm",general_relativity=True) +sim_nogr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator="whm",general_relativity=False) + +# Get the start and end date of the simulation so we can compare with the real solar system +start_date = sim_gr.ephemeris_date +tstop_d = sim_gr.param['TSTOP'] * sim_gr.param['TU2S'] / swiftest.JD2S + +stop_date = (datetime.datetime.fromisoformat(start_date) + datetime.timedelta(days=tstop_d)).isoformat() + +#Get the ephemerides of Mercury for the same timeframe as the simulation +obj = Horizons(id='1', location='@sun', + epochs={'start':start_date, 'stop':stop_date, + 'step':'10y'}) +el = obj.elements() +t = (el['datetime_jd']-el['datetime_jd'][0]) / 365.25 +varpi_obs = el['w'] + el['Omega'] + +# Compute the longitude of the periapsis +sim_gr.data['varpi'] = np.mod(sim_gr.data['omega'] + sim_gr.data['capom'],360) +sim_nogr.data['varpi'] = np.mod(sim_nogr.data['omega'] + sim_nogr.data['capom'],360) + +varpisim_gr= sim_gr.data['varpi'].sel(name="Mercury") +varpisim_nogr= sim_nogr.data['varpi'].sel(name="Mercury") +tsim = sim_gr.data['time'] + +dvarpi_gr = np.diff(varpisim_gr) * 3600 * 100 / tstep_out +dvarpi_nogr = np.diff(varpisim_nogr) * 3600 * 100 / tstep_out +dvarpi_obs = np.diff(varpi_obs) / np.diff(t) * 3600 * 100 + + +fig, ax = plt.subplots() + +ax.plot(t, varpi_obs, label="JPL Horizons",linewidth=2.5) +ax.plot(tsim, varpisim_gr, label="Swiftest WHM GR",linewidth=1.5) +ax.plot(tsim, varpisim_nogr, label="Swiftest WHM No GR",linewidth=1.5) +ax.set_xlabel('Time (y)') +ax.set_ylabel('Mercury $\\varpi$ (deg)') +ax.legend() +plt.savefig("whm_gr_mercury_precession.png",dpi=300) + +print('Mean precession rate for Mercury long. peri. (arcsec/100 y)') +print(f'JPL Horizons : {np.mean(dvarpi_obs)}') +print(f'Swiftest No GR : {np.mean(dvarpi_nogr)}') +print(f'Swiftest GR : {np.mean(dvarpi_gr)}') +print(f'Obs - Swiftest GR : {np.mean(dvarpi_obs - dvarpi_gr)}') +print(f'Obs - Swiftest No GR : {np.mean(dvarpi_obs - dvarpi_nogr)}') From 848cecc10babc8b0f8d40c42a4fab4c2a551aa87 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 23:01:19 -0500 Subject: [PATCH 12/13] Added GR test script and Notebook to helio_gr_test --- .../helio_gr_test/swiftest_relativity.ipynb | 177 ++++++++++++++++++ examples/helio_gr_test/swiftest_relativity.py | 60 ++++++ 2 files changed, 237 insertions(+) create mode 100644 examples/helio_gr_test/swiftest_relativity.ipynb create mode 100644 examples/helio_gr_test/swiftest_relativity.py diff --git a/examples/helio_gr_test/swiftest_relativity.ipynb b/examples/helio_gr_test/swiftest_relativity.ipynb new file mode 100644 index 000000000..6abb9524a --- /dev/null +++ b/examples/helio_gr_test/swiftest_relativity.ipynb @@ -0,0 +1,177 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import swiftest\n", + "from astroquery.jplhorizons import Horizons\n", + "import datetime\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sim_gr = swiftest.Simulation(param_file=\"param.gr.in\", output_file_name=\"bin.gr.nc\")\n", + "sim_gr.add_solar_system_body([\"Sun\",\"Mercury\",\"Venus\",\"Earth\",\"Mars\",\"Jupiter\",\"Saturn\",\"Uranus\",\"Neptune\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sim_nogr = swiftest.Simulation(param_file=\"param.nogr.in\", output_file_name=\"bin.nogr.nc\")\n", + "sim_nogr.add_solar_system_body([\"Sun\",\"Mercury\",\"Venus\",\"Earth\",\"Mars\",\"Jupiter\",\"Saturn\",\"Uranus\",\"Neptune\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "tstep_out = 10.0\n", + "sim_gr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator=\"whm\",general_relativity=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "sim_nogr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator=\"whm\",general_relativity=False)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the start and end date of the simulation so we can compare with the real solar system\n", + "start_date = sim_gr.ephemeris_date\n", + "tstop_d = sim_gr.param['TSTOP'] * sim_gr.param['TU2S'] / swiftest.JD2S\n", + "\n", + "stop_date = (datetime.datetime.fromisoformat(start_date) + datetime.timedelta(days=tstop_d)).isoformat()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Get the ephemerides of Mercury for the same timeframe as the simulation\n", + "obj = Horizons(id='1', location='@sun',\n", + " epochs={'start':start_date, 'stop':stop_date,\n", + " 'step':'10y'})\n", + "el = obj.elements()\n", + "t = (el['datetime_jd']-el['datetime_jd'][0]) / 365.25\n", + "varpi_obs = el['w'] + el['Omega']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Compute the longitude of the periapsis\n", + "sim_gr.data['varpi'] = np.mod(sim_gr.data['omega'] + sim_gr.data['capom'],360)\n", + "sim_nogr.data['varpi'] = np.mod(sim_nogr.data['omega'] + sim_nogr.data['capom'],360)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "varpisim_gr= sim_gr.data['varpi'].sel(name=\"Mercury\")\n", + "varpisim_nogr= sim_nogr.data['varpi'].sel(name=\"Mercury\")\n", + "tsim = sim_gr.data['time']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dvarpi_gr = np.diff(varpisim_gr) * 3600 * 100 / tstep_out\n", + "dvarpi_nogr = np.diff(varpisim_nogr) * 3600 * 100 / tstep_out\n", + "dvarpi_obs = np.diff(varpi_obs) / np.diff(t) * 3600 * 100" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "\n", + "ax.plot(t, varpi_obs, label=\"JPL Horizons\",linewidth=2.5)\n", + "ax.plot(tsim, varpisim_gr, label=\"Swiftest WHM GR\",linewidth=1.5)\n", + "ax.plot(tsim, varpisim_nogr, label=\"Swiftest WHM No GR\",linewidth=1.5)\n", + "ax.set_xlabel('Time (y)')\n", + "ax.set_ylabel('Mercury $\\\\varpi$ (deg)')\n", + "ax.legend()\n", + "plt.savefig(\"whm_gr_mercury_precession.png\",dpi=300)\n", + "print('Mean precession rate for Mercury long. peri. (arcsec/100 y)')\n", + "print(f'JPL Horizons : {np.mean(dvarpi_obs)}')\n", + "print(f'Swiftest No GR : {np.mean(dvarpi_nogr)}')\n", + "print(f'Swiftest GR : {np.mean(dvarpi_gr)}')\n", + "print(f'Obs - Swiftest GR : {np.mean(dvarpi_obs - dvarpi_gr)}')\n", + "print(f'Obs - Swiftest No GR : {np.mean(dvarpi_obs - dvarpi_nogr)}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "swiftest", + "language": "python", + "name": "swiftest" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/helio_gr_test/swiftest_relativity.py b/examples/helio_gr_test/swiftest_relativity.py new file mode 100644 index 000000000..a5f4e4371 --- /dev/null +++ b/examples/helio_gr_test/swiftest_relativity.py @@ -0,0 +1,60 @@ +#!/usr/bin/env python +import swiftest +from astroquery.jplhorizons import Horizons +import datetime +import numpy as np +import matplotlib.pyplot as plt + +sim_gr = swiftest.Simulation(param_file="param.gr.in", output_file_name="bin.gr.nc") +sim_gr.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"]) + +sim_nogr = swiftest.Simulation(param_file="param.nogr.in", output_file_name="bin.nogr.nc") +sim_nogr.add_solar_system_body(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"]) + +tstep_out = 10.0 +sim_gr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator="helio",general_relativity=True) +sim_nogr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator="helio",general_relativity=False) + +# Get the start and end date of the simulation so we can compare with the real solar system +start_date = sim_gr.ephemeris_date +tstop_d = sim_gr.param['TSTOP'] * sim_gr.param['TU2S'] / swiftest.JD2S + +stop_date = (datetime.datetime.fromisoformat(start_date) + datetime.timedelta(days=tstop_d)).isoformat() + +#Get the ephemerides of Mercury for the same timeframe as the simulation +obj = Horizons(id='1', location='@sun', + epochs={'start':start_date, 'stop':stop_date, + 'step':'10y'}) +el = obj.elements() +t = (el['datetime_jd']-el['datetime_jd'][0]) / 365.25 +varpi_obs = el['w'] + el['Omega'] + +# Compute the longitude of the periapsis +sim_gr.data['varpi'] = np.mod(sim_gr.data['omega'] + sim_gr.data['capom'],360) +sim_nogr.data['varpi'] = np.mod(sim_nogr.data['omega'] + sim_nogr.data['capom'],360) + +varpisim_gr= sim_gr.data['varpi'].sel(name="Mercury") +varpisim_nogr= sim_nogr.data['varpi'].sel(name="Mercury") +tsim = sim_gr.data['time'] + +dvarpi_gr = np.diff(varpisim_gr) * 3600 * 100 / tstep_out +dvarpi_nogr = np.diff(varpisim_nogr) * 3600 * 100 / tstep_out +dvarpi_obs = np.diff(varpi_obs) / np.diff(t) * 3600 * 100 + + +fig, ax = plt.subplots() + +ax.plot(t, varpi_obs, label="JPL Horizons",linewidth=2.5) +ax.plot(tsim, varpisim_gr, label="Swiftest Helio GR",linewidth=1.5) +ax.plot(tsim, varpisim_nogr, label="Swiftest Helio No GR",linewidth=1.5) +ax.set_xlabel('Time (y)') +ax.set_ylabel('Mercury $\\varpi$ (deg)') +ax.legend() +plt.savefig("helio_gr_mercury_precession.png",dpi=300) + +print('Mean precession rate for Mercury long. peri. (arcsec/100 y)') +print(f'JPL Horizons : {np.mean(dvarpi_obs)}') +print(f'Swiftest No GR : {np.mean(dvarpi_nogr)}') +print(f'Swiftest GR : {np.mean(dvarpi_gr)}') +print(f'Obs - Swiftest GR : {np.mean(dvarpi_obs - dvarpi_gr)}') +print(f'Obs - Swiftest No GR : {np.mean(dvarpi_obs - dvarpi_nogr)}') From 60cf506fce744ace346e3f1b33a9fd93121851d7 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 17 Nov 2022 23:05:59 -0500 Subject: [PATCH 13/13] Fixed typos in relativity test scripts --- .../helio_gr_test/swiftest_relativity.ipynb | 10 +- .../whm_gr_test/swiftest_relativity.ipynb | 894 +++++++++++++++++- 2 files changed, 891 insertions(+), 13 deletions(-) diff --git a/examples/helio_gr_test/swiftest_relativity.ipynb b/examples/helio_gr_test/swiftest_relativity.ipynb index 6abb9524a..6946ef658 100644 --- a/examples/helio_gr_test/swiftest_relativity.ipynb +++ b/examples/helio_gr_test/swiftest_relativity.ipynb @@ -41,7 +41,7 @@ "source": [ "%%capture\n", "tstep_out = 10.0\n", - "sim_gr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator=\"whm\",general_relativity=True)" + "sim_gr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator=\"helio\",general_relativity=True)" ] }, { @@ -51,7 +51,7 @@ "outputs": [], "source": [ "%%capture\n", - "sim_nogr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator=\"whm\",general_relativity=False)" + "sim_nogr.run(tstop=1000.0, dt=0.005, tstep_out=tstep_out, integrator=\"helio\",general_relativity=False)" ] }, { @@ -124,12 +124,12 @@ "fig, ax = plt.subplots()\n", "\n", "ax.plot(t, varpi_obs, label=\"JPL Horizons\",linewidth=2.5)\n", - "ax.plot(tsim, varpisim_gr, label=\"Swiftest WHM GR\",linewidth=1.5)\n", - "ax.plot(tsim, varpisim_nogr, label=\"Swiftest WHM No GR\",linewidth=1.5)\n", + "ax.plot(tsim, varpisim_gr, label=\"Swiftest helio GR\",linewidth=1.5)\n", + "ax.plot(tsim, varpisim_nogr, label=\"Swiftest helio No GR\",linewidth=1.5)\n", "ax.set_xlabel('Time (y)')\n", "ax.set_ylabel('Mercury $\\\\varpi$ (deg)')\n", "ax.legend()\n", - "plt.savefig(\"whm_gr_mercury_precession.png\",dpi=300)\n", + "plt.savefig(\"helio_gr_mercury_precession.png\",dpi=300)\n", "print('Mean precession rate for Mercury long. peri. (arcsec/100 y)')\n", "print(f'JPL Horizons : {np.mean(dvarpi_obs)}')\n", "print(f'Swiftest No GR : {np.mean(dvarpi_nogr)}')\n", diff --git a/examples/whm_gr_test/swiftest_relativity.ipynb b/examples/whm_gr_test/swiftest_relativity.ipynb index 6abb9524a..113e10f81 100644 --- a/examples/whm_gr_test/swiftest_relativity.ipynb +++ b/examples/whm_gr_test/swiftest_relativity.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -15,9 +15,448 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Creating the Sun as a central body\n", + "Fetching ephemerides data for Mercury from JPL/Horizons\n", + "Fetching ephemerides data for Venus from JPL/Horizons\n", + "Fetching ephemerides data for Earth from JPL/Horizons\n", + "Fetching ephemerides data for Mars from JPL/Horizons\n", + "Fetching ephemerides data for Jupiter from JPL/Horizons\n", + "Fetching ephemerides data for Saturn from JPL/Horizons\n", + "Fetching ephemerides data for Uranus from JPL/Horizons\n", + "Fetching ephemerides data for Neptune from JPL/Horizons\n", + "Writing initial conditions to file init_cond.nc\n", + "Writing parameter inputs to file /home/daminton/git_debug/swiftest/examples/whm_gr_test/param.gr.in\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:        (name: 9, time: 1)\n",
+       "Coordinates:\n",
+       "  * name           (name) <U32 'Sun' 'Mercury' 'Venus' ... 'Uranus' 'Neptune'\n",
+       "  * time           (time) float64 0.0\n",
+       "Data variables: (12/14)\n",
+       "    particle_type  (name) <U32 'Central Body' 'Massive Body' ... 'Massive Body'\n",
+       "    id             (name) int64 0 1 2 3 4 5 6 7 8\n",
+       "    a              (time, name) float64 nan 0.3871 0.7233 ... 9.532 19.24 30.04\n",
+       "    e              (time, name) float64 nan 0.2056 0.006718 ... 0.04796 0.008956\n",
+       "    inc            (time, name) float64 nan 7.003 3.394 ... 2.488 0.773 1.771\n",
+       "    capom          (time, name) float64 nan 48.3 76.6 ... 113.6 74.01 131.8\n",
+       "    ...             ...\n",
+       "    Gmass          (time, name) float64 39.48 6.554e-06 ... 0.001724 0.002034\n",
+       "    radius         (time, name) float64 0.00465 1.631e-05 ... 0.0001646\n",
+       "    j2rp2          (time, name) float64 4.754e-12 nan nan nan ... nan nan nan\n",
+       "    j4rp4          (time, name) float64 -2.247e-18 nan nan nan ... nan nan nan\n",
+       "    ntp            (time) int64 0\n",
+       "    npl            (time) int64 8
" + ], + "text/plain": [ + "\n", + "Dimensions: (name: 9, time: 1)\n", + "Coordinates:\n", + " * name (name) \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:        (name: 9, time: 1)\n",
+       "Coordinates:\n",
+       "  * name           (name) <U32 'Sun' 'Mercury' 'Venus' ... 'Uranus' 'Neptune'\n",
+       "  * time           (time) float64 0.0\n",
+       "Data variables: (12/14)\n",
+       "    particle_type  (name) <U32 'Central Body' 'Massive Body' ... 'Massive Body'\n",
+       "    id             (name) int64 0 1 2 3 4 5 6 7 8\n",
+       "    a              (time, name) float64 nan 0.3871 0.7233 ... 9.532 19.24 30.04\n",
+       "    e              (time, name) float64 nan 0.2056 0.006718 ... 0.04796 0.008956\n",
+       "    inc            (time, name) float64 nan 7.003 3.394 ... 2.488 0.773 1.771\n",
+       "    capom          (time, name) float64 nan 48.3 76.6 ... 113.6 74.01 131.8\n",
+       "    ...             ...\n",
+       "    Gmass          (time, name) float64 39.48 6.554e-06 ... 0.001724 0.002034\n",
+       "    radius         (time, name) float64 0.00465 1.631e-05 ... 0.0001646\n",
+       "    j2rp2          (time, name) float64 4.754e-12 nan nan nan ... nan nan nan\n",
+       "    j4rp4          (time, name) float64 -2.247e-18 nan nan nan ... nan nan nan\n",
+       "    ntp            (time) int64 0\n",
+       "    npl            (time) int64 8
" + ], + "text/plain": [ + "\n", + "Dimensions: (name: 9, time: 1)\n", + "Coordinates:\n", + " * name (name)