From c2ddf3e2988deda56aa0965396ea71f0c1963f59 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 12 Nov 2021 15:35:52 -0500 Subject: [PATCH 01/10] Tweaks to compiler flags --- Makefile.Defines | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/Makefile.Defines b/Makefile.Defines index a5417a18d..046d9d3f9 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -41,9 +41,6 @@ SWIFTEST_HOME = $(ROOT_DIR) USER_MODULES = COLLRESOLVE_HOME = $(ROOT_DIR)/collresolve/ -#NETCDF_FORTRAN_HOME = /home/daminton/local - - # Compiler definitions # DO NOT include in FFLAGS the "-c" option to compile object only @@ -69,8 +66,8 @@ GPRODUCTION = -O3 -ffree-line-length-none $(GPAR) INCLUDES = -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include -I$(MKLROOT)/include LINKS = -L$(MKLROOT)/lib/intel64 -L$(NETCDF_FORTRAN_HOME)/lib -L$(ADVISOR_2019_DIR)/lib64 -lswiftest -lnetcdf -lnetcdff -qopt-matmul $(PAR) -FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) -g -traceback -FFLAGS = $(IPRODUCTION) -fp-model=fast -g -traceback +FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) +FFLAGS = $(IPRODUCTION) -fp-model=fast FORTRAN = ifort AR = xiar From 125946f35741fb0655ec112f56a37c46b91058af Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 15 Nov 2021 15:26:58 -0500 Subject: [PATCH 02/10] Added more diagnostic info to the mass conservation failure --- src/io/io.f90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/io/io.f90 b/src/io/io.f90 index ebdd68531..c657ec058 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -72,6 +72,10 @@ module subroutine io_conservation_report(self, param, lterminal) if (lterminal) write(*, EGYTERMFMT) Lerror, Ecoll_error, Etotal_error, Merror if (abs(Merror) > 100 * epsilon(Merror)) then write(*,*) "Severe error! Mass not conserved! Halting!" + write(*,*) "Merror = ", Merror + write(*,*) "GMtot_now : ",GMtot_now + write(*,*) "GMtot_orig: ",system%GMtot_orig + write(*,*) "Difference: ",GMtot_now - system%GMtot_orig call pl%xv2el(cb) call self%write_hdr(param%nciu, param) call cb%write_frame(param%nciu, param) From 733e68ebad97b816844834900b237fc402698c22 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 24 Nov 2021 11:09:48 -0500 Subject: [PATCH 03/10] Updated NetCDF reads and potential energy loop --- src/netcdf/netcdf.f90 | 87 ++++++++++++++------------- src/util/util_get_energy_momentum.f90 | 13 +++- 2 files changed, 56 insertions(+), 44 deletions(-) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index d28fcc788..6c1e05e7a 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -250,7 +250,8 @@ module subroutine netcdf_initialize_output(self, param) call check( nf90_def_var(self%ncid, ORIGIN_TIME_VARNAME, self%out_type, self%id_dimid, self%origin_time_varid) ) !call check( nf90_def_var_chunking(self%ncid, self%origin_time_varid, NF90_CHUNKED, [self%id_chunk]) ) - call check( nf90_def_var(self%ncid, ORIGIN_TYPE_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], self%origin_type_varid) ) + call check( nf90_def_var(self%ncid, ORIGIN_TYPE_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], & + self%origin_type_varid) ) !call check( nf90_def_var_chunking(self%ncid, self%origin_type_varid, NF90_CHUNKED, [NAMELEN, self%id_chunk]) ) call check( nf90_def_var(self%ncid, ORIGIN_XHX_VARNAME, self%out_type, self%id_dimid, self%origin_xhx_varid) ) !call check( nf90_def_var_chunking(self%ncid, self%origin_xhx_varid, NF90_CHUNKED, [self%id_chunk]) ) @@ -543,55 +544,55 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr) ! Now read in each variable and split the outputs by body type if ((param%in_form == XV) .or. (param%in_form == XVEL)) then call check( nf90_get_var(iu%ncid, iu%xhx_varid, rtemp, start=[1, tslot]) ) - pl%xh(1,:) = pack(rtemp, plmask) - tp%xh(1,:) = pack(rtemp, tpmask) + if (npl > 0) pl%xh(1,:) = pack(rtemp, plmask) + if (ntp > 0) tp%xh(1,:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%xhy_varid, rtemp, start=[1, tslot]) ) - pl%xh(2,:) = pack(rtemp, plmask) - tp%xh(2,:) = pack(rtemp, tpmask) + if (npl > 0) pl%xh(2,:) = pack(rtemp, plmask) + if (ntp > 0) tp%xh(2,:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%xhz_varid, rtemp, start=[1, tslot]) ) - pl%xh(3,:) = pack(rtemp, plmask) - tp%xh(3,:) = pack(rtemp, tpmask) + if (npl > 0) pl%xh(3,:) = pack(rtemp, plmask) + if (ntp > 0) tp%xh(3,:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%vhx_varid, rtemp, start=[1, tslot]) ) - pl%vh(1,:) = pack(rtemp, plmask) - tp%vh(1,:) = pack(rtemp, tpmask) + if (npl > 0) pl%vh(1,:) = pack(rtemp, plmask) + if (ntp > 0) tp%vh(1,:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%vhy_varid, rtemp, start=[1, tslot]) ) - pl%vh(2,:) = pack(rtemp, plmask) - tp%vh(2,:) = pack(rtemp, tpmask) + if (npl > 0) pl%vh(2,:) = pack(rtemp, plmask) + if (ntp > 0) tp%vh(2,:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%vhz_varid, rtemp, start=[1, tslot]) ) - pl%vh(3,:) = pack(rtemp, plmask) - tp%vh(3,:) = pack(rtemp, tpmask) + if (npl > 0) pl%vh(3,:) = pack(rtemp, plmask) + if (ntp > 0) tp%vh(3,:) = pack(rtemp, tpmask) end if if ((param%in_form == EL) .or. (param%in_form == XVEL)) then call check( nf90_get_var(iu%ncid, iu%a_varid, rtemp, start=[1, tslot]) ) - pl%a(:) = pack(rtemp, plmask) - tp%a(:) = pack(rtemp, tpmask) + if (npl > 0) pl%a(:) = pack(rtemp, plmask) + if (ntp > 0) tp%a(:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%e_varid, rtemp, start=[1, tslot]) ) - pl%e(:) = pack(rtemp, plmask) - tp%e(:) = pack(rtemp, tpmask) + if (npl > 0) pl%e(:) = pack(rtemp, plmask) + if (ntp > 0) tp%e(:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%inc_varid, rtemp, start=[1, tslot]) ) - pl%inc(:) = pack(rtemp, plmask) - tp%inc(:) = pack(rtemp, tpmask) + 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]) ) - pl%capom(:) = pack(rtemp, plmask) - tp%capom(:) = pack(rtemp, tpmask) + 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]) ) - pl%omega(:) = pack(rtemp, plmask) - tp%omega(:) = pack(rtemp, tpmask) + 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]) ) - pl%capm(:) = pack(rtemp, plmask) - tp%capm(:) = pack(rtemp, tpmask) + if (npl > 0) pl%capm(:) = pack(rtemp, plmask) + if (ntp > 0) tp%capm(:) = pack(rtemp, tpmask) end if @@ -599,57 +600,59 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr) cb%Gmass = rtemp(1) cb%mass = cb%Gmass / param%GU - pl%Gmass(:) = pack(rtemp, plmask) - pl%mass(:) = pl%Gmass(:) / param%GU + if (npl > 0) then + pl%Gmass(:) = pack(rtemp, plmask) + pl%mass(:) = pl%Gmass(:) / param%GU - if (param%lrhill_present) then - call check( nf90_get_var(iu%ncid, iu%rhill_varid, rtemp, start=[1, tslot]) ) - pl%rhill(:) = pack(rtemp, plmask) + if (param%lrhill_present) then + call check( nf90_get_var(iu%ncid, iu%rhill_varid, rtemp, start=[1, tslot]) ) + pl%rhill(:) = pack(rtemp, plmask) + end if end if if (param%lclose) then call check( nf90_get_var(iu%ncid, iu%radius_varid, rtemp, start=[1, tslot]) ) cb%radius = rtemp(1) - pl%radius(:) = pack(rtemp, plmask) + if (npl > 0) pl%radius(:) = pack(rtemp, plmask) else cb%radius = param%rmin - pl%radius(:) = 0.0_DP + if (npl > 0) pl%radius(:) = 0.0_DP end if if (param%lrotation) then call check( nf90_get_var(iu%ncid, iu%Ip1_varid, rtemp, start=[1, tslot]) ) cb%Ip(1) = rtemp(1) - pl%Ip(1,:) = pack(rtemp, plmask) + if (npl > 0) pl%Ip(1,:) = pack(rtemp, plmask) call check( nf90_get_var(iu%ncid, iu%Ip2_varid, rtemp, start=[1, tslot]) ) cb%Ip(2) = rtemp(1) - pl%Ip(2,:) = pack(rtemp, plmask) + if (npl > 0) pl%Ip(2,:) = pack(rtemp, plmask) call check( nf90_get_var(iu%ncid, iu%Ip3_varid, rtemp, start=[1, tslot]) ) cb%Ip(3) = rtemp(1) - pl%Ip(3,:) = pack(rtemp, plmask) + if (npl > 0) pl%Ip(3,:) = pack(rtemp, plmask) call check( nf90_get_var(iu%ncid, iu%rotx_varid, rtemp, start=[1, tslot]) ) cb%rot(1) = rtemp(1) - pl%rot(1,:) = pack(rtemp, plmask) + if (npl > 0) pl%rot(1,:) = pack(rtemp, plmask) call check( nf90_get_var(iu%ncid, iu%roty_varid, rtemp, start=[1, tslot]) ) cb%rot(2) = rtemp(1) - pl%rot(2,:) = pack(rtemp, plmask) + if (npl > 0) pl%rot(2,:) = pack(rtemp, plmask) call check( nf90_get_var(iu%ncid, iu%rotz_varid, rtemp, start=[1, tslot]) ) cb%rot(3) = rtemp(1) - pl%rot(3,:) = pack(rtemp, plmask) + if (npl > 0) pl%rot(3,:) = pack(rtemp, plmask) end if if (param%ltides) then call check( nf90_get_var(iu%ncid, iu%k2_varid, rtemp, start=[1, tslot]) ) cb%k2 = rtemp(1) - pl%k2(:) = pack(rtemp, plmask) + if (npl > 0) pl%k2(:) = pack(rtemp, plmask) call check( nf90_get_var(iu%ncid, iu%Q_varid, rtemp, start=[1, tslot]) ) cb%Q = rtemp(1) - pl%Q(:) = pack(rtemp, plmask) + if (npl > 0) pl%Q(:) = pack(rtemp, plmask) end if call self%read_particle_info(iu, param, plmask, tpmask) @@ -936,7 +939,9 @@ module subroutine netcdf_write_frame_base(self, iu, param) select type(self) class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body call check( nf90_put_var(iu%ncid, iu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]) ) - if (param%lrhill_present) call check( nf90_put_var(iu%ncid, iu%rhill_varid, self%rhill(j), start=[idslot, tslot]) ) + if (param%lrhill_present) then + call check( nf90_put_var(iu%ncid, iu%rhill_varid, self%rhill(j), start=[idslot, tslot]) ) + end if if (param%lclose) call check( nf90_put_var(iu%ncid, iu%radius_varid, self%radius(j), start=[idslot, tslot]) ) if (param%lrotation) then call check( nf90_put_var(iu%ncid, iu%Ip1_varid, self%Ip(1, j), start=[idslot, tslot]) ) diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index 3474e2921..a86f5efd4 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -189,10 +189,17 @@ subroutine util_get_energy_potential_triangular(npl, lmask, GMcb, Gmass, mass, x pecb(i) = -GMcb * mass(i) / norm2(xb(:,i)) end do - pe = sum(pecb(1:npl), lmask(1:npl)) - do concurrent(i = 1:npl, j = 1:npl, (j > i) .and. lmask(i) .and. lmask(j)) - pe = pe - (Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j)) + pe = 0.0_DP + !$omp parallel do default(private) schedule(static)& + !$omp shared(npl, lmask, Gmass, mass, xb) & + !$omp reduction(-:pe) + do i = 1, npl + do concurrent(j = i+1:npl, lmask(i) .and. lmask(j)) + pe = pe - (Gmass(i) * mass(j)) / norm2(xb(:, i) - xb(:, j)) + end do end do + !$omp end parallel do + pe = pe + sum(pecb(1:npl), lmask(1:npl)) return end subroutine util_get_energy_potential_triangular From 8eac6a0f3d1f9b904f357ac48ba721e707aadd89 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 24 Nov 2021 11:14:46 -0500 Subject: [PATCH 04/10] Updated NetCDF reads from the profiling branch --- src/netcdf/netcdf.f90 | 87 +++++++++++++++++++++++-------------------- 1 file changed, 46 insertions(+), 41 deletions(-) diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index d28fcc788..6c1e05e7a 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -250,7 +250,8 @@ module subroutine netcdf_initialize_output(self, param) call check( nf90_def_var(self%ncid, ORIGIN_TIME_VARNAME, self%out_type, self%id_dimid, self%origin_time_varid) ) !call check( nf90_def_var_chunking(self%ncid, self%origin_time_varid, NF90_CHUNKED, [self%id_chunk]) ) - call check( nf90_def_var(self%ncid, ORIGIN_TYPE_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], self%origin_type_varid) ) + call check( nf90_def_var(self%ncid, ORIGIN_TYPE_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], & + self%origin_type_varid) ) !call check( nf90_def_var_chunking(self%ncid, self%origin_type_varid, NF90_CHUNKED, [NAMELEN, self%id_chunk]) ) call check( nf90_def_var(self%ncid, ORIGIN_XHX_VARNAME, self%out_type, self%id_dimid, self%origin_xhx_varid) ) !call check( nf90_def_var_chunking(self%ncid, self%origin_xhx_varid, NF90_CHUNKED, [self%id_chunk]) ) @@ -543,55 +544,55 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr) ! Now read in each variable and split the outputs by body type if ((param%in_form == XV) .or. (param%in_form == XVEL)) then call check( nf90_get_var(iu%ncid, iu%xhx_varid, rtemp, start=[1, tslot]) ) - pl%xh(1,:) = pack(rtemp, plmask) - tp%xh(1,:) = pack(rtemp, tpmask) + if (npl > 0) pl%xh(1,:) = pack(rtemp, plmask) + if (ntp > 0) tp%xh(1,:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%xhy_varid, rtemp, start=[1, tslot]) ) - pl%xh(2,:) = pack(rtemp, plmask) - tp%xh(2,:) = pack(rtemp, tpmask) + if (npl > 0) pl%xh(2,:) = pack(rtemp, plmask) + if (ntp > 0) tp%xh(2,:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%xhz_varid, rtemp, start=[1, tslot]) ) - pl%xh(3,:) = pack(rtemp, plmask) - tp%xh(3,:) = pack(rtemp, tpmask) + if (npl > 0) pl%xh(3,:) = pack(rtemp, plmask) + if (ntp > 0) tp%xh(3,:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%vhx_varid, rtemp, start=[1, tslot]) ) - pl%vh(1,:) = pack(rtemp, plmask) - tp%vh(1,:) = pack(rtemp, tpmask) + if (npl > 0) pl%vh(1,:) = pack(rtemp, plmask) + if (ntp > 0) tp%vh(1,:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%vhy_varid, rtemp, start=[1, tslot]) ) - pl%vh(2,:) = pack(rtemp, plmask) - tp%vh(2,:) = pack(rtemp, tpmask) + if (npl > 0) pl%vh(2,:) = pack(rtemp, plmask) + if (ntp > 0) tp%vh(2,:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%vhz_varid, rtemp, start=[1, tslot]) ) - pl%vh(3,:) = pack(rtemp, plmask) - tp%vh(3,:) = pack(rtemp, tpmask) + if (npl > 0) pl%vh(3,:) = pack(rtemp, plmask) + if (ntp > 0) tp%vh(3,:) = pack(rtemp, tpmask) end if if ((param%in_form == EL) .or. (param%in_form == XVEL)) then call check( nf90_get_var(iu%ncid, iu%a_varid, rtemp, start=[1, tslot]) ) - pl%a(:) = pack(rtemp, plmask) - tp%a(:) = pack(rtemp, tpmask) + if (npl > 0) pl%a(:) = pack(rtemp, plmask) + if (ntp > 0) tp%a(:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%e_varid, rtemp, start=[1, tslot]) ) - pl%e(:) = pack(rtemp, plmask) - tp%e(:) = pack(rtemp, tpmask) + if (npl > 0) pl%e(:) = pack(rtemp, plmask) + if (ntp > 0) tp%e(:) = pack(rtemp, tpmask) call check( nf90_get_var(iu%ncid, iu%inc_varid, rtemp, start=[1, tslot]) ) - pl%inc(:) = pack(rtemp, plmask) - tp%inc(:) = pack(rtemp, tpmask) + 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]) ) - pl%capom(:) = pack(rtemp, plmask) - tp%capom(:) = pack(rtemp, tpmask) + 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]) ) - pl%omega(:) = pack(rtemp, plmask) - tp%omega(:) = pack(rtemp, tpmask) + 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]) ) - pl%capm(:) = pack(rtemp, plmask) - tp%capm(:) = pack(rtemp, tpmask) + if (npl > 0) pl%capm(:) = pack(rtemp, plmask) + if (ntp > 0) tp%capm(:) = pack(rtemp, tpmask) end if @@ -599,57 +600,59 @@ module function netcdf_read_frame_system(self, iu, param) result(ierr) cb%Gmass = rtemp(1) cb%mass = cb%Gmass / param%GU - pl%Gmass(:) = pack(rtemp, plmask) - pl%mass(:) = pl%Gmass(:) / param%GU + if (npl > 0) then + pl%Gmass(:) = pack(rtemp, plmask) + pl%mass(:) = pl%Gmass(:) / param%GU - if (param%lrhill_present) then - call check( nf90_get_var(iu%ncid, iu%rhill_varid, rtemp, start=[1, tslot]) ) - pl%rhill(:) = pack(rtemp, plmask) + if (param%lrhill_present) then + call check( nf90_get_var(iu%ncid, iu%rhill_varid, rtemp, start=[1, tslot]) ) + pl%rhill(:) = pack(rtemp, plmask) + end if end if if (param%lclose) then call check( nf90_get_var(iu%ncid, iu%radius_varid, rtemp, start=[1, tslot]) ) cb%radius = rtemp(1) - pl%radius(:) = pack(rtemp, plmask) + if (npl > 0) pl%radius(:) = pack(rtemp, plmask) else cb%radius = param%rmin - pl%radius(:) = 0.0_DP + if (npl > 0) pl%radius(:) = 0.0_DP end if if (param%lrotation) then call check( nf90_get_var(iu%ncid, iu%Ip1_varid, rtemp, start=[1, tslot]) ) cb%Ip(1) = rtemp(1) - pl%Ip(1,:) = pack(rtemp, plmask) + if (npl > 0) pl%Ip(1,:) = pack(rtemp, plmask) call check( nf90_get_var(iu%ncid, iu%Ip2_varid, rtemp, start=[1, tslot]) ) cb%Ip(2) = rtemp(1) - pl%Ip(2,:) = pack(rtemp, plmask) + if (npl > 0) pl%Ip(2,:) = pack(rtemp, plmask) call check( nf90_get_var(iu%ncid, iu%Ip3_varid, rtemp, start=[1, tslot]) ) cb%Ip(3) = rtemp(1) - pl%Ip(3,:) = pack(rtemp, plmask) + if (npl > 0) pl%Ip(3,:) = pack(rtemp, plmask) call check( nf90_get_var(iu%ncid, iu%rotx_varid, rtemp, start=[1, tslot]) ) cb%rot(1) = rtemp(1) - pl%rot(1,:) = pack(rtemp, plmask) + if (npl > 0) pl%rot(1,:) = pack(rtemp, plmask) call check( nf90_get_var(iu%ncid, iu%roty_varid, rtemp, start=[1, tslot]) ) cb%rot(2) = rtemp(1) - pl%rot(2,:) = pack(rtemp, plmask) + if (npl > 0) pl%rot(2,:) = pack(rtemp, plmask) call check( nf90_get_var(iu%ncid, iu%rotz_varid, rtemp, start=[1, tslot]) ) cb%rot(3) = rtemp(1) - pl%rot(3,:) = pack(rtemp, plmask) + if (npl > 0) pl%rot(3,:) = pack(rtemp, plmask) end if if (param%ltides) then call check( nf90_get_var(iu%ncid, iu%k2_varid, rtemp, start=[1, tslot]) ) cb%k2 = rtemp(1) - pl%k2(:) = pack(rtemp, plmask) + if (npl > 0) pl%k2(:) = pack(rtemp, plmask) call check( nf90_get_var(iu%ncid, iu%Q_varid, rtemp, start=[1, tslot]) ) cb%Q = rtemp(1) - pl%Q(:) = pack(rtemp, plmask) + if (npl > 0) pl%Q(:) = pack(rtemp, plmask) end if call self%read_particle_info(iu, param, plmask, tpmask) @@ -936,7 +939,9 @@ module subroutine netcdf_write_frame_base(self, iu, param) select type(self) class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body call check( nf90_put_var(iu%ncid, iu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]) ) - if (param%lrhill_present) call check( nf90_put_var(iu%ncid, iu%rhill_varid, self%rhill(j), start=[idslot, tslot]) ) + if (param%lrhill_present) then + call check( nf90_put_var(iu%ncid, iu%rhill_varid, self%rhill(j), start=[idslot, tslot]) ) + end if if (param%lclose) call check( nf90_put_var(iu%ncid, iu%radius_varid, self%radius(j), start=[idslot, tslot]) ) if (param%lrotation) then call check( nf90_put_var(iu%ncid, iu%Ip1_varid, self%Ip(1, j), start=[idslot, tslot]) ) From 638c507611a76fe78106924c9e4e5a06b52b2aa3 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 24 Nov 2021 11:40:18 -0500 Subject: [PATCH 05/10] Removed Intel Advisor stuff from the Makefile.Defines --- Makefile.Defines | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/Makefile.Defines b/Makefile.Defines index 9f002eeda..fcc43b14a 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -45,13 +45,10 @@ COLLRESOLVE_HOME = $(ROOT_DIR)/collresolve/ # DO NOT include in FFLAGS the "-c" option to compile object only # this is done explicitly as needed in the Makefile -ADVIXE_FLAGS = -g -O2 -qopt-report=5 -vecabi=cmdtarget -simd -shared-intel -debug inline-debug-info -DTBB_DEBUG -DTBB_USE_THREADING_TOOLS -xhost -traceback -I$(ADVISOR_2019_DIR)/include/intel64 -parallel-source-info=2 -#Be sure to set the environment variable KMP_FORKJOIN_FRAMES=1 for OpenMP debuging in vtune - IDEBUG = -O0 -init=snan,arrays -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -qopt-matmul -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin STRICTREAL = -fp-model=precise -prec-div -prec-sqrt -assume protect-parens SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -fp-model no-except -fma -PAR = -qopenmp -parallel -parallel-source-info=2 +PAR = -qopenmp -parallel HEAPARR = -heap-arrays 4194304 OPTREPORT = -qopt-report=5 IPRODUCTION = -no-wrap-margin -O3 -qopt-prefetch=0 -qopt-matmul -sox $(PAR) $(SIMDVEC) #$(HEAPARR) @@ -64,7 +61,7 @@ GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporari GPRODUCTION = -O3 -ffree-line-length-none $(GPAR) INCLUDES = -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include -I$(MKLROOT)/include -LINKS = -L$(MKLROOT)/lib/intel64 -L$(NETCDF_FORTRAN_HOME)/lib -L$(ADVISOR_2019_DIR)/lib64 -lswiftest -lnetcdf -lnetcdff -qopt-matmul $(PAR) +LINKS = -L$(MKLROOT)/lib/intel64 -L$(NETCDF_FORTRAN_HOME)/lib -lswiftest -lnetcdf -lnetcdff -qopt-matmul $(PAR) FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) FFLAGS = $(IPRODUCTION) -fp-model=fast From 522cfd669a24151bb62e7a90839d14c91a085eb4 Mon Sep 17 00:00:00 2001 From: Carlisle Wishard Date: Sun, 5 Dec 2021 11:47:02 -0500 Subject: [PATCH 06/10] fixed index of massive body names --- src/discard/discard.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index ad5b73426..d791274ef 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -247,7 +247,7 @@ subroutine discard_pl_tp(tp, system, param) write(idstrj, *) pl%id(j) write(timestr, *) param%t write(*, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & - // " too close to massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstrj)) & + // " too close to massive body " // trim(adjustl(pl%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & // " at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) From 90bfc56e79bfbd96bd4d133e6f9e1db9cc556db3 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 8 Dec 2021 14:01:07 -0500 Subject: [PATCH 07/10] Put back debug Makefile.Defines --- Makefile.Defines | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/Makefile.Defines b/Makefile.Defines index fcc43b14a..90a2ceb58 100644 --- a/Makefile.Defines +++ b/Makefile.Defines @@ -41,14 +41,20 @@ SWIFTEST_HOME = $(ROOT_DIR) USER_MODULES = COLLRESOLVE_HOME = $(ROOT_DIR)/collresolve/ +#NETCDF_FORTRAN_HOME = /home/daminton/local + + # Compiler definitions # DO NOT include in FFLAGS the "-c" option to compile object only # this is done explicitly as needed in the Makefile +ADVIXE_FLAGS = -g -O2 -qopt-report=5 -vecabi=cmdtarget -simd -shared-intel -debug inline-debug-info -DTBB_DEBUG -DTBB_USE_THREADING_TOOLS -xhost -traceback -parallel-source-info=2 +#Be sure to set the environment variable KMP_FORKJOIN_FRAMES=1 for OpenMP debuging in vtune + IDEBUG = -O0 -init=snan,arrays -nogen-interfaces -no-pie -no-ftz -fpe-all=0 -g -traceback -mp1 -qopt-matmul -fp-model strict -fpe0 -debug all -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin STRICTREAL = -fp-model=precise -prec-div -prec-sqrt -assume protect-parens SIMDVEC = -simd -xhost -align all -assume contiguous_assumed_shape -vecabi=cmdtarget -fp-model no-except -fma -PAR = -qopenmp -parallel +PAR = -qopenmp -parallel -parallel-source-info=2 HEAPARR = -heap-arrays 4194304 OPTREPORT = -qopt-report=5 IPRODUCTION = -no-wrap-margin -O3 -qopt-prefetch=0 -qopt-matmul -sox $(PAR) $(SIMDVEC) #$(HEAPARR) @@ -60,11 +66,13 @@ GMEM = -fsanitize-address-use-after-scope -fstack-check -fsanitize=bounds-stri GWARNINGS = -Wall -Warray-bounds -Wimplicit-interface -Wextra -Warray-temporaries GPRODUCTION = -O3 -ffree-line-length-none $(GPAR) -INCLUDES = -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include -I$(MKLROOT)/include -LINKS = -L$(MKLROOT)/lib/intel64 -L$(NETCDF_FORTRAN_HOME)/lib -lswiftest -lnetcdf -lnetcdff -qopt-matmul $(PAR) +INCLUDES = -I$(SWIFTEST_HOME)/include -I$(NETCDF_FORTRAN_HOME)/include -I$(MKLROOT)/include +LINKS = -L$(MKLROOT)/lib/intel64 -L$(NETCDF_FORTRAN_HOME)/lib -lswiftest -lnetcdf -lnetcdff -qopt-matmul #$(PAR) -FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) -FFLAGS = $(IPRODUCTION) -fp-model=fast +FSTRICTFLAGS = $(IDEBUG) #$(SIMDVEC) $(PAR) $(HEAPARR) +FFLAGS = $(IDEBUG) #$(SIMDVEC) $(PAR) $(HEAPARR) +#FSTRICTFLAGS = $(IPRODUCTION) $(STRICTREAL) -g -traceback +#FFLAGS = $(IPRODUCTION) -fp-model=fast -g -traceback FORTRAN = ifort AR = xiar From 56a5b3e08101b5f04d9c03d3f88f4979f5d0b144 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 8 Dec 2021 17:08:42 -0500 Subject: [PATCH 08/10] Added the missing gr_pos_kick step in the close encounter integrations. Moved the vb2vh call to be consistent with the regular helio step sequence --- examples/symba_energy_momentum/param.sun.in | 11 ++++------- .../8pl_16tp_encounters/cb.in | 4 ++-- .../8pl_16tp_encounters/cb.swiftest.in | 4 ++-- .../8pl_16tp_encounters/param.swiftest.in | 5 +++-- src/symba/symba_step.f90 | 8 ++++++-- 5 files changed, 17 insertions(+), 15 deletions(-) diff --git a/examples/symba_energy_momentum/param.sun.in b/examples/symba_energy_momentum/param.sun.in index 351d8adb2..eaf0839f7 100644 --- a/examples/symba_energy_momentum/param.sun.in +++ b/examples/symba_energy_momentum/param.sun.in @@ -1,6 +1,6 @@ !Parameter file for the SyMBA-RINGMOONS test T0 0.0 -TSTOP 3.0e-2 +TSTOP 3.0e-1 DT 1e-3 CB_IN cb.in PL_IN sun.in @@ -8,10 +8,10 @@ TP_IN tp.in IN_TYPE ASCII ISTEP_OUT 1 ISTEP_DUMP 1 -BIN_OUT bin.sun.dat +BIN_OUT bin.sun.nc PARTICLE_OUT particle.sun.dat -OUT_TYPE REAL8 -OUT_FORM XV ! osculating element output +OUT_TYPE NETCDF_DOUBLE +OUT_FORM XVEL ! osculating element output OUT_STAT REPLACE ISTEP_DUMP 1 ! system dump cadence CHK_CLOSE yes ! check for planetary close encounters @@ -19,8 +19,6 @@ CHK_RMIN 0.005 CHK_RMAX 1e2 CHK_EJECT -1.0 ! ignore this check CHK_QMIN -1.0 ! ignore this check -ENC_OUT enc.sun.dat -DISCARD_OUT discard.sun.dat EXTRA_FORCE no ! no extra user-defined forces BIG_DISCARD no ! output all planets if anything discarded RHILL_PRESENT no ! Hill's sphere radii in input file @@ -30,6 +28,5 @@ MU2KG 1.98908e30 TU2S 3.1556925e7 DU2M 1.49598e11 ENERGY yes -ENERGY_OUT energy.sun.out ROTATION yes SEED 8 1230834 2346113 123409874 -123121105 -767545 -534058022 343309814 -12535638 diff --git a/examples/symba_swifter_comparison/8pl_16tp_encounters/cb.in b/examples/symba_swifter_comparison/8pl_16tp_encounters/cb.in index 2df47f957..e3318aec0 100644 --- a/examples/symba_swifter_comparison/8pl_16tp_encounters/cb.in +++ b/examples/symba_swifter_comparison/8pl_16tp_encounters/cb.in @@ -1,5 +1,5 @@ Sun 0.00029591220819207774 0.004650467260962157 -4.7535806948127355e-12 --2.2473967953572827e-18 +0.0 +0.0 diff --git a/examples/symba_swifter_comparison/8pl_16tp_encounters/cb.swiftest.in b/examples/symba_swifter_comparison/8pl_16tp_encounters/cb.swiftest.in index 2df47f957..e3318aec0 100644 --- a/examples/symba_swifter_comparison/8pl_16tp_encounters/cb.swiftest.in +++ b/examples/symba_swifter_comparison/8pl_16tp_encounters/cb.swiftest.in @@ -1,5 +1,5 @@ Sun 0.00029591220819207774 0.004650467260962157 -4.7535806948127355e-12 --2.2473967953572827e-18 +0.0 +0.0 diff --git a/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in b/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in index bbf9e7089..7f4d21c9b 100644 --- a/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in +++ b/examples/symba_swifter_comparison/8pl_16tp_encounters/param.swiftest.in @@ -32,6 +32,7 @@ FRAGMENTATION NO ROTATION NO TIDES NO ENERGY NO -GR NO +GR YES GMTINY 1e-12 -ENCOUNTER_CHECK SORTSWEEP +ENCOUNTER_CHECK TRIANGULAR +INTERACTION_LOOPS TRIANGULAR diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index eac112dcd..330f47a56 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -71,22 +71,26 @@ module subroutine symba_step_interp_system(self, param, t, dt) call pl%vh2vb(cb) call pl%lindrift(cb, dth, lbeg=.true.) call pl%kick(system, param, t, dth, lbeg=.true.) + if (param%lgr) call pl%gr_pos_kick(param, dth) call pl%drift(system, param, dt) call tp%vh2vb(vbcb = -cb%ptbeg) call tp%lindrift(cb, dth, lbeg=.true.) call tp%kick(system, param, t, dth, lbeg=.true.) + if (param%lgr) call tp%gr_pos_kick(param, dth) call tp%drift(system, param, dt) call system%recursive_step(param, t, 0) call pl%kick(system, param, t, dth, lbeg=.false.) - call pl%vb2vh(cb) + if (param%lgr) call pl%gr_pos_kick(param, dth) call pl%lindrift(cb, dth, lbeg=.false.) + call pl%vb2vh(cb) call tp%kick(system, param, t, dth, lbeg=.false.) - call tp%vb2vh(vbcb = -cb%ptend) + if (param%lgr) call tp%gr_pos_kick(param, dth) call tp%lindrift(cb, dth, lbeg=.false.) + call tp%vb2vh(vbcb = -cb%ptend) end select end select end select From 27ae7369614f412140fb02040973e9b112ccd8b2 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 9 Dec 2021 11:48:19 -0500 Subject: [PATCH 09/10] Added in override to gr_pos_kick for SyMBA, which allows the position kick operation to be masked during recursive steps so only bodies at the current recursion depth are kicked --- src/helio/helio_gr.f90 | 23 +++++++----- src/helio/helio_step.f90 | 12 +++--- src/modules/fraggle_classes.f90 | 2 +- src/modules/helio_classes.f90 | 24 ++++++------ src/modules/symba_classes.f90 | 23 ++++++++++++ src/modules/whm_classes.f90 | 18 +++++---- src/symba/symba_gr.f90 | 66 +++++++++++++++++++++++++++++++++ src/symba/symba_step.f90 | 17 ++++++--- src/whm/whm_gr.f90 | 22 ++++++----- src/whm/whm_step.f90 | 8 ++-- 10 files changed, 160 insertions(+), 55 deletions(-) create mode 100644 src/symba/symba_gr.f90 diff --git a/src/helio/helio_gr.f90 b/src/helio/helio_gr.f90 index 5d2936324..1ff8eb5d5 100644 --- a/src/helio/helio_gr.f90 +++ b/src/helio/helio_gr.f90 @@ -56,7 +56,7 @@ module pure subroutine helio_gr_kick_getacch_tp(self, param) end subroutine helio_gr_kick_getacch_tp - module pure subroutine helio_gr_p4_pl(self, param, dt) + module pure subroutine helio_gr_p4_pl(self, system, param, dt) !! author: David A. Minton !! !! Position kick to massive bodies due to p**4 term in the post-Newtonian correction @@ -65,11 +65,12 @@ module pure subroutine helio_gr_p4_pl(self, param, dt) !! Adapted from David A. Minton's Swifter routine routine gr_helio_p4.f90 implicit none ! Arguments - class(helio_pl), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(helio_pl), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i + integer(I4B) :: i if (self%nbody == 0) return @@ -82,7 +83,8 @@ module pure subroutine helio_gr_p4_pl(self, param, dt) return end subroutine helio_gr_p4_pl - module pure subroutine helio_gr_p4_tp(self, param, dt) + + module pure subroutine helio_gr_p4_tp(self, system, param, dt) !! author: David A. Minton !! !! Position kick to test particles due to p**4 term in the post-Newtonian correction @@ -91,11 +93,12 @@ module pure subroutine helio_gr_p4_tp(self, param, dt) !! Adapted from David A. Minton's Swifter routine routine gr_helio_p4.f90 implicit none ! Arguments - class(helio_tp), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(helio_tp), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i + integer(I4B) :: i if (self%nbody == 0) return diff --git a/src/helio/helio_step.f90 b/src/helio/helio_step.f90 index 039884596..fb14c9e75 100644 --- a/src/helio/helio_step.f90 +++ b/src/helio/helio_step.f90 @@ -35,7 +35,7 @@ module subroutine helio_step_pl(self, system, param, t, dt) implicit none ! Arguments class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Stepsize @@ -54,10 +54,10 @@ module subroutine helio_step_pl(self, system, param, t, dt) end if call pl%lindrift(cb, dth, lbeg=.true.) call pl%kick(system, param, t, dth, lbeg=.true.) - if (param%lgr) call pl%gr_pos_kick(param, dth) + if (param%lgr) call pl%gr_pos_kick(system, param, dth) call pl%drift(system, param, dt) call pl%kick(system, param, t + dt, dth, lbeg=.false.) - if (param%lgr) call pl%gr_pos_kick(param, dth) + if (param%lgr) call pl%gr_pos_kick(system, param, dth) call pl%lindrift(cb, dth, lbeg=.false.) call pl%vb2vh(cb) end select @@ -78,7 +78,7 @@ module subroutine helio_step_tp(self, system, param, t, dt) implicit none ! Arguments class(helio_tp), intent(inout) :: self !! Helio test particle data structure - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Stepsize @@ -97,10 +97,10 @@ module subroutine helio_step_tp(self, system, param, t, dt) end if call tp%lindrift(cb, dth, lbeg=.true.) call tp%kick(system, param, t, dth, lbeg=.true.) - if (param%lgr) call tp%gr_pos_kick(param, dth) + if (param%lgr) call tp%gr_pos_kick(system, param, dth) call tp%drift(system, param, dt) call tp%kick(system, param, t + dt, dth, lbeg=.false.) - if (param%lgr) call tp%gr_pos_kick(param, dth) + if (param%lgr) call tp%gr_pos_kick(system, param, dth) call tp%lindrift(cb, dth, lbeg=.false.) call tp%vb2vh(vbcb = -cb%ptend) end select diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index cd648c04b..042ff5a78 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -150,7 +150,7 @@ module subroutine fraggle_placeholder_step(self, system, param, t, dt) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none class(fraggle_fragments), intent(inout) :: self !! Helio massive body particle object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Stepsiz diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index 2d15565b2..78a4bdc34 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -122,20 +122,22 @@ module pure subroutine helio_gr_kick_getacch_tp(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine helio_gr_kick_getacch_tp - module pure subroutine helio_gr_p4_pl(self, param, dt) - use swiftest_classes, only : swiftest_parameters + module pure subroutine helio_gr_p4_pl(self, system, param, dt) + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system implicit none - class(helio_pl), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(helio_pl), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size end subroutine helio_gr_p4_pl - module pure subroutine helio_gr_p4_tp(self, param, dt) - use swiftest_classes, only : swiftest_parameters + module pure subroutine helio_gr_p4_tp(self, system, param, dt) + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system implicit none - class(helio_tp), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(helio_tp), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size end subroutine helio_gr_p4_tp module subroutine helio_kick_getacch_pl(self, system, param, t, lbeg) @@ -191,7 +193,7 @@ module subroutine helio_step_pl(self, system, param, t, dt) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none class(helio_pl), intent(inout) :: self !! Helio massive body particle object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Stepsize diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 4a584eda9..99c527143 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -75,6 +75,7 @@ module symba_classes procedure :: discard => symba_discard_pl !! Process massive body discards procedure :: drift => symba_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level procedure :: encounter_check => symba_encounter_check_pl !! Checks if massive bodies are going through close encounters with each other + procedure :: gr_pos_kick => symba_gr_p4_pl !! Position kick due to p**4 term in the post-Newtonian correction procedure :: accel_int => symba_kick_getacch_int_pl !! Compute direct cross (third) term heliocentric accelerations of massive bodiess, with no mutual interactions between bodies below GMTINY procedure :: accel => symba_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies procedure :: setup => symba_setup_pl !! Constructor method - Allocates space for the input number of bodies @@ -113,6 +114,7 @@ module symba_classes contains procedure :: drift => symba_drift_tp !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level procedure :: encounter_check => symba_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body + procedure :: gr_pos_kick => symba_gr_p4_tp !! Position kick due to p**4 term in the post-Newtonian correction procedure :: accel => symba_kick_getacch_tp !! Compute heliocentric accelerations of test particles procedure :: setup => symba_setup_tp !! Constructor method - Allocates space for the input number of bodies procedure :: append => symba_util_append_tp !! Appends elements from one structure to another @@ -272,6 +274,7 @@ module subroutine symba_drift_tp(self, system, param, dt) end subroutine symba_drift_tp module function symba_encounter_check_pl(self, param, system, dt, irec) result(lany_encounter) + use swiftest_classes, only : swiftest_nbody_system implicit none class(symba_pl), intent(inout) :: self !! SyMBA test particle object class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters @@ -282,6 +285,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l end function symba_encounter_check_pl module function symba_encounter_check(self, param, system, dt, irec) result(lany_encounter) + use swiftest_classes, only : swiftest_parameters implicit none class(symba_encounter), intent(inout) :: self !! SyMBA pl-pl encounter list object class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters @@ -292,6 +296,7 @@ module function symba_encounter_check(self, param, system, dt, irec) result(lany end function symba_encounter_check module function symba_encounter_check_tp(self, param, system, dt, irec) result(lany_encounter) + use swiftest_classes, only : swiftest_parameters implicit none class(symba_tp), intent(inout) :: self !! SyMBA test particle object class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters @@ -301,6 +306,24 @@ module function symba_encounter_check_tp(self, param, system, dt, irec) result(l logical :: lany_encounter !! Returns true if there is at least one close encounter end function symba_encounter_check_tp + module pure subroutine symba_gr_p4_pl(self, system, param, dt) + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system + implicit none + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size + end subroutine symba_gr_p4_pl + + module pure subroutine symba_gr_p4_tp(self, system, param, dt) + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system + implicit none + class(symba_tp), intent(inout) :: self !! SyMBA test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size + end subroutine symba_gr_p4_tp + module function symba_collision_casedisruption(system, param, colliders, frag) result(status) use fraggle_classes, only : fraggle_colliders, fraggle_fragments implicit none diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index 066975d41..a7cd2f49c 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -170,20 +170,22 @@ module pure subroutine whm_gr_kick_getacch_tp(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine whm_gr_kick_getacch_tp - module pure subroutine whm_gr_p4_pl(self, param, dt) + module pure subroutine whm_gr_p4_pl(self, system, param, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(whm_pl), intent(inout) :: self !! WHM massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size end subroutine whm_gr_p4_pl - module pure subroutine whm_gr_p4_tp(self, param, dt) + module pure subroutine whm_gr_p4_tp(self, system, param, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(whm_tp), intent(inout) :: self !! WHM test particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(whm_tp), intent(inout) :: self !! WHM test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size end subroutine whm_gr_p4_tp !> Reads WHM massive body object in from file diff --git a/src/symba/symba_gr.f90 b/src/symba/symba_gr.f90 new file mode 100644 index 000000000..1b0246aa8 --- /dev/null +++ b/src/symba/symba_gr.f90 @@ -0,0 +1,66 @@ +submodule(symba_classes) s_symba_gr + use swiftest +contains + + module pure subroutine symba_gr_p4_pl(self, system, param, dt) + !! author: David A. Minton + !! + !! Position kick to massive bodies due to p**4 term in the post-Newtonian correction + !! Based on Saha & Tremaine (1994) Eq. 28 + !! + !! Adapted from David A. Minton's Swifter routine routine gr_symba_p4.f90 + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size + ! Internals + integer(I4B) :: i + + if (self%nbody == 0) return + + associate(pl => self, npl => self%nbody) + select type(system) + class is (symba_nbody_system) + do concurrent(i = 1:npl, pl%lmask(i) .and. pl%levelg(i) == system%irec ) + call gr_p4_pos_kick(param, pl%xh(:, i), pl%vb(:, i), dt) + end do + end select + end associate + + return + end subroutine symba_gr_p4_pl + + + module pure subroutine symba_gr_p4_tp(self, system, param, dt) + !! author: David A. Minton + !! + !! Position kick to test particles due to p**4 term in the post-Newtonian correction + !! Based on Saha & Tremaine (1994) Eq. 28 + !! + !! Adapted from David A. Minton's Swifter routine routine gr_symba_p4.f90 + implicit none + ! Arguments + class(symba_tp), intent(inout) :: self !! SyMBA test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size + ! Internals + integer(I4B) :: i + + if (self%nbody == 0) return + + associate(tp => self, ntp => self%nbody) + select type(system) + class is (symba_nbody_system) + do concurrent(i = 1:ntp, tp%lmask(i) .and. tp%levelg(i) == system%irec) + call gr_p4_pos_kick(param, tp%xh(:, i), tp%vb(:, i), dt) + end do + end select + end associate + + return + end subroutine symba_gr_p4_tp + +end submodule s_symba_gr diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 330f47a56..f10850e7f 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -71,24 +71,24 @@ module subroutine symba_step_interp_system(self, param, t, dt) call pl%vh2vb(cb) call pl%lindrift(cb, dth, lbeg=.true.) call pl%kick(system, param, t, dth, lbeg=.true.) - if (param%lgr) call pl%gr_pos_kick(param, dth) + if (param%lgr) call pl%gr_pos_kick(system, param, dth) call pl%drift(system, param, dt) call tp%vh2vb(vbcb = -cb%ptbeg) call tp%lindrift(cb, dth, lbeg=.true.) call tp%kick(system, param, t, dth, lbeg=.true.) - if (param%lgr) call tp%gr_pos_kick(param, dth) + if (param%lgr) call tp%gr_pos_kick(system, param, dth) call tp%drift(system, param, dt) call system%recursive_step(param, t, 0) call pl%kick(system, param, t, dth, lbeg=.false.) - if (param%lgr) call pl%gr_pos_kick(param, dth) + if (param%lgr) call pl%gr_pos_kick(system, param, dth) call pl%lindrift(cb, dth, lbeg=.false.) call pl%vb2vh(cb) call tp%kick(system, param, t, dth, lbeg=.false.) - if (param%lgr) call tp%gr_pos_kick(param, dth) + if (param%lgr) call tp%gr_pos_kick(system, param, dth) call tp%lindrift(cb, dth, lbeg=.false.) call tp%vb2vh(vbcb = -cb%ptend) end select @@ -185,7 +185,10 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci) call plplenc_list%kick(system, dth, irecp, -1) call pltpenc_list%kick(system, dth, irecp, -1) end if - + if (param%lgr) then + call pl%gr_pos_kick(system, param, dth) + call tp%gr_pos_kick(system, param, dth) + end if call pl%drift(system, param, dtl) call tp%drift(system, param, dtl) @@ -198,6 +201,10 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci) call plplenc_list%kick(system, dth, irecp, -1) call pltpenc_list%kick(system, dth, irecp, -1) end if + if (param%lgr) then + call pl%gr_pos_kick(system, param, dth) + call tp%gr_pos_kick(system, param, dth) + end if if (param%lclose) then lplpl_collision = plplenc_list%collision_check(system, param, t+dtl, dtl, ireci) diff --git a/src/whm/whm_gr.f90 b/src/whm/whm_gr.f90 index bbe6ca30b..12ae82a35 100644 --- a/src/whm/whm_gr.f90 +++ b/src/whm/whm_gr.f90 @@ -61,7 +61,7 @@ module pure subroutine whm_gr_kick_getacch_tp(self, param) end subroutine whm_gr_kick_getacch_tp - module pure subroutine whm_gr_p4_pl(self, param, dt) + module pure subroutine whm_gr_p4_pl(self, system, param, dt) !! author: David A. Minton !! !! Position kick to massive bodies due to p**4 term in the post-Newtonian correction @@ -70,11 +70,12 @@ module pure subroutine whm_gr_p4_pl(self, param, dt) !! Adapted from David A. Minton's Swifter routine routine gr_whm_p4.f90 implicit none ! Arguments - class(whm_pl), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(whm_pl), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i + integer(I4B) :: i associate(pl => self, npl => self%nbody) if (npl == 0) return @@ -87,7 +88,7 @@ module pure subroutine whm_gr_p4_pl(self, param, dt) end subroutine whm_gr_p4_pl - module pure subroutine whm_gr_p4_tp(self, param, dt) + module pure subroutine whm_gr_p4_tp(self, system, param, dt) !! author: David A. Minton !! !! Position kick to test particles due to p**4 term in the post-Newtonian correction @@ -96,11 +97,12 @@ module pure subroutine whm_gr_p4_tp(self, param, dt) !! Adapted from David A. Minton's Swifter routine routine gr_whm_p4.f90 implicit none ! Arguments - class(whm_tp), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(whm_tp), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i + integer(I4B) :: i associate(tp => self, ntp => self%nbody) if (ntp == 0) return diff --git a/src/whm/whm_step.f90 b/src/whm/whm_step.f90 index d194e2c02..244ebbbb1 100644 --- a/src/whm/whm_step.f90 +++ b/src/whm/whm_step.f90 @@ -50,9 +50,9 @@ module subroutine whm_step_pl(self, system, param, t, dt) dth = 0.5_DP * dt call pl%kick(system, param, t, dth,lbeg=.true.) call pl%vh2vj(cb) - if (param%lgr) call pl%gr_pos_kick(param, dth) + if (param%lgr) call pl%gr_pos_kick(system, param, dth) call pl%drift(system, param, dt) - if (param%lgr) call pl%gr_pos_kick(param, dth) + if (param%lgr) call pl%gr_pos_kick(system, param, dth) call pl%j2h(cb) call pl%kick(system, param, t + dt, dth, lbeg=.false.) end associate @@ -85,9 +85,9 @@ module subroutine whm_step_tp(self, system, param, t, dt) associate(tp => self, cb => system%cb, pl => system%pl) dth = 0.5_DP * dt call tp%kick(system, param, t, dth, lbeg=.true.) - if (param%lgr) call tp%gr_pos_kick(param, dth) + if (param%lgr) call tp%gr_pos_kick(system, param, dth) call tp%drift(system, param, dt) - if (param%lgr) call tp%gr_pos_kick(param, dth) + if (param%lgr) call tp%gr_pos_kick(system, param, dth) call tp%kick(system, param, t + dt, dth, lbeg=.false.) end associate end select From bf82f65ee7f003645e2b38dbff1ad20be5809476 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 10 Dec 2021 14:35:18 -0500 Subject: [PATCH 10/10] Updated Notebook with fixes to the angles from radians to degree, and removed the Swifter comparisons --- .../helio_gr_test/swiftest_relativity.ipynb | 67 +++++-------------- 1 file changed, 18 insertions(+), 49 deletions(-) diff --git a/examples/helio_gr_test/swiftest_relativity.ipynb b/examples/helio_gr_test/swiftest_relativity.ipynb index db7a4926e..db2b358ed 100644 --- a/examples/helio_gr_test/swiftest_relativity.ipynb +++ b/examples/helio_gr_test/swiftest_relativity.ipynb @@ -21,56 +21,33 @@ "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": [ + "env: HDF5_USE_FILE_LOCKING=FALSE\n", "Reading Swiftest file param.swiftest.in\n", - "Reading in time 1.000e+03\n", - "Creating Dataset\n", + "\n", + "Creating Dataset from NetCDF file\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" + "%env HDF5_USE_FILE_LOCKING=FALSE\n", + "sim = swiftest.Simulation(param_file=\"param.swiftest.in\")" ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ - "swifterdat['varpi'] = swifterdat['omega'] + swifterdat['capom']\n", - "swiftestdat['varpi'] = swiftestdat['omega'] + swiftestdat['capom']" + "sim.ds = sim.ds.swap_dims({\"id\" : \"name\"})\n", + "sim.ds['varpi'] = sim.ds['omega'] + sim.ds['capom']" ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -84,29 +61,27 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ - "varpiswiftest = swiftestdat['varpi'].sel(id=1) * 180.0 / np.pi\n", - "varpiswifter = swifterdat['varpi'].sel(id=1) * 180.0 / np.pi\n", - "tsim = swiftestdat['time']" + "tsim = sim.ds['time']\n", + "varpiswiftest = sim.ds.sel(name='Mercury')['varpi']" ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 6, "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, + "execution_count": 9, "metadata": {}, "outputs": [ { @@ -115,16 +90,13 @@ "text": [ "Mean precession rate for Mercury long. peri. (arcsec/100 y)\n", "JPL Horizons : 571.3210506300043\n", - "Swifter GR : 571.6183105524942\n", - "Swiftest GR : 571.5701928436062\n", - "Obs - Swifter : -0.2972599224899675\n", - "Obs - Swiftest : -0.2491422136018948\n", - "Swiftest - Swifter: -0.04811770888807132\n" + "Swiftest GR : 571.519655556524\n", + "Obs - Swiftest : -0.1986049265197031\n" ] }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEGCAYAAAB2EqL0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA45klEQVR4nO3deZzN1RvA8c9zZ7XvLZKtpKyDsa9ZkiUKhZKlhGijjRbh91OS9pJQKLJniUQ7ShhLdj/Zaigh65gxM/c+vz/upTEG9zIz31me9+t1X3O/536X51zM45zz/Z4jqooxxhhzKS6nAzDGGJM5WMIwxhjjF0sYxhhj/GIJwxhjjF8sYRhjjPFLsNMBpKXChQtryZIlnQ7DGGMyjTVr1hxS1SIpfZalE0bJkiWJiopyOgxjjMk0RGTvhT6zLiljjDF+sYRhjDHGL5YwjDHG+CVLj2GkJCEhgejoaOLi4pwOJdsJDw+nWLFihISEOB2KMeYyZLuEER0dTZ48eShZsiQi4nQ42YaqcvjwYaKjoylVqpTT4RhjLkO265KKi4ujUKFClizSmYhQqFAha9kZk4llu4QBWLJwiH3vxmRu2TJhGGNMVjXt23cYOaUnHrc71c9tCcMBuXPnZs+ePeTIkYOIiAjKlStHnz598Hg87NmzhwoVKlz0+CFDhjBq1KhzykqWLMmhQ4cCiqNly5YcPXo00PCNMRlQ/OkTPDmhIcOjx7Ek9heOnjyc6tfIdoPeGckNN9zA+vXrSUxMpHHjxsydO5eqVaum+XVVFVXlyy+/TPNrGWPS3uJvxzB55zjWh8UTeTI/XesMpGC+q1L9OtbCyACCg4OpU6cOv/32W6qc74033qBChQpUqFCBt956C4A9e/Zwyy230LdvX6pWrcoff/xxtlUyZswYIiIiiIiIoFSpUtx6660ATJ06lYoVK1KhQgWeffbZs+fPnTs3zz//PJUrV6ZWrVocOHAAgJkzZ1KhQgUqV65MgwYNUqUuxpgL27LrZx4bX5+not9nfVg8jeOvZUzvH7i1Sqs0uV62bmEM/WIzW/YfT9Vzliual5fuKB/QMadOneLbb79l2LBhfh/z5ptvMnny5LPb+/fvB2DNmjVMmDCBlStXoqrUrFmThg0bUqBAAbZv386ECRMYPXr0Oefq06cPffr0ISEhgcaNGzNgwAD279/Ps88+y5o1ayhQoAC33XYbc+fO5c477yQmJoZatWoxfPhwnnnmGcaNG8cLL7zAsGHDWLx4Mdddd511dRmThk4nxPJV1Eze3DKSf4Kh8vEiPNSgLw0qtkdcadcOsBaGg3bu3ElERAR169alVatWtGjRwu9j+/fvz/r168++ihYtCsDy5cu56667yJUrF7lz56Zdu3YsW7YMgBIlSlCrVq0LnvPxxx+ncePG3HHHHaxevZpGjRpRpEgRgoODue+++1i6dCkAoaGhtG7dGoBq1aqxZ88eAOrWrUv37t0ZN24c7jQYcDPGwI/rF9Luk9q88L/XOBws9M17D691/ZKGle9O02QB2byFEWhLILWdGcNITap6wc9y5cp1wc8mTpzI3r17ee+99y55npCQkLO3yAYFBZGYmAjAmDFjWLlyJQsXLiQiIoL169dTqFChy6mGMSaZU3HHmf/TJIZHj4VgqHUyP41L1aNzy8HpFoO1MLKYBg0aMHfuXE6dOkVMTAxz5syhfv36Fz1mzZo1jBo1ismTJ+Py/Q+lZs2a/Pjjjxw6dAi3283UqVNp2LDhRc+zc+dOatasybBhwyhcuDB//PFHqtXLmOxKVdm8YxENptZhePRYglW5x9WGYV2W0LnlK+kaS7q1MESkLDA9SVFpYDDwPTAGyA3sAe5T1fMGFkTkduBtIAgYr6oj0jrmtJCYmEhYWNhF99m+fTvFihU7u/3mm29y9913+3X+qlWr0r17d2rUqAFAz549qVKlytluo5S89957/PPPP2cHuyMjIxk/fjyvvPIKt956K6pKy5Ytadu27UWv/fTTT7Njxw5UlSZNmlC5cmW/YjbGpMyjHvqPa8V3YdGEoTQ9mpc21R7m1nr3OxKPXKzrIc0uKhIE7ANqArOAp1T1RxF5ACilqi+msP//gGZANLAa6KyqWy52ncjISE2+gNLWrVu55ZZbUq0ugfr111956KGHWLVqlWMxOMnp79+YzOKTuY8w/p8fOBIk5PB4aBtch+e6jE3zGRNEZI2qRqb0mVNjGE2Anaq619fyWOor/xpYDLyYbP8awG+qugtARKYBbYGLJoyMZsyYMbzzzjtnb3U1xpjk9vyxlncXvcCSsD/Ig1L3aFFeeWAuBfLkdDo0xxJGJ2Cq7/0moA0wD7gbuD6F/a8DknaIR+NtnZxHRHoBvQCKFy+eSuGmjjO3rxpjTHKJiW4++/K/jD00g2NhLoJUGVZ+KHWr3EWO0CCnwwMcGPQWkVC8CWKmr+gBoJ+IrAHyAPEpHZZCWYp9aao6VlUjVTWySJEU1zE3xpgMZcP2BTwzoRWvHZnFsSAXzWJrMqP+OzSt2SHDJAtwpoXRAlirqgcAVHUbcBuAiNwEpPSIYjTntjyKAfvTOE5jjElTiW4Pm7Z/SbdVA/GECtcmKF2K9aNNgwfJnzPU6fDO40TC6My/3VGIyFWq+reIuIAX8N4xldxqoIyIlMI7WN4JuDc9gjXGmNSmHg9HTkTzn+mP8U3YThChaUxFBrQdzPXX3ex0eBeUrglDRHLivdOpd5LiziLSz/f+c2CCb9+ieG+fbamqiSLyCN4B8SDgY1XdnI6hG2NMqvngi958cPQXCIPyp11UL3gX/e9/CZcrY68Zk65jGKp6SlULqeqxJGVvq+pNvtdA9d3nq6r7VbVlkv2+9O1zg6oOT8+4U9vw4cMpX748lSpVIiIigpUrV/p13ODBg/nmm28AWLZsGeXLlyciIoIVK1akysyzBw4c4N5776V06dJUq1aN2rVrM2fOHAB++OEH8uXLR5UqVbj55pt56qmnrvh6xmQ3v/25gX7jOjPxnxUAdEisycRuK3myw5AMnywgm08N4oQVK1awYMEC1q5dS1hYGIcOHSI+PqVx/vMlnZxwypQpPPXUU/To0YOJEycSFRVFy5YtL3L0uRITEwkO/vePX1W588476datG5999hkAe/fuZf78+Wf3qV+/PgsWLCA2NpYqVapw1113UbduXb+vaUx2lZAYz4BPW7BMDuD2jVW8WmEEt1Zv7XRoAbGEkc7+/PNPChcufPZp78KFCwOwatUqRowYweeff868efPo1KkTx44dw+PxUK5cOXbt2kX37t1p3bo1R48eZcaMGSxevJglS5bw008/ERsby/Llyxk0aBCtW7fm0UcfZePGjSQmJjJkyBDatm3LxIkTWbhwIXFxccTExPDdd9+djeu7774jNDT0nNt+S5QowaOPPnpeHc4s/LRv3740/raMyfwW/vQGc7bOZ2XIYa5PcNMoZ2ee7vI8EpRx7n7yV/ZOGIsGwl8bU/ec11SEFheeteS2225j2LBh3HTTTTRt2pSOHTvSsGFDqlatyrp16wBvd1OFChVYvXo1iYmJ1Kx57iMnPXv2ZPny5bRu3ZoOHTqcbWGcmTjwueeeo3Hjxnz88cccPXqUGjVq0LRpU8DbwtmwYQMFCxY855ybN2/2e/GmI0eOsGPHDlvzwpiL+POfI3z81SCmJfwEIRAZF8K7968md85wp0O7bDb5YDrLnTs3a9asYezYsRQpUoSOHTsyceJEgoODufHGG9m6dSurVq1iwIABLF26lGXLll1y8sDklixZwogRI4iIiKBRo0bExcXx+++/A9CsWbPzkkVK+vXrR+XKlalevfrZsmXLllGpUiWuueYaWrduzTXXXBNY5Y3JDjwe9u79mf4zmjIt4SeuTnAz+KrujH9wRaZOFpDdWxgXaQmkpaCgIBo1akSjRo2oWLEikyZNonv37tSvX59FixYREhJC06ZN6d69O263+7z1uy9FVZk9ezZly5Y9p3zlypUXnOK8fPnyzJ49++z2+++/z6FDh4iM/HdKmTNjGP/73/+oV68ed911FxEREQHFZkxWph4PL39yG9PkAOSA2+Ju4sEWb1GueEoTWGQ+1sJIZ9u3b2fHjh1nt9evX0+JEiUA79Tkb731FrVr16ZIkSIcPnyYbdu2Ub78xdftyJMnDydOnDi73bx5c959992za1qc6eq6mMaNGxMXF8cHH3xwtuzUqVMp7nvTTTcxaNAgXn311Uue15jsYsW6OXQdU4tpcoDKsSG8UOR+Xu89O8skC7CEke5OnjxJt27dKFeuHJUqVWLLli0MGTIE8K5BceDAgbNjA5UqVaJSpUqXnJ3y1ltvZcuWLURERDB9+nRefPFFEhISqFSpEhUqVODFF5PP5Xg+EWHu3Ln8+OOPlCpViho1atCtW7cLJoU+ffqwdOlSdu/eHdgXYEwWc+LkP/xnUjseX/8CW3OcoqW7BB90+4WOLZ9xOrRU58j05uklI05vnt3Z92+yki9XfMjILe9yOFgoFxNOr1ojaVL1VqfDuiIZcXpzY4zJtGLi4hk9txufnN4EwULH0Fr06/A+BXJlvPmfUpMlDGOM8ZPH42bWsjcZ/dsEDge7qBgTSrdqT9C8pjMr4KU3SxjGGOMHTTjN05NuZUnICQh20dFVjztaDaFyiaudDi3dWMIwxphL+OTLV5gbPYUdYUK5uNx0K/8QLes94HRY6c4ShjHGXMC67YsYtPxZ9gUrhAl14osyost8CuQKczo0R1jCMMaYZE7EHubNeYP46tTPhKDUPViSAR3e5qaSNzodmqPsOQwHZKTpzY8ePcro0aMv+LlNeW6ymykLh9Fpyq3MPL2CE0HC/QWaM3rAgmyfLMASRrpLOr35hg0b+Oabb7j+ev+eBB02bNjZSQTPTG++fv16tm/fniYJ48yU5w0aNGDXrl2sWbOGadOmER0dfXaf+vXrs27dOtatW8eCBQv46aefLisOY5y2Zet8RkzrwYhDM/k9RKl3rBhfNHiXnu3eyBRrVaSHdEsYIlJWRNYneR0XkSdEJEJEfvGVRYlIjQsc319ENovIJhGZKiKZchavlKY3L1q0KKtWraJdu3YAzJs3jxw5chAfH09cXBylS5cGoHv37syaNYvx48czY8YMhg0bRufOnRk8eDDTp08/+6R3TEwMDzzwANWrV6dKlSrMmzcP8M5IW6NGDSIiIqhUqRI7duxg4MCB7Ny5k4iICJ5++ulzYrUpz012kJjo5pufp9H1l+eYcjqKHB4PTxZ7iSHd51CyVCOnw8tQ0m0MQ1W3AxEAIhKEd23uOcA4YKiqLhKRlsBIoFHSY0XkOuAxoJyqxorIDLzrek+8kpheXfUq2/7ZdiWnOM/NBW/m2RrPXvBzJ6c3HzNmDI8//jj33Xcf8fHxuN1uRowYwaZNm1i/fv15sdqU5ybLizvG0M8eZm7QRnAJzY5fT5tq99CoTgenI8uQnBr0bgLsVNW9IqJAXl95PmD/BY4JBnKISAKQ8yL7ZWhnpjdftmwZ33//PR07dmTEiBF07949xenN3W73ZU1vPn/+/LOz3J6Z3rx27doMHz6c6Oho2rVrR5kyZQI6b79+/Vi+fDmhoaGsXr0a+HfK8+3btzNw4ECb8txkGp8veoaX/l4EvnWM7g+tyxMPjyY02HrqL8SphNEJmOp7/wSwWERG4e0iq5N8Z1Xd5/v8dyAWWKKqS1I6sYj0AnoBFC9e/KJBXKwlkJacmt78lltuoWbNmixcuJDmzZszfvz4s91dKbEpz01WdCr2GE9Pas/SHAcAaHjqeh5r/R43XX/hfwvGK91TqYiEAm2Amb6ih4H+qno90B/4KIVjCgBtgVJAUSCXiHRJ6fyqOlZVI1U1skiRImlRhSvi5PTmu3btonTp0jz22GO0adOGDRs2nHdsUjbluclqPv3iadpNqcPSHAfI4/YwuEQf3nv4S0sWfnKi7dUCWKuqB3zb3YDPfe9nAikNejcFdqvqQVVN8O1/XkskM3ByevPp06dToUIFIiIi2LZtG127dqVQoULUrVuXChUqnDfobVOem6xicdQkeo9tysh/vmJfiIu6x8uwtNtG7m7Uz+nQMpV0n95cRKYBi1V1gm97K/Cwqv4gIk2AkapaLdkxNYGPgep4u6QmAlGq+u7FrmXTm2c89v2b9JTg9jD7+3G8Gv0uiSJclejhlfqzqViyDDlCg5wOL0PKMNObi0hOoBnQO0nxQ8DbIhIMxOEbfxCRosB4VW2pqitFZBawFkgE1gFj0zN2Y0zm4fa4mb38Zb7dvpSfQ/8iGOgUcg+PtO9LvtyFnA4v00rXhKGqp4BCycqWA9VS2Hc/0DLJ9kvAS2kdozEmc3N7lPHfDua9/fMhFCJig2lR5nk6NGxnd0BdoWw5l5SqXnJcwKS+rLy6o3GeJzGefUd28+KsPqzLcZDrEt3cmbM13e4ZSo4cOZ0OL0vIdgkjPDycw4cPU6hQIUsa6UhVOXz4MOHhmfIBfZPBqcfDI5MasCw4BnJCrdhcdKk+hIbVWjgdWpaS7RJGsWLFiI6O5uDBg06Hku2Eh4dTrFgxp8MwWcyCH0Ywd/scVoafokS8h1YFOvNwtxecDitLynYJIyQkhFKlSjkdhjHmCh05foR35/RlJpsgHGqezs273b4nR5i1YtNKtksYxphMTpWodZ/y9uq3WR8eT7lTLnpHvknj6o2djizLs4RhjMk0TscnMOyzVsyXPyEcWrkrMKDLBK7KY62K9GAJwxiTKWzbvZIhS/qwOTyR6jGh3HXzg9zRuK/TYWUrljCMMRnasWP7+c/nXVnsOkB4qNJGKvJ0t4nkz6brajvJEoYxJsOa+/0bfLhzIvuDPVQ9npe2lR6lXcPOToeVbVnCMMZkOLv2bWP8N0/zBXsgBDqHVGbQI5Pt2SmHWcIwxmQYsXHHeHvBY0yJWQtAxPEcPNvqEyqUvtnhyAxYwjDGZBBHT8Uy8LPm/BQSA0DvXHdya5MnKV8sv7OBmbMsYRhjHHXkyC6+WDaRdw9/TlyIUCemED1q9qZWNRuryGgsYRhjHLNtx2K6LH+S0y4hFGjjqcwjncdybQGbLDAjsoRhjEl3sadPMPDTu/kpKBoXSpNDZejSpD+REQ2cDs1chCUMY0y6+uyLAbx3cDEnQlyEKHTOUYsne3zkdFjGD+mWMESkLDA9SVFpYDDwAzAGCMe7ml5fVV2VwvH5gfFABUCBB1R1RdpGbYxJLet/ncFnayaxKOh3CHJR49i1vNFzIflyhjgdmvFTuiUMVd0ORACISBCwD5gDjAOGquoiEWkJjAQapXCKt4GvVLWDiIQC1slpTCagHg+Tv3yV9w5O5lSQC1HllZsG0aDqPeQJt2SRmTjVJdUE2Kmqe0VEgby+8nzA/uQ7i0heoAHQHUBV44H49AnVGHO5NmyZy7z1nzMjYR24XLSJr0WnOndQsXwbp0MzlyHghCEiuYA4VXVfwXU7AVN9758AFovIKMAF1Elh/9LAQWCCiFQG1gCPq2pMCvH1AnoBFC9e/ApCNMZcrrgEN6tXT6Lf/95ARcjn9tC3xJO0b9iNsOAgp8Mzl+mSK6KLiEtE7hWRhSLyN7AN+FNENovIayJSJpAL+rqT2gAzfUUPA/1V9XqgP5DS6FcwUBX4QFWrADHAwJTOr6pjVTVSVSOLFCkSSGjGmFSgqgyc0Iq+O94EoPmpsrxa4z3ubfKAJYtMzp8WxvfAN8AgYJOqegBEpCBwKzBCROao6mQ/r9kCWKuqB3zb3YDHfe9n4h3YTi4aiFbVlb7tWVwgYRhjnDPli/68dvhr3GFCydPK3SV70vW2J5wOy6QSfxJGU1VNSF6oqv8As4HZIhLIyFVn/u2OAu+YRUO8d0s1BnakcK2/ROQPESnrGzxvAmwJ4JrGmDT0zYaZjP35DXaFHcftctHsxI0M7vop+fPmdjo0k4oumTBSShaXsw+AiOQEmgG9kxQ/BLwtIsFAHL7xBxEpCoxX1Za+/R4Fpvi6tHYBPfy5pjEm7Xjcbl6f3oc5cT/jDoPrE10MbTiGSjfVdTo0kwb8HvQWkQEpFB8D1qjqen/OoaqngELJypYD1VLYdz/QMsn2eiDS33iNMWlHPR7emdOVTYd380vYcUJFaBvUicE9XnA6NJOGArlLKtL3+sK33QpYDfQRkZmqOjK1gzPGZDynTicy/svnGX/yVwiDiLgQXmjzNTddU9Dp0EwaCyRhFAKqqupJABF5Ce/gcwO8t7lawjAmC/Mknmbtxs94c/VYNoSd5JoED31LPk2bRp0JCrblUrODQBJGcc59WC4BKKGqsSJyOnXDMsZkNEMmt2SO/E1QqNL01HV0b/oalctUcjosk44CSRifAb+IyDzf9h3AVN+DfHbHkjFZ1HfL32HyxsmszhnLLXFBtL2uG/e17u90WMYBficMVf2PiHwJ1AME6KOqUb6P70uL4IwxzomNi2XE9E7M0524cwoN3dcwpPN0Cue1sYrsKpC7pAS4BcinqsNEpLiI1EhpZlljTOalbjeLf36LMdumsDM0gcoxuXiowTgaVrTup+wukC6p0YAH78N1w4ATeB/cq54GcRljHHAkJoYRU9vyZcgBCIU7KMFLvefZlB4GCCxh1FTVqiKyDkBVj/geojPGZAGrtn7JS8sHEh2qRJzMR+eIB2hRuzviuuSUcyabCCRhJPjWsVAAESmCt8VhjMnETsYc49FptxEVfApCoWtQQ7reN5Kr89uSM+ZcgSSMd/AueHS1iAwHOgD2WKcxmdj4Lwby0cEFnAwWap/IT5MyHejY/PFLH2iypUDukpoiImvwTvwHcKeqbk2bsIwxaWnzjqW89uPTrAk5RU6Bu4MqM7DPJ4QGW/eTubBLJowLzCEF0EJEWqjqG6kckzEmjRw9vo9hn/fga/mToGCl2pGreaztFKqWusbp0Ewm4E8LI4/vZ1m8d0TN923fASxNi6CMMalv8c9TmbB5FJtD4xFVHi/YljvvfokCuezeFeMff6Y3HwogIkvwziV1wrc9hH9XzTPGZFDbt3/F7FWTmOneSGKoUPdYCYZ3GkWhq252OjSTyVzJXFLxQMlUjcYYk2o8HmXL7k30/PlJYlwuQhV65etA+zZPUahALqfDM5lQIAnjU2CViMzBe2vtXcCkNInKGHNF1OPh2Y/v4KuQ38HlosWRsrSr141akXc4HZrJxAK5S2q4iCwC6vuKeqjqOn+PF5GywPQkRaWBwXiXZh0DhAOJQN8LTTfiew4kCtinqq39vbYx2cmSH19n4O4JJIQI1yQoLfLU4YnHPsTlEqdDM5mcP3dJiaoqgKquBdZebJ8L8a3FHeHbPwjYh/e5jnHAUFVdJCIt8a6r0egCp3kc2ArkvVTcxmQ3B/7Zw9DPe7EmaB8JLhcNYkrxXKfPuK6grattUoc/N11/LyKPikjxpIUiEioijUVkEtAtwOs2AXaq6l683VtnEkA+YH9KB4hIMbyr/I0P8FrGZHnvzX6YO+a1ZlnIn5xyuXi04B2833e+JQuTqvzpkrodeADv2helgKNADrzJZgnwpr9reifRCZjqe/8EsFhERvnOWecCx7wFPMO/t/mmSER6Ab0AihcvfrFdjcn0lq37hC9+nc0i2QUu4dYT1Xm151vkCLdGuEl9comepHN3FgkBCgOxqnr0si7onbBwP1BeVQ+IyDvAj6o6W0TuAXqpatNkx7QGWqpqXxFpBDzlzxhGZGSkRkVFXWo3YzKdk7HxvD/3aSbHfwdAbreHMfWnULZkRcJDbGZZc/lEZI2qRqb0WSB3SaGqCcCfVxhPC2Ctqh7wbXfDOzYB3uc6Uupyqgu08Y1xhAN5RWSyqna5wliMyVRUlS27v+bDH0bzfdBOglVpL81oUactlctEOB2eyeICShippDP/dkeBt7XREO/dUo2BHckPUNVBwCCAJC0MSxYmW/n7RBwLlg7lzUMLIAhuOq30q/Yat1a9He/6ZsakrXRNGCKSE2gG9E5S/BDwtogEA3H4xh9EpCgwXlVbpmeMxmREO/5YxTNf9ua38ETKnE6kmqshD7YbxjWFr3I6NJONBLJE64Oq+tGVXExVTwGFkpUtB6qlsO9+4Lxkoao/4G2NGJMtDJ7QlrmyEw0X6p7KQ9faL1InooXTYZlsKJAWxusich/eh+tWAVNVdXPahGWMmfDFo3y3fyXrw2MpGe+hdo72DOozzLqfjGMCSRiHgf8CoXgfwJshIu+o6odpEZgx2dWxk3FMm/887yX8AOFQIy6Md+7/llw58zkdmsnmAkkYx1T1O9/7r0TkbWAlYAnDmFSgbjeLfnqVaVtnsi48kWvjlccjPqRV9bpOh2YMcBmD3iLyLN5nMfIBJ1I9ImOyIbdHGTOzG2NO/wrh0Iqb6dd+HNcXzO90aMacdTl3Sc3GO7VHW+Dl1A3HmOznx6iJfBj1LhtzxFP6dBCPlB9Is7qdnA7LmPMEkjAKiMj1qvob8JuIjAPWAQvTJjRjsjaPx82gT+/gS/4gV5iHxp7iPNN+CtcVKuh0aMakKJCEkRf4QUQOAVuA/IA7LYIyJqub+e3LfLp7FrtDEqh+PB+d675Os6o1nQ7LmIsKJGHcCmwCauJd31ux1oUxAdl/6ACvfN6JH8IOQQjclngVrz68hOBgm//JZHyBLKC0wfd2he9ljPGT253IktWjeXvTOPaFQeSxwjzceDDVyzVEXP6sMmCM85yYS8qYbGXX/v08v7ANm0JPQwg8FN6Q+9q/TqHcYU6HZkxALGEYk4Y+WjKSt/78FEKh2YmCNLzpTtre1t/psIy5LIHMJfUIMEVVj6RhPMZkCb/tXckTXz/I3hDh2gQ3rfM0o0fHUeQJD3E6NGMuWyAtjGuA1SKyFvgYWHypdbyNyW7c7kQGTGrBd0F/QYhQ80gJejb5L7UqRjgdmjFXLJBB7xdE5EXgNqAH8J6IzAA+UtWdaRWgMZnF1IX/4Yt9n7MxLBGAh3LW5LFutgS9yToCXXFPReQv4C+8s9YWAGaJyNeq+kxaBGhMRrdz1w98+OPLLHL9CWEQebwQb/ZYTH4b1DZZTCBjGI/hXU71EN5lVJ9W1QQRceFdJe+iCUNEygLTkxSVBgbjXdtiDN6lVxOBvqq6Ktmx1wOf4O0W8wBjVfVtf2M3Ji14PMqCZZ/x4Y5X+D1ECFbljQovU6Xc7eTPGep0eMakOr8Shngn4K8MtFPVvUk/U1WPiLS+1DlUdTveadERkSBgHzAHGAcMVdVFvjW7RwKNkh2eCDypqmtFJA+wxteq2eJP/Maktu07v2P6Tx8zU3+FEKHpyYo81Kw35W5q6HRoxqQZvxKGryuqSvJkkeTzrQFetwmwU1X3iojinXYEvDPg7k/h/H8Cf/renxCRrcB1eKcoMSbdqCo7dnxNp58HkChCqEfpU7gNd7d/ybqgTJYXyBjGChGprqqrU+G6nYCpvvdPAItFZBTgAupc7EARKQlUwbsWhzHp5njMQQZ81oWVoftBhNrHyvFY2xFUKFHK6dCMSReBziXVR0T2ADGA4G18VArkgiISCrQBBvmKHgb6q+psEbkH+AhoeoFjc+OdXv0JVT1+gX16Ab0AihcvHkhoxlzQazN68ElsFIRCyXilUb66DOg6xpZLNdmK+PsohYiUSKn8Qt1UFzlPW6Cfqt7m2z4G5Pd1ewnelf3ypnBcCLAA7/Mfb/hzrcjISI2KigokPGPOsXLzbD5c8S7rgg/hBhrH1GN4z/fIFWaTJJisSUTWqGpkSp8F8re+2wXKhwUYT2f+7Y4C75hFQ7x3SzXGe8fVOXyJ5CNgq7/JwpgrEZ8QzzMT2vJD6B+4Q4TcbuX1Ki9Tp0pbp0MzxjGBJIyYJO/DgdZAQIPdIpITaAb0TlL8EPC2iAQDcfi6k0SkKDBeVVsCdYH7gY0ist533HOq+mUg1zfGH98ufZnP/7eYpWH/AMI90p6BnZ8gJEd+p0MzxlF+d0mdd6BIGDBfVZunbkipx7qkTCD2HTnBnCX/5cM47/9DKsa5eKPTT1ydN5eNVZhsI7W6pJLLiffhO2MyN1X27VvFEwv7si08nsKJHjoV6sa97XuRJ3dup6MzJsMI5EnvjXhX2QMIAooQ+PiFMRnOiE/aMoXdEA4t466mY4OXqHpLfafDMibDCaSFkfRp7kTggKompnI8xqSblVEfMSZqHGvDT3JjXDAtrurAQ12ft+4nYy4gkNlqA7p91piMKi7BzbuzBjAt7lvc4VAroQCv3DuXgnkLOR2aMRma34sJi8gkEcmfZLuAiHycJlEZk0a+XTaKOz6pyifx31Es3kX/Yq/w4UPLLFkY44dAuqQqqerRMxuqekREqqR+SMakvth4N+Omd2WcZwMEw51SnKe7TCNvrjxOh2ZMphFIwnCJSIEzS7SKSMEAjzcm3bkTE5jz41A+3TWPXaFQ5lQ4vSKHcXv1Fk6HZkymE8gv/NeBn0VkFt67pe4BhqdJVMakAo/bzeOfNuNH12FyB3m4g1vo2X4spa8q6HRoxmRKgayH8T0QhXf6DsG7NoZNL24ypE+/eYVP9szgr5BEGpwoRKda/6V+ZD2nwzImUwtkPYy5qloNW4PCZGD/2xPF8MW9WRseDyHQ3F2E//ZaTHhoiNOhGZPpBdIl9UsqrodhTKqKjTvGhIVPMun4CgiFyodL81yHsZQrfq3ToRmTZaT7ehjGpLYvl0/lky2j2BwWDy4Xj+dtxgPdXsflsgfwjElNgSQMu63EZChHDm7l27XTGbF/FqfDhDrHr6V3/b5UjbjT6dCMyZICSRi/A/cBpVV1mIgUB64B7Alwk+727d9Mx6/u4ViQi9yq9CvShbva97d1tY1JQ4EkjNGAB+9dUsOAE3iXS62eBnEZkyL1eBg6+R5m63YIctHoRAT31nuY2hEXXQreGJMKAkkYNVW1qoisg7NPeoemUVzGnOfL70fy2q5JHAp2UcDtoXWuqjzT7VOnwzIm2wgkYSSISBC+Kc5FpAjeFodfRKQsMD1JUWlgMN6lWcfgXcUvEeirqqtSOP524G28U6uPV9URAcRuMrGjR3bz3Oe9WBb8FwS7qHW8CP3bz6JcUXsAz5j0FEjCeAeYA1wtIsOBDsCL/h6sqtuBCABf4tnnO984YKiqLhKRlsBIoFHSY337v493eddoYLWIzLcHB7O+sfP+w+cHp7MvxHvH07NXd6RLtxccjsqY7CmQ6c2niMgaoImvqK2qbrvM6zYBdqrqXhFRIK+vPB+wP4X9awC/qeouABGZBrTFHiLMsr5f8xFT101kRchRCBGqHbuR17u+Q6H81zsdmjHZ1iUThojMT17k+9lcRFDVNpdx3U7AVN/7J4DFIjIK73TrKY1eXgf8kWQ7Gqh5gXh7Ab0AihcvfhmhGScluD0s3biSgRveJC5EuCbBw9OVXqFy2WYUypfD6fCMydb8aWHUxvvLeiqwkn8TxmXxDZS3AQb5ih4G+qvqbBG5B/gIaJr8sBROpSmUoapjgbEAkZGRKe5jMh6PO5HvV77J/K3f8l3wPnAJdya2odft3bm+WBmnwzPG4F/CuAbv2EFn4F5gITBVVTdf5jVbAGtV9YBvuxvwuO/9TGB8CsdEA0n7IoqRcteVyaQ++WoQrx/6CoKhUmwQda5uQZ+2/yXIntY2JsO4ZMJQVTfwFfCViIThTRw/iMgwVX33Mq7ZmX+7o8D7i78h3rulGgM7UjhmNVBGRErhHSzvhDd5mUxu2+7veO6bZ9gVEkcQcLunA092fpoi+XI5HZoxJhl/pzcPA1rh/WVfEu8dU58HejERyYm3tdI7SfFDwNsiEgzE4Rt/EJGieG+fbamqiSLyCLAY7221H19BC8dkAKrK8A9vY3qOvyAUysUF0aP8Y9xe70GnQzPGXICoXrybX0QmARWARcA0Vd2UHoGlhsjISI2KinI6DJPMzMVP88Pe5SwNOwnAPdqKF7q9gnfZFWOMk0RkjapGpvSZPy2M+/HOTnsT8FiSf9RnZqvNe6EDjUnqyMkYPpjdn6msgDC4IV6Y0OlbCuQp4nRoxhg/+DOG4UqPQEwWpsrKdR8zOmo0a8PiKZTo4ZEyI2lX73ZcriCnozPG+CmQJ72NCVhcgpuPpz/AB+61EAbNE0rR/85xXFf4aqdDM8YEyBKGSTN79/7CoK8eYVNYHNfHh3B/8YfoeNvDtrCRMZmUJQyT6hITExg1pStfuDcQHyo0kZI832kihfMWdjo0Y8wVsIRhUtXMpSMZvutT3AI3ng6nTemn6dG8k9NhGWNSgSUMkypi407z2sy+zPSsAoF2lGZI7zmIy+6ZMCarsIRhrkhCfCxjF/Ri7pEo/goOpnxMME80+oRa5So6HZoxJpVZwjCXLSEhkf5TWvCj6zBhLhftQ2pzf7sR3HCNLWxkTFZkCcMEzONxM+3H4byzewYxQULjE1fxQMMhVK5Y3+nQjDFpyBKGCcj+/ZvovbATe0KFQqo0k9IM6jmbnOEhTodmjEljljCMX+Jij/DyjK4sS9zJ0RAXtQ5XoWerodQsW8rp0Iwx6cQShrmkRUsn8Om2t9kY5obgIO4Lq8LAAZ84HZYxJp1ZwjAXtG/PMqb99D4TEzdDGFQ7UYS3us8jf+48TodmjHGAJQyTouXrv2HYmsf5M9iFS5UhNwygec2u5Ay1vzLGZFf2r9+cIy72CB8v/g8fHPsagl3UO16JHg26U6NyM6dDM8Y4LN0ShoiUBaYnKSoNDAZqA2V9ZfmBo6oakcLx/YGegAIbgR6qGpeGIWc7v6yfSb91Q4l3CeEeDz0K1KV3lzEEBdnT2saYdEwYqrodiAAQkSC8a3PPUdW3zuwjIq8Dx5IfKyLXAY8B5VQ1VkRm4F3Xe2KaB54NnE44Re9JLVgT8g+4hIYnStOh4X9pVNGe1jbG/MupLqkmwE5V3XumQLxL+d0DNL7AMcFADhFJAHIC+9M8ymxg7Pz/MOXQdP4JEUrEK3dd04IHu73mdFjGmAzIqYTRCZiarKw+cEBVdyTfWVX3icgo4HcgFliiqktSOrGI9AJ6ARQvXjxVg85Kvosaz+vr3+P3EDcECdWPXsuIB+ZzVZ5wp0MzxmRQ6d45LSKhQBtgZrKPOnN+EjlzTAGgLVAKKArkEpEuKe2rqmNVNVJVI4sUsbWiUzJuziu8uOFNfg9xk9/tYVzlEYzuu8iShTHmopxoYbQA1qrqgTMFIhIMtAOqXeCYpsBuVT3o2/9zoA4wOY1jzTLU42HSgj789Pdmfgk5DkEu2iQ0Yeh9QwnOkc/p8IwxmYATCSOllkRTYJuqRl/gmN+BWiKSE2+XVBMgKu1CzFpiTicy74e3eP3ICgiBirEuekS8SOMa7Qmy5VKNMX5K14Th+4XfDOid7KPzxjREpCgwXlVbqupKEZkFrAUSgXXA2HQIOVNTj4e/9q/imYWPsT48lvxuD10KPMLdd95LwXzWqjDGBEZU1ekY0kxkZKRGRWXfhsjrn93LxISNANSNzctdlXrRvE43h6MyxmRkIrJGVSNT+sye9M6Cfln1Ph9u+JRNQSfIr0rr0NY8+dAIgu0BPGPMFbCEkYUkuj2MmtKPOe6lnApxccvpEF6+fRI3Fq/sdGjGmCzAEkYW8d3Kd5i4cQrrQk5RIkFoUewl+rZoj/d5SGOMuXKWMDK5Y6fiGTXzQeayHkKgzul8vNVtMTnCcjkdmjEmi7GEkYlt27mUl755gi3hCVSNCeWucs/QtsHdiMvGKowxqc8SRiZ06vRpRszowhfurQSHKu1dETzS5UMK583tdGjGmCzMEkYmM+nboYyKngVA1Zhwmt/wJPfe3tnhqIwx2YEljEzi7yN/MXJOHxYH7QTg3qCbeLbPTFx2q6wxJp1YwsjgEhNieXVmJxbG/caJIBcRJ8Lp13gytcqVvfTBxhiTiixhZGB/HDzEO1905KugvyHIRc+c9eh45xtcky+H06EZY7IhSxgZUHzcMeb/+Aav7Z/FqSAXNU8WYVib0RS9/manQzPGZGOWMDKYkyf+ptP0xuwNEYp4lC5569Ht7lHkzWnPVRhjnGUJI6PwuPnvp3ey0L2LkyEuGsRUplO9QdSvVN7pyIwxBrCEkSEsi5rOq+uHszdECXEJraQ4I/ra2lDGmIzFEoaDovdF8daS5/hW9pMYIpQ/UYAX75pL+esKOh2aMcacxxKGQ6Z+/QGT9r7PvhABhH55m9Kzy+s2BbkxJsNKt4QhImWB6UmKSgODgdrAmYcK8gNHVTUihePzA+OBCoACD6jqirSLOG1s37mE77csZvThxWiIUO/EzTzZuj83lqzjdGjGGHNR6ZYwVHU7EAEgIkHAPmCOqr51Zh8ReR04doFTvA18paodRCQUyJmmAaey2PhEZq/4ktE7B3EiyEUuVZ4v3Y/6kQ+SP2eo0+EZY8wlOdUl1QTYqap7zxSId+GGe4DGyXcWkbxAA6A7gKrGA/HpEmkqiDv1D09Mv5efXfsgyEXLmAbcFtGCJnVaOx2aMcb4zamE0QmYmqysPnBAVXeksH9p4CAwQUQqA2uAx1U1JvmOItIL6AVQvHjxVA06UKrKjG8+4L/7PwAXlIsT2hVvTccWLzsalzHGXA5R1fS9oLc7aT9QXlUPJCn/APhNVV9P4ZhI4BegrqquFJG3geOq+uLFrhUZGalRUVGpWwE//bplDk+ueJEDwd4V7+qcLM1THSZR5ur8jsRjjDH+EJE1qhqZ0mdOtDBaAGuTJYtgoB1Q7QLHRAPRqrrStz0LGJimUV4mVWX4J71Z4v6JI8Euyp4WHir3GM3r9XQ6NGOMuSJOJIzOnN8d1RTYpqrRKR2gqn+JyB8iUtY3eN4E2JLGcQZs2sJHWbzvZ6LC4iHIxW2x1RnV+yNbV9sYkyWka8IQkZxAM6B3so/OG9MQkaLAeFVt6St6FJji69LaBfRI43D9luj2MPGLUbx97AcIgxtOKx93/JqC+a51OjRjjEk16ZowVPUUUCiF8u4plO0HWibZXg+k2K/mpK0bP+O9X95laehJwjxKr8ID6HFbZ0LCbApyY0zWYk96XyZVZcrCZ3n18CIIhXqxV/Fw05epdGNNp0Mzxpg0YQnjMuz83xIG/vA8u0NiKeJR2hV8iF73PUZosE3rYYzJuixhBGjkxC5MZT2JYUKN+Dz0azCSqmXrOx2WMcakOUsYfvph5fsM2/whB4OUkrEhtL9xMN2b3uV0WMYYk24sYVzCPydjGf7ZfSwJ2wFBUCcuN693XULunHmcDs0YY9KVJYyL2LD+U/6z6nW2hbmJOJmDTlWG07JWE8RlYxXGmOzHEkYKjpyM48NZXZnKFgiFrlKJ3j0mkDdnmNOhGWOMYyxhJHMkJp7206txMNhFhdhwOtzQn/bN7nM6LGOMcZwljGTyh0EluYbiOa7j0XvHExJia1UYYwxYwjiPBIfy1gPfOh2GMcZkODZ6a4wxxi+WMIwxxvjFEoYxxhi/WMIwxhjjF0sYxhhj/GIJwxhjjF8sYRhjjPGLJQxjjDF+EVV1OoY0IyIHgb2XeXhh4FAqhpMZWJ2zvuxWX7A6B6qEqhZJ6YMsnTCuhIhEqWqGW0M8LVmds77sVl+wOqcm65IyxhjjF0sYxhhj/GIJ48LGOh2AA6zOWV92qy9YnVONjWEYY4zxi7UwjDHG+MUShjHGGL9YwkhGRG4Xke0i8puIDHQ6ntQiIteLyPcislVENovI477ygiLytYjs8P0skOSYQb7vYbuINHcu+ssnIkEisk5EFvi2s3R9AUQkv4jMEpFtvj/v2lm53iLS3/d3epOITBWR8KxYXxH5WET+FpFNScoCrqeIVBORjb7P3hER8TsIVbWX7wUEATuB0kAo8CtQzum4Uqlu1wJVfe/zAP8DygEjgYG+8oHAq7735Xz1DwNK+b6XIKfrcRn1HgB8BizwbWfp+vrqMgno6XsfCuTPqvUGrgN2Azl82zOA7lmxvkADoCqwKUlZwPUEVgG1AQEWAS38jcFaGOeqAfymqrtUNR6YBrR1OKZUoap/qupa3/sTwFa8/9ja4v0Fg+/nnb73bYFpqnpaVXcDv+H9fjINESkGtALGJynOsvUFEJG8eH+xfASgqvGqepSsXe9gIIeIBAM5gf1kwfqq6lLgn2TFAdVTRK4F8qrqCvVmj0+SHHNJljDOdR3wR5LtaF9ZliIiJYEqwErgalX9E7xJBbjKt1tW+C7eAp4BPEnKsnJ9wds6PghM8HXFjReRXGTReqvqPmAU8DvwJ3BMVZeQReubgkDreZ3vffJyv1jCOFdKfXlZ6r5jEckNzAaeUNXjF9s1hbJM812ISGvgb1Vd4+8hKZRlmvomEYy32+IDVa0CxODtqriQTF1vX599W7zdLkWBXCLS5WKHpFCWaeobgAvV84rqbwnjXNHA9Um2i+Ft3mYJIhKCN1lMUdXPfcUHfM1UfD//9pVn9u+iLtBGRPbg7VpsLCKTybr1PSMaiFbVlb7tWXgTSFatd1Ngt6oeVNUE4HOgDlm3vskFWs9o3/vk5X6xhHGu1UAZESklIqFAJ2C+wzGlCt+dEB8BW1X1jSQfzQe6+d53A+YlKe8kImEiUgoog3ewLFNQ1UGqWkxVS+L9c/xOVbuQRet7hqr+BfwhImV9RU2ALWTdev8O1BKRnL6/403wjs9l1fomF1A9fd1WJ0Sklu/76prkmEtzeuQ/o72AlnjvINoJPO90PKlYr3p4m54bgPW+V0ugEPAtsMP3s2CSY573fQ/bCeBOioz2Ahrx711S2aG+EUCU7896LlAgK9cbGApsAzYBn+K9MyjL1ReYinecJgFvS+HBy6knEOn7rnYC7+Gb8cOfl00NYowxxi/WJWWMMcYvljCMMcb4xRKGMcYYv1jCMMYY4xdLGMYYY/xiCcOYSxCRQiKy3vf6S0T2+d6fFJHRaXTNJ0Sk6yX2mSYiZdLi+sakxG6rNSYAIjIEOKmqo9LwGsHAWryzCydeZL+GQBdVfSitYjEmKWthGHOZRKRRknU2hojIJBFZIiJ7RKSdiIz0rTvwlW9aljNrEfwoImtEZPGZaR2SaQysVdVEEblBRNYmuWYZETkzP9YyoKkvwRiT5ixhGJN6bsA7nXpbYDLwvapWBGKBVr6k8S7QQVWrAR8Dw1M4T11gDYCq7gSOiUiE77MewETfZx6801ZXTqP6GHMO+5+JMalnkaomiMhGvItxfeUr3wiUBMoCFYCvfYucBeGd6iG5a/HOh3TGeKCHiAwAOnLu+g1/452l1d9ZeY25bJYwjEk9p8H7P38RSdB/Bwg9eP+tCbBZVWtf4jyxQHiS7dnAS8B3wBpVPZzks3Df/sakOeuSMib9bAeKiEht8E43LyLlU9hvK3DjmQ1VjQMWAx8AE5LtexOwOW3CNeZcljCMSSfqXfa3A/CqiPyKd8bgOinsugjvMqtJTcE72/CSMwUicjUQq74V14xJa3ZbrTEZkIjMAZ5R1R2+7aeAfKr6YpJ9+gPHVfUjh8I02YyNYRiTMQ3EO/i9w5c8bsB7u21SR/Gu/2BMurAWhjHGGL/YGIYxxhi/WMIwxhjjF0sYxhhj/GIJwxhjjF8sYRhjjPHL/wGwFteyg9Ky4gAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEGCAYAAAB2EqL0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA2zUlEQVR4nO3dd3hUZfbA8e9JoVchSjNAkB5IgNB7VSJNRCkWEBVRVhF+1lUsuO6qyyogrghIWVFEQFApglgAESkBpEekBwQp0gIh7fz+mAFDCDADmdxkcj7Pk4e579w797wj5vCW+76iqhhjjDFXE+B0AMYYY3IGSxjGGGM8YgnDGGOMRyxhGGOM8YglDGOMMR4JcjoAXypZsqRWqFDB6TCMMSbHiImJOaKqIRm959cJo0KFCqxZs8bpMIwxJscQkT2Xe8+6pIwxxnjEEoYxxhiPWMIwxhjjEb8ew8hIUlIScXFxJCQkOB1KrpMvXz7KlStHcHCw06EYY65BrksYcXFxFC5cmAoVKiAiToeTa6gqR48eJS4ujooVKzodjjHmGuS6LqmEhARKlChhySKLiQglSpSwlp0xOViuSxiAJQuH2PduTM6WKxOGMcb4q19XLuCXT1/FF1tXWMJwQKFChdi9ezf58+cnMjKSGjVqMHDgQFJTU9m9ezfh4eFXvP6VV15hxIgRF5VVqFCBI0eOeBVHdHQ0x48f9zZ8Y0x2lHyOAxPvpcqCXpSM/ZizZ05l+i0sYTioUqVKrF+/ng0bNrBlyxbmzJmTJfdVVVJTU5k/fz7FihXLknsaY3wkNZW9303k4IjGlNn7FR8G9OBY3yUUKFgk029lCSMbCAoKokmTJvz222+Z8nlvv/024eHhhIeHM3LkSAB2795N9erVeeyxx6hbty779u270CoZO3YskZGRREZGUrFiRVq3bg3AtGnTqFWrFuHh4Tz77LMXPr9QoUK88MILRERE0KhRIw4dOgTAjBkzCA8PJyIighYtWmRKXYwxlxd/ZC+xH/YndOkQSiXs4Mfi3ejz3DhqVSjtk/vlumm1ab361Wa2HDiZqZ9Zo0wRXu5c06trzpw5w7fffsvw4cM9vuadd95h6tSpF44PHDgAQExMDJMmTWLlypWoKg0bNqRly5YUL16c2NhYJk2axH//+9+LPmvgwIEMHDiQpKQk2rRpw9ChQzlw4ADPPvssMTExFC9enA4dOjBnzhy6detGfHw8jRo14vXXX+eZZ55h/PjxvPjiiwwfPpyFCxdStmxZ6+oyxpfOnebw5u9h7hCqph5mmnag1f3DaFYxHAJ81w6wFoaDduzYQWRkJE2bNuX222+nY8eOHl87ZMgQ1q9ff+GnTJkyAPz444/ccccdFCxYkEKFCtG9e3eWLVsGQPny5WnUqNFlP3Pw4MG0adOGzp07s3r1alq1akVISAhBQUHcc889LF26FIA8efLQqVMnAOrVq8fu3bsBaNq0Kf369WP8+PGkpKRcy1dijLmKuAMHODS6LSFf3ktI6mHmR7xLkycmU7pSbZ8mC8jlLQxvWwKZ7fwYRma60syIggULXva9yZMns2fPHsaMGXPVzwkODr4wRTYwMJDk5GQAxo4dy8qVK5k3bx6RkZGsX7+eEiVKXEs1jDEZ2L5yHpUX9AHg4zw9aNjubqIb3Jpl97cWhp9p0aIFc+bM4cyZM8THxzN79myaN29+xWtiYmIYMWIEU6dOJcD9L5SGDRuyZMkSjhw5QkpKCtOmTaNly5ZX/JwdO3bQsGFDhg8fTsmSJdm3b1+m1cuY3EyP7+XoG7WpvKAPKQjzSg+i+YDR3JKFyQKysIUhIlWB6WmKwoCXgO+BsUAhYDdwj6peMrAgIrcBo4BAYIKqvuHrmH0hOTmZvHnzXvGc2NhYypUrd+H4nXfe4a677vLo8+vWrUu/fv1o0KABAA899BB16tS50G2UkTFjxnDs2LELg91RUVFMmDCBf/3rX7Ru3RpVJTo6mq5du17x3k8//TTbt29HVWnbti0REREexWyMuQxVFn/wFO0OTqAEMCU1moiew7i9Zg1HwhFfPNxx1ZuKBAL7gYbATOApVV0iIv2Biqo6LIPzfwXaA3HAaqC3qm650n2ioqI0/QZKW7dupXr16plWF2/98ssvPPzww6xatcqxGJzk9PdvTE6xf+Vsgr57hZvO7Qbg57L9aPDgSAICfLtigojEqGpURu85NYbRFtihqnvcLY+l7vJvgIXAsHTnNwB+U9WdACLyKdAVuGLCyG7Gjh3L6NGjL0x1NcaY9BJ+38baz/9Dk8OfATBH2tJsyFQaFSngcGTOJYxewDT3601AF+AL4C7g5gzOLwuk7RCPw9U6uYSIDAAGAISGhmZSuJnj/PRVY4xJLzVVWbf2ZyosfIAmSb8DsC16Bq1rtaZo/uyxJUCWD3qLSB5cCWKGu6g/MEhEYoDCQGJGl2VQlmFfmqqOU9UoVY0KCclwH3NjjMlWkncu46cP/496c2+jRNLvfFb0AeLvW0C1Bh2yTbIAZ1oYHYG1qnoIQFW3AR0ARKQKcHsG18RxccujHHDAx3EaY4xPpaYqSQc2kPd/nWgG/B5Qmt+iXqJN8+4ULJzP6fAu4cS02t781R2FiNzo/jMAeBHXjKn0VgOVRaSiu4XSC/gyC2I1xpjMpwqnD/PTew+j49sCMKbQEwQ9tpTm0X0omQ2TBWRxC0NECuCa6fRImuLeIjLI/fpzYJL73DK4ps9Gq2qyiPwN14B4IDBRVTdnYejGGJNpTi79L0W+/zvNgG35arG70j0MvPMRggKz96NxWRqdqp5R1RKqeiJN2ShVreL+eU7d83xV9YCqRqc5b777nEqq+npWxp3ZXn/9dWrWrEnt2rWJjIxk5cqVHl330ksvsXjxYgCWLVtGzZo1iYyMZMWKFcyfP/+aYjl+/Pgla0uldejQIfr06UNYWBj16tWjcePGzJ49G4AffviBokWLUqdOHapVq8ZTTz11TTEYk1skxa1l6bsDkO9fA2Bq2Rcp/39LuO3uR7N9sgB70jvLrVixgrlz57J27Vo2bNjA4sWLufnmjCaGXWr48OG0a9cOgI8//pinnnqK9evXExsb65OEoap069aNFi1asHPnTmJiYvj000+Ji4u7cE7z5s1Zt24d69atY+7cuSxfvvya4jDGn6kqK+e8R8rEaFocnc5JKcLGbt9w78NPkz9PoNPhecwSRhb7/fffKVmy5IWnvUuWLEmZMmVYtWoV3bt3B+CLL74gf/78JCYmkpCQQFhYGAD9+vVj5syZTJgwgc8++4zhw4fTu3dvXnrpJaZPn05kZCTTp08nPj6e/v37U79+ferUqcMXX3wBwObNm2nQoAGRkZHUrl2b7du389xzz11YBPHpp5++KNbvvvuOPHnyXDQVuHz58jz++OOX1Ov8ZlD79+/3yfdmTE51atdqVr0/gIbr/05SCkysOpayL/9KrcgGTofmtVy9+CALnoODGzP3M0vVgo6XX7WkQ4cODB8+nCpVqtCuXTt69uxJy5YtqVu3LuvWrQNc3U3h4eGsXr2a5ORkGja8+JGThx56iB9//JFOnTrRo0cPJk+ezJo1ay4sHPj3v/+dNm3aMHHiRI4fP06DBg1o164dY8eOZfDgwdxzzz0kJiaSkpLCG2+8waZNmzJcBHHz5s3UrVvXo2r/+eefbN++3fbBMMbtz/hEfln2Fc1+fpiGpLAxTwSlHp1L/+KZv7FRVrEWRhYrVKgQMTExjBs3jpCQEHr27MnkyZMJCgrilltuYevWraxatYqhQ4eydOlSli1bdtXFA9NbtGgRb7zxBpGRkbRq1YqEhAT27t1L48aN+ec//8mbb77Jnj17yJ8/v1efO2jQICIiIqhfv/6FsmXLllG7dm1KlSpFp06dKFWqlFefaYw/Sok/xtrRPWn1c3/+0GKsaDaJGs98R0gOThaQ21sYV2gJ+FJgYCCtWrWiVatW1KpViylTptCvXz+aN2/OggULCA4Opl27dvTr14+UlJRL9u++GlVl1qxZVK1a9aLy6tWr07BhQ+bNm8ett97KhAkTLnR3ZaRmzZrMmjXrwvF7773HkSNHiIr6a5mZ5s2bM3fuXH799VeaNWvGHXfcQWRkpFfxGuNPdi4YTdjKYbQF5hS6m2rd/07jsIpOh5UprIWRxWJjY9m+ffuF4/Xr11O+fHnAtTT5yJEjady4MSEhIRw9epRt27ZRs+aV9+0oXLgwp079teH7rbfeyrvvvnthT4vzXV07d+4kLCyMJ554gi5durBhw4ZLrk2rTZs2JCQk8P77718oO3PmTIbnVqlSheeff54333zTg2/BGP/zx6+rWPVaa8JWDmNjQFVWNJ1A1/8bRzU/SRZgCSPLnT59mr59+1KjRg1q167Nli1beOWVVwDXHhSHDh26MA5Qu3ZtateufWGzostp3bo1W7ZsuTDoPWzYMJKSkqhduzbh4eEMG+Zay3H69OmEh4cTGRnJtm3buP/++ylRogRNmzYlPDz8kkFvEWHOnDksWbKEihUr0qBBA/r27XvZpDBw4ECWLl3Krl27rvNbMibnSDyXwC+fvkrhTzrRIGUtscVaEDrkOxq3v+uq/+/mNI4sb55VsuPy5rmdff/Gb6hy8qeJFPlmKABLNZJCd75L3dq1HQ7s+mTH5c2NMSbHOpecwpGv36TsGldre0XYYKp2fZ6bino3kSSnsYRhjDGeSknm8Iav+fPLF6miu/iWBlTq/jKNazdzOrIskSsThqr6Xd9iTuDP3Z8mF0hO5I/3buXGP9dSQPOxosz93HDrS1SokHu2Uch1CSNfvnwcPXqUEiVKWNLIQqrK0aNHyZcve67CacyV/PrNREr89Bo36jEmFHiIqG6DaFzl8lPS/VWuSxjlypUjLi6Ow4cPOx1KrpMvXz7KlSvndBjGeCzh91h2T3uSaid/AmBZkU7cMeCflCiU1+HInJHrEkZwcDAVK/rPvGhjTObThBNsmPs+FTePISz1DJM0mmYDx9C8TAmnQ3NUrksYxhhzJZvWLKPAN88QcW6L67jOS/S+/UnyBeecVWV9xRKGMcYA547sZsv6n6m57DHySAqfBXbijv7PEF6mNth4J5CFCUNEqgLT0xSFAS8BP+DaljUfkAw8pqqrMrh+CPAQoMBG4AFVTfBx2MYYP5eaquzYd4BS/2tJnZTjpCIsifg3TVr2IfiGQk6Hl61kWcJQ1VggEkBEAoH9wGxgPPCqqi4QkWjgLaBV2mtFpCzwBFBDVc+KyGe49vWenFXxG2P8UGoqSya/TOu9owGYlPdeWnbsSctIW6Y/I051SbUFdqjqHhFR4Pyav0WBA5e5JgjILyJJQIErnGeMMVf157ovKPzFA7QmhZMUZEfoXfS+720bq7gCpxJGL2Ca+/WTwEIRGYFrMcQm6U9W1f3u9/cCZ4FFqrooow8WkQHAAIDQ0NDMj9wYk6Ol/rmXjRMGEhHv2k74s4J9aHD/P6lzU3GHI8v+sny1WhHJA3QBZriLHgWGqOrNwBDgwwyuKQ50BSoCZYCCInJvRp+vquNUNUpVo0JCcs8TmMaYq9u/dAqHP+hCRPxyTlKI2PaTufvp96lgycIjTrQwOgJrVfWQ+7gvMNj9egYwIYNr2gG7VPUwgIh8jqslMtXHsRpj/MCp1dPY+eMMIk58C8AHefrS/9mRVA20HR684UTC6M1f3VHgGotoiWu2VBtgewbX7AUaiUgBXF1SbYE1GZxnjDEXpKQqm9YuJ2LeQCKA3VKOP3t+SZ8KoQRbsvBaliYM9y/89sAjaYofBkaJSBCQgHv8QUTKABNUNVpVV4rITGAtrqm364BxWRm7MSYHUSUpdiHrFk2lwbGvSNYAplR9nz53dKNCflvP7Frlug2UjDH+LTVVObX6Y4ouGATA90HNSWz0BC1btrUZUB6wDZSMMblC0rG9LJ7wAk3iF3OcgnwbOZLo2+8kfx5LFJnBEoYxxi9snT6M6ltH0xFYVaAZBW97mTsjGjgdll+xhGGMydFOblrIvvn/oeaZlezjJlbVeZM7u97hdFh+yRKGMSZHij+XzM/zptByw9PUJIXYAnUpP2gOdxYs6nRofssShjEmxzm3fQm/zf4Xbc+sYLuW42TnidSr18BWlfUxSxjGmBwjJVX59fPXqL7pP0QAP5a8i7Deb1O5RJGrXmuunyUMY0yOcHrPWuI+epTqydtYGVCHgm2foVnTaKfDylUsYRhjsrXUsyfY8enTVNwzk/IayJqy91Cl9wiKFy7gdGi5jiUMY0y2dfiHDyj2w9+pqKnMSG1Dqc7DaN0g0umwci1LGMaYbOfkkd85uOBNquyYBMC6ms9yd4/nCQywQW0nWcIwxmQf506zf+lk5KfRVNFDzElpQs1HP6JOmZJOR2awhGGMySbOnkvm2IQ+lD28hNOan5+rPkVYs8eoXKaE06EZN0sYxhhnqbJl8RRK/fgiZeUUH+Z/gA79X6VRiD2Al91YwjDGOCbleBwn32tLjaSDHJYbWFLqAW7r+Tplixd0OjSTAUsYxpisl5LMz+OfoMrBr7iBk4ymFy3uf4WWYaWdjsxcgSUMY0yW2r9qDoW/foJGqScAiK38ME/cM8LhqIwnsixhiEhVYHqaojDgJVxbs44F8uHaTe8xVV2VwfXFcO33HQ4o0F9VV/g2amNMZknYu47YReOJiPsYgFkBHWj71MdULZDH4ciMp7IsYahqLBAJICKBwH5gNjAeeFVVF4hINPAW0CqDjxgFfK2qPUQkD2CPeRqTA6gq635aRI3F9xOhCQBs7TiDduGtKGrJIkdxqkuqLbBDVfeIiALnVw4rChxIf7KIFAFaAP0AVDURSMyaUI0x1yp5909sX/U1dbeMAmDqjU/RtWM01SvWczgycy28ThgiUhBIUNWU67hvL2Ca+/WTwEIRGQEEAE0yOD8MOAxMEpEIIAYYrKrxGcQ3ABgAEBoaeh0hGmOu1bnkFE7sXMONn3SkOnCUouxpNJwebe+zfbVzsICrnSAiASLSR0TmicgfwDbgdxHZLCL/FpHK3tzQ3Z3UBZjhLnoUGKKqNwNDgA8zuCwIqAu8r6p1gHjguYw+X1XHqWqUqkaFhIR4E5oxJjOcO8XWd3tw4ycdSFHhvRue5Vj/n6l7Wz9LFjmcJy2M74HFwPPAJlVNBRCRG4DWwBsiMltVp3p4z47AWlU95D7uCwx2v56Ba2A7vTggTlVXuo9ncpmEYYxxzoGVsyizoD+RwA/UpWirxxnUurvTYZlM4knCaKeqSekLVfUYMAuYJSLBXtyzN391R4FrzKIlrtlSbYDtGdzroIjsE5Gq7sHztsAWL+5pjPGhUxvns/fLf1A+cQcITCnQl86DRnBDQRvU9idXTRgZJYtrOQdARAoA7YFH0hQ/DIwSkSAgAff4g4iUASao6vkdUh4HPnZ3ae0EHvDknsYY31FVls37mIh1L3Jjciqx+SIo3WcMfctXcTo04wMeD3qLyNAMik8AMaq63pPPUNUzQIl0ZT8Cl0yZUNUDQHSa4/VAlKfxGmN86+iPkzjx0yRanPmFJA1kRtgIHurb3+mwjA95M0sqyv3zlfv4dmA1MFBEZqjqW5kdnDEm+0lISmH9z9/S6NsnKQEsy9uCkN4f0D/UlvXwd94kjBJAXVU9DSAiL+MafG6Ba5qrJQxj/Fzq4d+ImfJ3Gp76hj8oxrZWH9CsRQcCAq864dL4AW8SRigXPyyXBJRX1bMici5zwzLGZDc7Z79G2C8jaKgBLCrYmWrdnqFFlXCnwzJZyJuE8Qnws4h84T7uDExzP8hnM5aM8VPHNyzg2JcvEJa8g40B1TjSeBgd292OiG2Xmtt4nDBU9TURmQ80AwQYqKpr3G/f44vgjDHOSUpOJnbmcKptG0MxUogt3orQvpOoVewGp0MzDvFmlpQA1YGiqjpcREJFpEFGK8saY3IwVeLXfEzSgmGEpx5jgTahRM/RNKjh1aIOxg950yX1XyAV18N1w4FTuB7cq++DuIwxDjidkMSOT4YSsfd/AKyr/ASterxM/ry2dY7xLmE0VNW6IrIOQFX/dD9EZ4zxA6d3reL0R/cSkXqIWdKeWj1eoE7NOk6HZbIRbxJGknsfCwUQkRBcLQ5jTA6WnJzMr5Meocb+mSRpIRZVHkaTTo9RuphtOWMu5k3CGI1rw6ObROR1oAfwok+iMsZkib0L3uaGlf+mBmeYFtCJmzq/TIc6tqyHyZg3s6Q+FpEYXAv/AXRT1a2+CcsY40unDu9l5+x/EHFgOscpTEylQXTvNZy8wTZWYS7vqn87LrOGFEBHEemoqm9nckzGGB/RhJNsnT6MGrsmU1MD+FKbcnP/KdSrYHvHmKvz5J8Thd1/VsU1I+pL93FnYKkvgjLGZL7Y2C3o/GeocWIZAJubvEOTpr0pWSivw5GZnMKT5c1fBRCRRbjWkjrlPn6Fv3bNM8ZkU+dOHGTj4k8I3/BP8kkS04K7cdeTI4koWNzp0EwOcz1rSSUCFTI1GmNMpklNVX4/9if5x7YgKvkQSQSyqvrztGg/iKCCha/+Acak403C+AhYJSKzcU2tvQOY4pOojDHXJzWVmPcfov7hWQD8I99Q7ryrHw0qlXc4MJOTeTNL6nURWQA0dxc9oKrrPL1eRKoC09MUhQEv4dqadSyQD0gGHrvcciPu50DWAPtVtZOn9zYmN/lz40KKz7qb+sCvQVUIqNCE5/u8RGCALRZoro8ns6REVRVAVdcCa690zuW49+KOdJ8fCOzH9VzHeOBVVV0gItG49tVodZmPGQxsBYpcLW5jcpvk/b9w4JNBlD69BQSmFriPVg+9SbkbCjodmvETnux68r2IPC4ioWkLRSSPiLQRkSlAXy/v2xbYoap7cHVvnU8ARYEDGV0gIuVw7fI3wct7GeP3fp0/BhnfitD4jQRKKpubjebeZ8ZYsjCZypMuqduA/rj2vqgIHAfy40o2i4B3PN3TO41ewDT36yeBhSIywv2ZTS5zzUjgGf6a5pshERkADAAIDQ290qnG5HgJG78gdtUiIvZNBWBMocEMGvISNQPtATyT+eQqPUkXnywSDJQEzqrq8Wu6oWvBwgNATVU9JCKjgSWqOktE7gYGqGq7dNd0AqJV9TERaQU85ckYRlRUlK5Zs+ZqpxmT45w5l8iqeZNpteFpAA5Skvi+iyhbrgL5ggMdjs7kZCISo6pRGb3n1T9DVDUJ+P064+kIrFXVQ+7jvrjGJsD1XEdGXU5NgS7uMY58QBERmaqq915nLMbkPKcPs2PiY7Q6togzmpdPqo+hY+vWVLrJntY2vuVEu7U3f3VHgau10RLXbKk2wPb0F6jq88DzAGlaGJYsTK5y5PQ5Tm9aSIWv76MWsCxfa0rc/iIPhtez7VJNlsjShCEiBYD2wCNpih8GRolIEJCAe/xBRMoAE1Q1OitjNCY7Sjx+kP2j7yAidQvHtDBflxtMu55/48Yi+Z0OzeQiHo9hiMiDqvqhj+PJVDaGYfzBdzPeo/7WN8ibEs/q/M0o1+UFytdo6HRYxk9l1hjGf0TkHlwP160Cpqnq5swI0BhzqaM/TubUkndpk/QbB/QG5lZ+k173PGTdT8Yx3iSMo8A/gDy4HsD7TERGq+oHvgjMmNzqxNkkfv5mBreuHUwJYHv+2pR/fB69C9jzqsZZ3iSME6r6nfv11yIyClgJWMIwJjOokrBuOge+fpdbEzdxSIuxvdNsmkXVAWtVmGzA60FvEXkW17MYRYFTmR6RMblQaqoSN/efhK59i+rA8vKPUqnTUzQLKel0aMZccC2zpGbhWtqjK/DPzA3HmNzn5OZFJM4aSGjqUZZLHQp2f5emtWo5HZYxl/AmYRQXkZtV9TfgNxEZD6wD5vkmNGP8XGoq2z96gsq7PgJgY6k7uaXnKG4qbntVmOzJm4RRBPhBRI4AW4BiQIovgjLGr6ly7PsxnPvpAyon72O6dKRsj3/RrGZFpyMz5oq8SRitgU1AQ1z7eyvWujDGK8ePHmbLJ8/Q5OjnAKwtdx89HhhNYKAnC0cb4yxvNlDa4H65wv1jjPGUKvE7f+bEJ4/QJGUP07iNhgPepW6ZG52OzBiP2RrIxvhY3NGTxE+6k6qnVxGowfwYNZL2rftQslBep0MzxiuWMIzxoa1LZ1Lp20fII8lMy9+HsM5P0axGJafDMuaaeJwwRORvwMeq+qcP4zHGL5w5spe4yQ9S/fQqjkox9lZ7iE5dn6Vw/jxOh2bMNfOmhVEKWC0ia4GJwMKr7eNtTK6Tkszmcf2peegLqgCfSjSVer1F/ao3Ox2ZMdfNm0HvF0VkGNABeAAYIyKfAR+q6g5fBWhMTrHxpwUE/fgfap5ZDcCu8Cfo1eM1h6MyJvN4u+OeishB4CCuVWuLAzNF5BtVfcYXARqT3SUf2sZvc9+m1r7pAMyWtrQe+hEVCxd0ODJjMpc3YxhP4NpO9QiubVSfVtUkEQnAtUveFROGiFQFpqcpCgNewrXT3lhcW68mA4+p6qp0194M/A9Xt1gqME5VR3kauzG+kJqqrN64kYoL7qNawm6SNJC93T6nVeXGFLMZUMYPeZQwxLUAfwTQXVX3pH1PVVNFpNPVPkNVY3Eti46IBAL7gdnAeOBVVV3g3rP7LaBVusuTgf9T1bUiUhiIcbdqtngSvzGZLWnPKmK/+ZCGcZ8C8EHRwdzVsy+VytgMKOO/PEoY7q6oOumTRZr3t3p537bADlXdIyKKa9kRcK2AeyCDz/8d+N39+pSIbAXK4lqixJgso6okHool76T2hANrAiMpGtmFnm0GUaxgPqfDM8anvBnDWCEi9VV1dSbctxcwzf36SWChiIwAAoAmV7pQRCoAdXDtxWFM1jl5gD/e7cBNSfsAeD3fUHo98CSVbirqcGDGZA1v9vTegmsNqd1APCC4Gh+1vbqhSB5crYiaqnpIREYDS1R1lojcDQxQ1XaXubYQsAR4XVU/v8w5A4ABAKGhofX27MmwUWSMV377agS3xLxGogayKagGUrsnkV3+ZtulGr9zpT29vUkY5TMqv1w31RU+pyswSFU7uI9PAMXc3V6Ca2e/S/aiFJFgYC6u5z/e9uReUVFRumbNGm/CM+YiibGL2TX/HSofX85p8jOl2GP0f+x5Cua1RRKMf7pSwvDmb33fy5QP9zKe3vzVHQWu1kZLXLOl2uCacXURdyL5ENjqabIw5nokJ6cwd8Kr3H7wPaqSzDEpypl7vuLxyhFOh2aMY7xJGPFpXucDOgFeDXaLSAGgPfBImuKHgVEiEgQk4O5OEpEywARVjQaaAvcBG0Vkvfu6v6vqfG/ub4wnzu1YzuaFE+j2h6vX88PK7/HgPfdyg8NxGeM0j7ukLrlQJC/wparemrkhZR7rkjLeOHgigZ9XLqfTT3cRRArr8tSj7MDPCSle1MYqTK6RWV1S6RXA9fCdMTmexh/h6Lud6Za8hVOan9UNRtKkXXfy5bXFAo05z5snvTfi2mUPIBAIwfvxC2OynZiZ/6bepn9QTYUfQnpTud2DtKlW3+mwjMl2vGlhpH2aOxk4pKrJmRyPMVnmzMavSJ49iHqpJ1inlTnc8Hnad+xu3U/GXIY3q9XaAw3GLyQkJhEzeyT1tv2bInqOTcXbUu3BSdQpZA/gGXMlHu88LyJTRKRYmuPiIjLRJ1EZ4yNnlozm7BuVabr1H6xPqcTc1gsIH/w5+S1ZGHNV3nRJ1VbV4+cPVPVPEamT+SEZk/kSEpPZMvtN6m59iwLAxogXqR39JAXyBjsdmjE5hjcJI0BEip/folVEbvDyemOynirxaz9D5w6hrsaziIbccO8Uom4p7XRkxuQ43vzC/w/wk4jMxDVb6m7gdZ9EZUwm0NRUtv7vSWrsnsIRLcrmCgOo3ulFbg6x7idjroU3+2F8D6zBtXyH4Nobw5YXN9mPKge+GU3CqsnUSN7JV3k6UqLHKJpUucnpyIzJ0bzZD2OOqtbD9qAw2djJw3HsmvwwEfE/kUgQW27uRYf7x5A32MYqjLle3nRJ/ZyJ+2EYk7mSz7F7yUfcvOxpqmgQ/0q9j+6PvU6N0tb9ZExm8SZhtAYGishurmM/DGMy2y+/rCVowVBqJqwDYHeTf/Fs+wcJCLAH8IzJTN4kjI4+i8KYa5GUwL71i6g6tz/5JIlPg7vR/JF3qF7S1pU1xhe8SRh7gXuAMFUdLiKhQCnAngA3We70yWMkv9uQm5P+IEGC2RHxNB3bD6VooQJOh2aM3/ImYfwXSMU1S2o4cAqYBdgqbSbrqLL1w4epHjcDgI8KPUB41yHUqZzhhpDGmEzkTcJoqKp1RWQdXHjS29Z+Nlnmj3ULCJw3mOrJh9gvpUiu1pX7er7ldFjG5BreJIwkEQnEvcS5iITganF4RESqAtPTFIUBL+HamnUsrl38koHHVHVVBtffBozCtbT6BFV9w4vYTQ6mR7azb/KDhJ7+BYDpgZ2JfOg9qtoMKGOylDcJYzQwG7hJRF4HegDDPL1YVWOBSAB34tnv/rzxwKuqukBEooG3gFZpr3Wf/x6u7V3jgNUi8qU9OOj/li/4hCprXiU05SAJGsze2ybSs3EXp8MyJlfyZnnzj0UkBmjrLuqqqtuu8b5tgR2qukdEFCjiLi8KHMjg/AbAb6q6E0BEPgW6Yg8R+q0zG74g7vsPafrnEgDGBd/LA8+MpIo9gGeMY66aMETky/RF7j9vFRFU9Vr+udcLmOZ+/SSwUERG4FpuvUkG55cF9qU5jgMaXibeAcAAgNDQ0GsIzTgpMTmVjVs2U+vzB6lCEqukNkW7vkm3sLoEW7IwxlGetDAa4/plPQ1YyV8J45q4B8q7AM+7ix4FhqjqLBG5G/gQaJf+sgw+SjMoQ1XHAeMAoqKiMjzHZFMH1rF/+nPUO/EzqSq8WWkyD3aPpmShvE5HZozBs4RRCtfYQW+gDzAPmKaqm6/xnh2Btap6yH3cFxjsfj0DmJDBNXHAzWmOy5Fx15XJofYu/YjQ7/5GReDr4LYUbfIgT7WMJtCe1jYm27hqwlDVFOBr4GsRyYsrcfwgIsNV9d1ruGdv/uqOAtcv/pa4Zku1AbZncM1qoLKIVMQ1WN4LV/IyOVzqoW1sn/ggZRJ2cEbyMr3si3Tu/Yi1KozJhjxd3jwvcDuuX/YVcM2Y+tzbm4lIAVytlUfSFD8MjBKRICAB9/iDiJTBNX02WlWTReRvwEJc02onXkcLx2QDqsqiD56h1cHJVCWRHYEVCIx+iweibnU6NGPMZYjqlbv5RWQKEA4sAD5V1U1ZEVhmiIqK0jVr1jgdhknn1C9fceibd7jldAzHtBDzw17knvsH4tp2xRjjJBGJUdWojN7zpIVxH67VaasAT6T5n/r8arVFLnehMWmdOJPE3Pmz6b7pb9zCObYG1yRs6ELuzV/Y6dCMMR7wZAwjICsCMf4taccyfp0zkntOLeakFmBV9GJaNrRlyIzJSbx50tsYryUkJrHxi1HU3/wa9YFlN/ahZvfnaVnKnpExJqexhGF85tyuFciUrtTnHOsCwjnR/h1aNKxvGxsZk0NZwjCZLiU5idWfvErdne+Th2Q2ht5LtV5vkb9AQadDM8ZcB0sYJlOdWPIeQT/8g0Z6hm+0Pufa/oNOLRo5HZYxJhNYwjCZIiHhLCs/+zctd/4HgLU1nqP9Xc+BTZU1xm9YwjDXR5VjP3+ELnqFlnqUb1PrEtJ/OnUr3Oh0ZMaYTGYJw1yzpJRUdkx8mGr7Z7JLS7O59htUbt6H0BttYyNj/JElDHNNjq34iLOLXqOaHuLzfN2p1W8ULUrZM5zG+DNLGMYrp/88SNzYO6l2bhO/awl2hvUhuuc75MuXz+nQjDE+ZgnDeCYpgW3TnqPazklUVmFUnv407f0cUWE3OR2ZMSaLWMIwV7Vl7TJ04YvUPLeeJALZVXMQg+9+zemwjDFZzBKGuayUP37l1x9nUfaXdykq8cwLbEOL/5tGlQLW/WRMbmQJw2To183rKPpVP6on7OYkBdhx21Ta1LmN/HltX21jcitLGOZiCSfZueg9qqx9A4DRhYfQpc8gKpUOcTgwY4zTsixhiEhVYHqaojDgJaAxUNVdVgw4rqqRGVw/BHgIUGAj8ICqJvgw5FznaOxyik7rTBgpHJQQgho9wuPt/w8JsBXujTFZmDBUNRaIBBCRQFx7c89W1ZHnzxGR/wAn0l8rImWBJ4AaqnpWRD7Dta/3ZJ8HnhucOcb+d2+l7NlfAfiw0ADCuw6lYeXSDgdmjMlOnOqSagvsUNU95wvEtZXf3UCby1wTBOQXkSSgAHDA51HmAqvmTaL8un9TNnk/mwKqkq/1UzzY/G6nwzLGZENOJYxewLR0Zc2BQ6q6Pf3JqrpfREYAe4GzwCJVXZTRB4vIAGAAQGiobdJzOWc2zeOPr16lbsJ2giSVKdKV6CHjCSmc1+nQjDHZVJZ3TotIHqALMCPdW725NImcv6Y40BWoCJQBCorIvRmdq6rjVDVKVaNCQmygNj1VZf6CL5FZD1LhXCy/SBXi+q6i1wuTLFkYY67IiRZGR2Ctqh46XyAiQUB3oN5lrmkH7FLVw+7zPweaAFN9HKtfORkzk0PLJhN9fBkAk2/+B337/w2xJciNMR5wImFk1JJoB2xT1bjLXLMXaCQiBXB1SbUF1vguRP9y+lwyW2N+oP6iBykC/BTcmNK3P8/9tVtYsjDGeCxLE4b7F3574JF0b10ypiEiZYAJqhqtqitFZCawFkgG1gHjsiDkHE9P/8HRUbdSP2kn8ZqPbxuMp2Xr2yhaII/ToRljchhRVadj8JmoqChdsyb3NkTWfzGayHXDAFhWqCPlW/cltF5Hh6MyxmRnIhKjqlEZvWdPevuhlN0/cfqT/kQm/s4BQlhT7Wluv3sAgQHW/WSMuXaWMPxIcnIKi6e9Q5OdIymqp9iQvyHl+0+kS0g5p0MzxvgBSxh+IinmI/Yu+YjbTq5kR2ppvm32P+7ocLlnII0xxnuWMHK4E/HnWD/9VVrufY9KwLcl76P1o6OoFBjodGjGGD9jCSMHO7t7NQFTutFST7M4oAl57p5I22q2/pMxxjcsYeRACYnJbJz1L+rGvo1oIKurP039zk9TtKA9qW2M8R1LGDnM4WWTKPLt09QniW+0AWc7jqJLoxpOh2WMyQUsYeQQJ06eYNWnb9Bm//uclMLsjRhKu67PIgE2VmGMyRqWMLK71BT2LXibkNVv0Z5ElmltCt03jTq32FRZY0zWsoSRje0/cow/vnyVOnsnc0SLsrvBK1Rv0Z+ShfM7HZoxJheyhJEdpaawb/knlFw8lLKSyMI87Yl64mOqFbJBbWOMcyxhZDNJCfEceacpN5/bRZzciEbeS8v2Q8hnM6CMMQ6zhJFdqPLbxIe5Zd8MSgMzivajfKfnaGD7ahtjsglLGNnA/s3LCZg9gFuS44gnP4dvuYu77h3ldFjGGHMRSxgOSjm4hX0zn6fCkR8A+JQO1B04gSqlijobmDHGZMAShkNWL55J5eVDqaAnOKfB7GjzPnc370GALUFujMmmsixhiEhVYHqaojDgJaAxUNVdVgw4rqqRGVxfDJgAhAMK9FfVFb6L2Dd032r2rf6K+htcXU7jCj1K38GvUSM42OHIjDHmyrIsYahqLBAJICKBwH5gtqqOPH+OiPwHOHGZjxgFfK2qPUQkD1DApwFnsoTEJJavWE7jH3oTqmeI0aqU7DmGPmF1yGvJwhiTAzjVJdUW2KGqe84XiIgAdwOXbOIgIkWAFkA/AFVNBBKzJNLMoMqO8ffT9vB8AEbc9BbN299B+VtudDgwY4zxnFMJoxcwLV1Zc+CQqm7P4Pww4DAwSUQigBhgsKrGpz9RRAYAAwBCQ0MzNWhvpaYq6374nMilD1GTVGLyRFGu/SCeqt/d0biMMeZaBGT1Dd3dSV2AGene6s2lSeS8IKAu8L6q1gHigecyOlFVx6lqlKpGhYSEZFLU3kuMW8/x1ypQb2l/EsjDjAK9uOHhOdxkycIYk0M50cLoCKxV1UPnC0QkCOgO1LvMNXFAnKqudB/P5DIJw2mqysyp79N81yhK6XG2BtegcOd/cVftVk6HZowx18WJhJFRS6IdsE1V4zK6QFUPisg+EanqHjxvC2zxcZxeO7X6E/Z//yF3xscQIMq0G4fQ69GXcQ3PGGNMzpalCUNECgDtgUfSvXXJmIaIlAEmqGq0u+hx4GN3l9ZO4AEfh+uxpOQU5n2/lNuXP041ktkQHE7VofPpXcAewDPG+I8sTRiqegYokUF5vwzKDgDRaY7XA1E+DO+a6M4lHJn5NN3OxHJOg1jQcg63tWqJBGT58JAxxviUPel9jTQ5kV+Xz6Hq9w9TGviu2J1EdH2cjhXrOB2aMcb4hCWMa3Bu508E/q8zVUnmKMX4ufFYOrS7leBAa1UYY/yXJQxvqLJs0vPU3zuBIJJZH9KZst3/ye2lnX3ewxhjsoIlDA8lrv4fB74fR/MzG/kuJZLk2/5Nh6YNnA7LGGOyjCWMqzh28jQ/ThlGl6MTKaOBLAnpRfOB7xEcZF+dMSZ3sd96V3Bu1wrOfvwQXZLjWChNKNBzPC2rl3M6LGOMcYQljAwcO7iH3+a/S729Eyik+Vle7x2a33ofBfLaqrLGmNzLEkY6fx47StDYxjQgnu/ztCLpthF0qFvZ6bCMMcZxljDSKV60KMurPkGFAgm06vQCEmitCmOMAUsYlwoMomnvbLmuoTHGOMqeNDPGGOMRSxjGGGM8YgnDGGOMRyxhGGOM8YglDGOMMR6xhGGMMcYjljCMMcZ4xBKGMcYYj4iqOh2Dz4jIYWDPNV5eEjiSieHkBFZn/5fb6gtWZ2+VV9WQjN7w64RxPURkjapmuz3Efcnq7P9yW33B6pyZrEvKGGOMRyxhGGOM8YgljMsb53QADrA6+7/cVl+wOmcaG8MwxhjjEWthGGOM8YglDGOMMR6xhJGOiNwmIrEi8puI+M1OSiJys4h8LyJbRWSziAx2l98gIt+IyHb3n8XTXPO8+3uIFZFbnYv+2olIoIisE5G57mO/ri+AiBQTkZkiss3937uxP9dbRIa4/05vEpFpIpLPH+srIhNF5A8R2ZSmzOt6ikg9Ednofm+0iIjHQaiq/bh/gEBgBxAG5AF+AWo4HVcm1a00UNf9ujDwK1ADeAt4zl3+HPCm+3UNd/3zAhXd30ug0/W4hnoPBT4B5rqP/bq+7rpMAR5yv84DFPPXegNlgV1AfvfxZ0A/f6wv0AKoC2xKU+Z1PYFVQGNAgAVAR09jsBbGxRoAv6nqTlVNBD4FujocU6ZQ1d9Vda379SlgK67/2bri+gWD+89u7tddgU9V9Zyq7gJ+w/X95BgiUg64HZiQpthv6wsgIkVw/WL5EEBVE1X1OP5d7yAgv4gEAQWAA/hhfVV1KXAsXbFX9RSR0kARVV2hruzxvzTXXJUljIuVBfalOY5zl/kVEakA1AFWAjep6u/gSirAje7T/OG7GAk8A6SmKfPn+oKrdXwYmOTuipsgIgXx03qr6n5gBLAX+B04oaqL8NP6ZsDbepZ1v05f7hFLGBfLqC/Pr+Ydi0ghYBbwpKqevNKpGZTlmO9CRDoBf6hqjKeXZFCWY+qbRhCubov3VbUOEI+rq+JycnS93X32XXF1u5QBCorIvVe6JIOyHFNfL1yuntdVf0sYF4sDbk5zXA5X89YviEgwrmTxsap+7i4+5G6m4v7zD3d5Tv8umgJdRGQ3rq7FNiIyFf+t73lxQJyqrnQfz8SVQPy13u2AXap6WFWTgM+BJvhvfdPztp5x7tfpyz1iCeNiq4HKIlJRRPIAvYAvHY4pU7hnQnwIbFXVt9O89SXQ1/26L/BFmvJeIpJXRCoClXENluUIqvq8qpZT1Qq4/jt+p6r34qf1PU9VDwL7RKSqu6gtsAX/rfdeoJGIFHD/HW+La3zOX+ubnlf1dHdbnRKRRu7v6/4011yd0yP/2e0HiMY1g2gH8ILT8WRivZrhanpuANa7f6KBEsC3wHb3nzekueYF9/cQixczKbLbD9CKv2ZJ5Yb6RgJr3P+t5wDF/bnewKvANmAT8BGumUF+V19gGq5xmiRcLYUHr6WeQJT7u9oBjMG94ocnP7Y0iDHGGI9Yl5QxxhiPWMIwxhjjEUsYxhhjPGIJwxhjjEcsYRhjjPGIJQxjrkJESojIevfPQRHZ7359WkT+66N7Piki91/lnE9FpLIv7m9MRmxarTFeEJFXgNOqOsKH9wgC1uJaXTj5Cue1BO5V1Yd9FYsxaVkLw5hrJCKt0uyz8YqITBGRRSKyW0S6i8hb7n0HvnYvy3J+L4IlIhIjIgvPL+uQThtgraomi0glEVmb5p6VReT8+ljLgHbuBGOMz1nCMCbzVMK1nHpXYCrwvarWAs4Ct7uTxrtAD1WtB0wEXs/gc5oCMQCqugM4ISKR7vceACa730vFtWx1hI/qY8xF7F8mxmSeBaqaJCIbcW3G9bW7fCNQAagKhAPfuDc5C8S11EN6pXGth3TeBOABERkK9OTi/Rv+wLVKq6er8hpzzSxhGJN5zoHrX/4ikqR/DRCm4vp/TYDNqtr4Kp9zFsiX5ngW8DLwHRCjqkfTvJfPfb4xPmddUsZknVggREQag2u5eRGpmcF5W4Fbzh+oagKwEHgfmJTu3CrAZt+Ea8zFLGEYk0XUte1vD+BNEfkF14rBTTI4dQGubVbT+hjXasOLzheIyE3AWXXvuGaMr9m0WmOyIRGZDTyjqtvdx08BRVV1WJpzhgAnVfVDh8I0uYyNYRiTPT2Ha/B7uzt5VMI13Tat47j2fzAmS1gLwxhjjEdsDMMYY4xHLGEYY4zxiCUMY4wxHrGEYYwxxiOWMIwxxnjk/wFAaj+Z8eXNpAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -139,18 +111,15 @@ "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)}')" + "plt.savefig(f'helio_gr_test.png', dpi=300,facecolor='white', transparent=False)" ] }, {